google-research
352 строки · 17.5 Кб
1{2"cells": [3{4"cell_type": "markdown",5"metadata": {6"id": "iVeYZpy3fF0N"7},8"source": [9"##### Copyright 2021 Google LLC. Licensed under the Apache License, Version 2.0 (the \"License\");"10]11},12{13"cell_type": "code",14"execution_count": null,15"metadata": {16"id": "0eXL156ae-iT"17},18"outputs": [],19"source": [20"# Licensed under the Apache License, Version 2.0 (the \"License\");\n",21"# you may not use this file except in compliance with the License.\n",22"# You may obtain a copy of the License at\n",23"#\n",24"# https://www.apache.org/licenses/LICENSE-2.0\n",25"#\n",26"# Unless required by applicable law or agreed to in writing, software\n",27"# distributed under the License is distributed on an \"AS IS\" BASIS,\n",28"# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.\n",29"# See the License for the specific language governing permissions and\n",30"# limitations under the License"31]32},33{34"cell_type": "markdown",35"metadata": {36"id": "LbGgdE4mj1hd"37},38"source": [39"### Step 1. Prepare a compressed CSV file using [Open Buildings](https://sites.research.google/open-buildings/) data [takes 1-15 minutes depending on the region]"40]41},42{43"cell_type": "code",44"execution_count": null,45"metadata": {46"cellView": "form",47"id": "qP6ADuzRdZTF"48},49"outputs": [],50"source": [51"#@markdown First, select a region from either the [Natural Earth low res](https://www.naturalearthdata.com/downloads/110m-cultural-vectors/110m-admin-0-countries/) (fastest), [Natural Earth high res](https://www.naturalearthdata.com/downloads/10m-cultural-vectors/10m-admin-0-countries/) or [World Bank high res](https://datacatalog.worldbank.org/dataset/world-bank-official-boundaries) shapefiles:\n",52"region_border_source = 'Natural Earth (Low Res 110m)' #@param [\"Natural Earth (Low Res 110m)\", \"Natural Earth (High Res 10m)\", \"World Bank (High Res 10m)\"]\n",53"region = 'GHA (Ghana)' #@param [\"\", \"AGO (Angola)\", \"BDI (Burundi)\", \"BEN (Benin)\", \"BFA (Burkina Faso)\", \"BGD (Bangladesh)\", \"BRN (Brunei)\", \"BTN (Bhutan)\", \"BWA (Botswana)\", \"CAF (Central African Republic)\", \"CIV (C\\u00f4te d'Ivoire)\", \"CMR (Cameroon)\", \"COD (Democratic Republic of the Congo)\", \"COG (Republic of the Congo)\", \"COM (Comoros)\", \"CPV (Cape Verde)\", \"DJI (Djibouti)\", \"DZA (Algeria)\", \"EGY (Egypt)\", \"ERI (Eritrea)\", \"ETH (Ethiopia)\", \"GAB (Gabon)\", \"GHA (Ghana)\", \"GIN (Guinea)\", \"GMB (The Gambia)\", \"GNB (Guinea-Bissau)\", \"GNQ (Equatorial Guinea)\", \"IDN (Indonesia)\", \"IOT (British Indian Ocean Territory)\", \"KEN (Kenya)\", \"KHM (Cambodia)\", \"LAO (Laos)\", \"LBR (Liberia)\", \"LKA (Sri Lanka)\", \"LSO (Lesotho)\", \"MDG (Madagascar)\", \"MDV (Maldives)\", \"MOZ (Mozambique)\", \"MRT (Mauritania)\", \"MUS (Mauritius)\", \"MWI (Malawi)\", \"MYS (Malaysia)\", \"MYT (Mayotte)\", \"NAM (Namibia)\", \"NER (Niger)\", \"NGA (Nigeria)\", \"NPL (Nepal)\", \"PHL (Philippines)\", \"REU (R\\u00e9union)\", \"RWA (Rwanda)\", \"SDN (Sudan)\", \"SEN (Senegal)\", \"SGP (Singapore)\", \"SHN (Saint Helena, Ascension and Tristan da Cunha)\", \"SLE (Sierra Leone)\", \"SOM (Somalia)\", \"STP (S\\u00e3o Tom\\u00e9 and Pr\\u00edncipe)\", \"SWZ (Eswatini)\", \"SYC (Seychelles)\", \"TGO (Togo)\", \"THA (Thailand)\", \"TLS (Timor-Leste)\", \"TUN (Tunisia)\", \"TZA (Tanzania)\", \"UGA (Uganda)\", \"VNM (Vietnam)\", \"ZAF (South Africa)\", \"ZMB (Zambia)\", \"ZWE (Zimbabwe)\"]\n",54"#@markdown **or** specify an area of interest in [WKT format](https://en.wikipedia.org/wiki/Well-known_text_representation_of_geometry) (assumes crs='EPSG:4326'); this [tool](https://arthur-e.github.io/Wicket/sandbox-gmaps3.html) might be useful.\n",55"your_own_wkt_polygon = '' #@param {type:\"string\"}\n",56"#@markdown Select type of data to download here:\n",57"data_type = 'polygons' #@param [\"polygons\", \"points\"]\n",58"\n",59"!sudo apt-get install swig\n",60"!pip install s2geometry pygeos geopandas\n",61"\n",62"import functools\n",63"import glob\n",64"import gzip\n",65"import multiprocessing\n",66"import os\n",67"import shutil\n",68"import tempfile\n",69"from typing import List, Optional, Tuple\n",70"\n",71"import geopandas as gpd\n",72"from IPython import display\n",73"import pandas as pd\n",74"import s2geometry as s2\n",75"import shapely\n",76"import tensorflow as tf\n",77"import tqdm.notebook\n",78"\n",79"BUILDING_DOWNLOAD_PATH = ('gs://open-buildings-data/v3/'\n",80" f'{data_type}_s2_level_6_gzip_no_header')\n",81"\n",82"def get_filename_and_region_dataframe(\n",83" region_border_source: str, region: str,\n",84" your_own_wkt_polygon: str) -\u003e Tuple[str, gpd.geodataframe.GeoDataFrame]:\n",85" \"\"\"Returns output filename and a geopandas dataframe with one region row.\"\"\"\n",86"\n",87" if your_own_wkt_polygon:\n",88" filename = f'open_buildings_v3_{data_type}_your_own_wkt_polygon.csv.gz'\n",89" region_df = gpd.GeoDataFrame(\n",90" geometry=gpd.GeoSeries.from_wkt([your_own_wkt_polygon]),\n",91" crs='EPSG:4326')\n",92" if not isinstance(region_df.iloc[0].geometry,\n",93" shapely.geometry.polygon.Polygon) and not isinstance(\n",94" region_df.iloc[0].geometry,\n",95" shapely.geometry.multipolygon.MultiPolygon):\n",96" raise ValueError(\"`your_own_wkt_polygon` must be a POLYGON or \"\n",97" \"MULTIPOLYGON.\")\n",98" print(f'Preparing your_own_wkt_polygon.')\n",99" return filename, region_df\n",100"\n",101" if not region:\n",102" raise ValueError('Please select a region or set your_own_wkt_polygon.')\n",103"\n",104" if region_border_source == 'Natural Earth (Low Res 110m)':\n",105" url = ('https://www.naturalearthdata.com/http//www.naturalearthdata.com/'\n",106" 'download/110m/cultural/ne_110m_admin_0_countries.zip')\n",107" !wget -N {url}\n",108" display.clear_output()\n",109" region_shapefile_path = os.path.basename(url)\n",110" source_name = 'ne_110m'\n",111" elif region_border_source == 'Natural Earth (High Res 10m)':\n",112" url = ('https://www.naturalearthdata.com/http//www.naturalearthdata.com/'\n",113" 'download/10m/cultural/ne_10m_admin_0_countries.zip')\n",114" !wget -N {url}\n",115" display.clear_output()\n",116" region_shapefile_path = os.path.basename(url)\n",117" source_name = 'ne_10m'\n",118" elif region_border_source == 'World Bank (High Res 10m)':\n",119" url = ('https://development-data-hub-s3-public.s3.amazonaws.com/ddhfiles/'\n",120" '779551/wb_countries_admin0_10m.zip')\n",121" !wget -N {url}\n",122" !unzip -o {os.path.basename(url)}\n",123" display.clear_output()\n",124" region_shapefile_path = 'WB_countries_Admin0_10m'\n",125" source_name = 'wb_10m'\n",126"\n",127" region_iso_a3 = region.split(' ')[0]\n",128" filename = (f'open_buildings_v3_{data_type}_'\n",129" f'{source_name}_{region_iso_a3}.csv.gz')\n",130" region_df = gpd.read_file(region_shapefile_path).query(\n",131" f'ISO_A3 == \"{region_iso_a3}\"').dissolve(by='ISO_A3')[['geometry']]\n",132" print(f'Preparing {region} from {region_border_source}.')\n",133" return filename, region_df\n",134"\n",135"\n",136"def get_bounding_box_s2_covering_tokens(\n",137" region_geometry: shapely.geometry.base.BaseGeometry) -\u003e List[str]:\n",138" region_bounds = region_geometry.bounds\n",139" s2_lat_lng_rect = s2.S2LatLngRect_FromPointPair(\n",140" s2.S2LatLng_FromDegrees(region_bounds[1], region_bounds[0]),\n",141" s2.S2LatLng_FromDegrees(region_bounds[3], region_bounds[2]))\n",142" coverer = s2.S2RegionCoverer()\n",143" # NOTE: Should be kept in-sync with s2 level in BUILDING_DOWNLOAD_PATH.\n",144" coverer.set_fixed_level(6)\n",145" coverer.set_max_cells(1000000)\n",146" return [cell.ToToken() for cell in coverer.GetCovering(s2_lat_lng_rect)]\n",147"\n",148"\n",149"def s2_token_to_shapely_polygon(\n",150" s2_token: str) -\u003e shapely.geometry.polygon.Polygon:\n",151" s2_cell = s2.S2Cell(s2.S2CellId_FromToken(s2_token, len(s2_token)))\n",152" coords = []\n",153" for i in range(4):\n",154" s2_lat_lng = s2.S2LatLng(s2_cell.GetVertex(i))\n",155" coords.append((s2_lat_lng.lng().degrees(), s2_lat_lng.lat().degrees()))\n",156" return shapely.geometry.Polygon(coords)\n",157"\n",158"\n",159"def download_s2_token(\n",160" s2_token: str, region_df: gpd.geodataframe.GeoDataFrame) -\u003e Optional[str]:\n",161" \"\"\"Downloads the matching CSV file with polygons for the `s2_token`.\n",162"\n",163" NOTE: Only polygons inside the region are kept.\n",164" NOTE: Passing output via a temporary file to reduce memory usage.\n",165"\n",166" Args:\n",167" s2_token: S2 token for which to download the CSV file with building\n",168" polygons. The S2 token should be at the same level as the files in\n",169" BUILDING_DOWNLOAD_PATH.\n",170" region_df: A geopandas dataframe with only one row that contains the region\n",171" for which to keep polygons.\n",172"\n",173" Returns:\n",174" Either filepath which contains a gzipped CSV without header for the\n",175" `s2_token` subfiltered to only contain building polygons inside the region\n",176" or None which means that there were no polygons inside the region for this\n",177" `s2_token`.\n",178" \"\"\"\n",179" s2_cell_geometry = s2_token_to_shapely_polygon(s2_token)\n",180" region_geometry = region_df.iloc[0].geometry\n",181" prepared_region_geometry = shapely.prepared.prep(region_geometry)\n",182" # If the s2 cell doesn't intersect the country geometry at all then we can\n",183" # know that all rows would be dropped so instead we can just return early.\n",184" if not prepared_region_geometry.intersects(s2_cell_geometry):\n",185" return None\n",186" try:\n",187" # Using tf.io.gfile.GFile gives better performance than passing the GCS path\n",188" # directly to pd.read_csv.\n",189" with tf.io.gfile.GFile(\n",190" os.path.join(BUILDING_DOWNLOAD_PATH, f'{s2_token}_buildings.csv.gz'),\n",191" 'rb') as gf:\n",192" # If the s2 cell is fully covered by country geometry then can skip\n",193" # filtering as we need all rows.\n",194" if prepared_region_geometry.covers(s2_cell_geometry):\n",195" with tempfile.NamedTemporaryFile(mode='w+b', delete=False) as tmp_f:\n",196" shutil.copyfileobj(gf, tmp_f)\n",197" return tmp_f.name\n",198" # Else take the slow path.\n",199" # NOTE: We read in chunks to save memory.\n",200" csv_chunks = pd.read_csv(\n",201" gf, chunksize=2000000, dtype=object, compression='gzip', header=None)\n",202" tmp_f = tempfile.NamedTemporaryFile(mode='w+b', delete=False)\n",203" tmp_f.close()\n",204" for csv_chunk in csv_chunks:\n",205" points = gpd.GeoDataFrame(\n",206" geometry=gpd.points_from_xy(csv_chunk[1], csv_chunk[0]),\n",207" crs='EPSG:4326')\n",208" # sjoin 'within' was faster than using shapely's 'within' directly.\n",209" points = gpd.sjoin(points, region_df, predicate='within')\n",210" csv_chunk = csv_chunk.iloc[points.index]\n",211" csv_chunk.to_csv(\n",212" tmp_f.name,\n",213" mode='ab',\n",214" index=False,\n",215" header=False,\n",216" compression={\n",217" 'method': 'gzip',\n",218" 'compresslevel': 1\n",219" })\n",220" return tmp_f.name\n",221" except tf.errors.NotFoundError:\n",222" return None\n",223"\n",224"# Clear output after pip install.\n",225"display.clear_output()\n",226"filename, region_df = get_filename_and_region_dataframe(region_border_source,\n",227" region,\n",228" your_own_wkt_polygon)\n",229"# Remove any old outputs to not run out of disk.\n",230"for f in glob.glob('/tmp/open_buildings_*'):\n",231" os.remove(f)\n",232"# Write header to the compressed CSV file.\n",233"with gzip.open(f'/tmp/{filename}', 'wt') as merged:\n",234" if data_type == \"polygons\":\n",235" merged.write(','.join([\n",236" 'latitude', 'longitude', 'area_in_meters', 'confidence', 'geometry',\n",237" 'full_plus_code'\n",238" ]) + '\\n')\n",239" else:\n",240" merged.write(','.join([\n",241" 'latitude', 'longitude', 'area_in_meters', 'confidence', \n",242" 'full_plus_code'\n",243" ]) + '\\n')\n",244"download_s2_token_fn = functools.partial(download_s2_token, region_df=region_df)\n",245"s2_tokens = get_bounding_box_s2_covering_tokens(region_df.iloc[0].geometry)\n",246"# Downloads CSV files for relevant S2 tokens and after filtering appends them\n",247"# to the compressed output CSV file. Relies on the fact that concatenating\n",248"# gzipped files produces a valid gzip file.\n",249"# NOTE: Uses a pool to speed up output preparation.\n",250"with open(f'/tmp/{filename}', 'ab') as merged:\n",251" with multiprocessing.Pool(4) as e:\n",252" for fname in tqdm.notebook.tqdm(\n",253" e.imap_unordered(download_s2_token_fn, s2_tokens),\n",254" total=len(s2_tokens)):\n",255" if fname:\n",256" with open(fname, 'rb') as tmp_f:\n",257" shutil.copyfileobj(tmp_f, merged)\n",258" os.unlink(fname)"259]260},261{262"cell_type": "markdown",263"metadata": {264"id": "fwxfj3B1qUWu"265},266"source": [267"### Step 2. Access the generated CSV\n",268"\n",269"We currently support two options:\n",270"\n",271"* Upload to your Google Drive (requires authentication) and then download manually from Google Drive [recommended, very fast]\n",272"* Download directly from the Colab [not recommended, very slow]\n",273"\n",274"When asked to do the authentication please follow the url, paste the generated token into the text field and press enter.\n",275"\n"276]277},278{279"cell_type": "code",280"execution_count": null,281"metadata": {282"id": "6p5iY0WexcJC"283},284"outputs": [],285"source": [286"#@title Upload result to Google Drive [fast, requires authentication]\n",287"\n",288"from google.colab import auth\n",289"from googleapiclient.http import MediaFileUpload\n",290"from googleapiclient.discovery import build\n",291"\n",292"auth.authenticate_user()\n",293"\n",294"drive_service = build('drive', 'v3')\n",295"file_metadata = {\n",296" 'name': filename,\n",297" 'mimeType': 'application/gzip'\n",298"}\n",299"media = MediaFileUpload(f'/tmp/{filename}',\n",300" mimetype='application/gzip',\n",301" resumable=True)\n",302"created = drive_service.files().create(body=file_metadata,\n",303" media_body=media,\n",304" fields='id').execute()\n",305"print('Upload to Google Drive done.')\n",306"print(f'Please download {file_metadata[\"name\"]} manually from Google Drive.')\n",307"print(f'Search link: https://drive.google.com/drive/search?q={filename}')"308]309},310{311"cell_type": "code",312"execution_count": null,313"metadata": {314"id": "pm-hsiYIt7x2"315},316"outputs": [],317"source": [318"#@title Download result directly [slow]\n",319"\n",320"from google.colab import files\n",321"\n",322"files.download(f'/tmp/{filename}')"323]324}325],326"metadata": {327"colab": {328"collapsed_sections": [],329"last_runtime": {330"build_target": "//learning/deepmind/public/tools/ml_python:ml_notebook",331"kind": "private"332},333"name": "Open Buildings - download region polygons or points.",334"private_outputs": true,335"provenance": [336{337"file_id": "1fK7dY741MBONdcfVSLnzEh8vJlBaT89f",338"timestamp": 1626347718873339}340]341},342"kernelspec": {343"display_name": "Python 3",344"name": "python3"345},346"language_info": {347"name": "python"348}349},350"nbformat": 4,351"nbformat_minor": 0352}353