examples/example_building_footprints.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import geopandas as gpd\n",
    "import shapely.geometry\n",
    "import mercantile\n",
    "from tqdm import tqdm\n",
    "import os\n",
    "import tempfile\n",
    "import fiona"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 1 - Define our area of interest (AOI)\n",
    "\n",
    "We define our area of interest (or AOI) as a GeoJSON geometry, then use the `shapely` library to get the bounding box.\n",
    "\n",
    "**Note**: the coordinate reference system for the GeoJSON should be \"EPSG:4326\", i.e. in global lat/lon format."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Geometry copied from https://geojson.io\n",
    "aoi_geom = {\n",
    "    \"coordinates\": [\n",
    "        [\n",
    "            [-122.16484503187519, 47.69090474454916],\n",
    "            [-122.16484503187519, 47.6217555345674],\n",
    "            [-122.06529607517405, 47.6217555345674],\n",
    "            [-122.06529607517405, 47.69090474454916],\n",
    "            [-122.16484503187519, 47.69090474454916],\n",
    "        ]\n",
    "    ],\n",
    "    \"type\": \"Polygon\",\n",
    "}\n",
    "aoi_shape = shapely.geometry.shape(aoi_geom)\n",
    "minx, miny, maxx, maxy = aoi_shape.bounds\n",
    "\n",
    "output_fn = \"example_building_footprints.geojson\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 2 - Determine which tiles intersect our AOI"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "quad_keys = set()\n",
    "for tile in list(mercantile.tiles(minx, miny, maxx, maxy, zooms=9)):\n",
    "    quad_keys.add(int(mercantile.quadkey(tile)))\n",
    "quad_keys = list(quad_keys)\n",
    "print(f\"The input area spans {len(quad_keys)} tiles: {quad_keys}\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 3 - Download the building footprints for each tile that intersects our AOI and crop the results\n",
    "\n",
    "This is where most of the magic happens. We download all the building footprints for each tile that intersects our AOI, then only keep the footprints that are _contained_ by our AOI.\n",
    "\n",
    "*Note*: this step might take awhile depending on how many tiles your AOI covers and how many buildings footprints are in those tiles."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv(\n",
    "    \"https://minedbuildings.blob.core.windows.net/global-buildings/dataset-links.csv\"\n",
    ")\n",
    "\n",
    "idx = 0\n",
    "combined_rows = []\n",
    "\n",
    "with tempfile.TemporaryDirectory() as tmpdir:\n",
    "    # Download the GeoJSON files for each tile that intersects the input geometry\n",
    "    tmp_fns = []\n",
    "    for quad_key in tqdm(quad_keys):\n",
    "        rows = df[df[\"QuadKey\"] == quad_key]\n",
    "        if rows.shape[0] == 1:\n",
    "            url = rows.iloc[0][\"Url\"]\n",
    "\n",
    "            df2 = pd.read_json(url, lines=True)\n",
    "            df2[\"geometry\"] = df2[\"geometry\"].apply(shapely.geometry.shape)\n",
    "\n",
    "            gdf = gpd.GeoDataFrame(df2, crs=4326)\n",
    "            fn = os.path.join(tmpdir, f\"{quad_key}.geojson\")\n",
    "            tmp_fns.append(fn)\n",
    "            if not os.path.exists(fn):\n",
    "                gdf.to_file(fn, driver=\"GeoJSON\")\n",
    "        elif rows.shape[0] > 1:\n",
    "            raise ValueError(f\"Multiple rows found for QuadKey: {quad_key}\")\n",
    "        else:\n",
    "            raise ValueError(f\"QuadKey not found in dataset: {quad_key}\")\n",
    "\n",
    "    # Merge the GeoJSON files into a single file\n",
    "    for fn in tmp_fns:\n",
    "        with fiona.open(fn, \"r\") as f:\n",
    "            for row in tqdm(f):\n",
    "                row = dict(row)\n",
    "                shape = shapely.geometry.shape(row[\"geometry\"])\n",
    "\n",
    "                if aoi_shape.contains(shape):\n",
    "                    if \"id\" in row:\n",
    "                        del row[\"id\"]\n",
    "                    row[\"properties\"] = {\"id\": idx}\n",
    "                    idx += 1\n",
    "                    combined_rows.append(row)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 4 - Save the resulting footprints to file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "schema = {\"geometry\": \"Polygon\", \"properties\": {\"id\": \"int\"}}\n",
    "\n",
    "with fiona.open(output_fn, \"w\", driver=\"GeoJSON\", crs=\"EPSG:4326\", schema=schema) as f:\n",
    "    f.writerecords(combined_rows)"
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}