extract_geometry_tiles.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import re\n",
    "import requests\n",
    "import mercantile\n",
    "from rasterio.io import MemoryFile\n",
    "from rasterio.features import sieve\n",
    "from rasterio import features\n",
    "from shapely.geometry import shape, mapping, Polygon\n",
    "from shapely.ops import cascaded_union\n",
    "from shapely.ops import transform\n",
    "import numpy as np\n",
    "from affine import Affine\n",
    "from pyproj import Transformer\n",
    "import json\n",
    "import matplotlib.pyplot as plt\n",
    "from descartes.patch import PolygonPatch\n",
    "from scipy.ndimage.filters import gaussian_filter\n",
    "from ipywidgets import interact, interactive, FloatSlider, interact_manual, HTML\n",
    "import warnings\n",
    "warnings.simplefilter(\"ignore\")\n",
    "\n",
    "\n",
    "style = \"\"\"\n",
    "    <style>\n",
    "        .output_scroll {\n",
    "            height: unset !important;\n",
    "            border-radius: unset !important;\n",
    "            -webkit-box-shadow: unset !important;\n",
    "            box-shadow: unset !important;\n",
    "        }\n",
    "    </style>\n",
    "    \"\"\"\n",
    "display(HTML(style))\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Enter URL: E.g. https://mapproxy.osm.ch/tiles/AGIS2019/EPSG900913/{zoom}/{x}/{y}.png?origin=nw\")\n",
    "print(\"\")\n",
    "url = input()\n",
    "print(\"\")\n",
    "print(\"Enter Zoom level: E.g. 12\")\n",
    "print(\"\")\n",
    "zoom = int(input())\n",
    "print(\"\")\n",
    "print(\"Enter extent: (Format: min_lon,min_lat,max_lon,max_lat) E.g. 7.66646,47.08713,8.45118,47.63340\")\n",
    "print(\"\")\n",
    "bbox_str = input()\n",
    "\n",
    "min_lon, min_lat, max_lon, max_lat = map(float, bbox_str.split(','))\n",
    "\n",
    "print(\"Extent:\")\n",
    "print(\"West: {}\".format(min_lon))\n",
    "print(\"East: {}\".format(max_lon))\n",
    "print(\"South: {}\".format(min_lat))\n",
    "print(\"North: {}\".format(max_lat))\n",
    "\n",
    "parameters = {}\n",
    "# {z} instead of {zoom}\n",
    "if '{z}' in url:\n",
    "    print('{z} found instead of {zoom} in tile url')\n",
    "    exit()\n",
    "\n",
    "if '{apikey}' in url:\n",
    "    print(\"Not possible to check URL, apikey is required.\")\n",
    "    exit\n",
    "\n",
    "if \"{switch:\" in url:\n",
    "    match = re.search(r'switch:?([^}]*)', url)\n",
    "    switches = match.group(1).split(',')\n",
    "    tms_url = url.replace(match.group(0), 'switch')\n",
    "    parameters['switch'] = switches[0]\n",
    "\n",
    "\n",
    "def process_tile(tile):\n",
    "    query_url = url\n",
    "    if '{-y}' in url:\n",
    "        y = 2 ** tile.z - 1 - tile.y\n",
    "        query_url = query_url.replace('{-y}', str(y))\n",
    "    elif '{!y}' in url:\n",
    "        y = 2 ** (tile.z - 1) - 1 - tile.y\n",
    "        query_url = query_url.replace('{!y}', str(y))\n",
    "    else:\n",
    "        query_url = query_url.replace('{y}', str(tile.y))\n",
    "    parameters['x'] = tile.x\n",
    "    parameters['zoom'] = tile.z\n",
    "    query_url = query_url.format(**parameters)\n",
    "    print(\"Request tile url:\", query_url)\n",
    "\n",
    "    data = None\n",
    "    for i in range(3):\n",
    "        try:\n",
    "            r = requests.get(query_url)\n",
    "            data = r.content\n",
    "        except Exception as e:\n",
    "            print(str(e))\n",
    "        break\n",
    "    if data is None:\n",
    "        return []\n",
    "\n",
    "    bounds = mercantile.bounds(*tile)\n",
    "\n",
    "    try:\n",
    "        with MemoryFile(data) as memfile:\n",
    "            with memfile.open() as dataset:\n",
    "                pixel_x = (bounds.east - bounds.west) / dataset.width\n",
    "                pixel_y = (bounds.north - bounds.south) / dataset.height\n",
    "                geotransform = (bounds.west, pixel_x, 0, bounds.north, 0, -pixel_y)\n",
    "\n",
    "                data = np.zeros(shape=(dataset.height, dataset.width))\n",
    "                for band in range(dataset.count):\n",
    "                    data += dataset.read(band + 1)\n",
    "\n",
    "                # Convert extrema to 0\n",
    "                data[data < 2] = 0\n",
    "                data[data > 254] = 0\n",
    "\n",
    "                # Blur data to avoid holes\n",
    "                data = gaussian_filter(data, sigma=3)\n",
    "\n",
    "                # Convert other pixels to 1\n",
    "                data[data > 0] = 1\n",
    "\n",
    "                # Filter small areas\n",
    "                data = sieve(data.astype(np.uint8),\n",
    "                             size=int(dataset.height * dataset.width * 0.001),\n",
    "                             connectivity=8)\n",
    "\n",
    "                # Create shapes for all areas with value = 1\n",
    "                shapes = list(features.shapes(data.astype(np.uint8), transform=Affine.from_gdal(*geotransform)))\n",
    "                buffer_dist = max(pixel_x, pixel_y) * 5\n",
    "                geoms = [shape(s[0]).buffer(buffer_dist).buffer(-buffer_dist)\n",
    "                         for s in shapes if int(s[1]) == 1]\n",
    "                return geoms\n",
    "    except Exception as e:\n",
    "        print(str(e))\n",
    "    return []\n",
    "\n",
    "\n",
    "def simplify(geom, distance=200):\n",
    "    \"\"\" Simplify geometry in epsg:3857\"\"\"\n",
    "    transformer = Transformer.from_crs(\"epsg:4326\", \"epsg:3857\")\n",
    "    transformer_back = Transformer.from_crs(\"epsg:3857\", \"epsg:4326\")\n",
    "    geom_3857 = transform(transformer.transform, geom)\n",
    "    geom_3857_simplified = geom_3857.simplify(distance, preserve_topology=False)\n",
    "    geom_simplified = transform(transformer_back.transform, geom_3857_simplified)\n",
    "    return geom_simplified\n",
    "\n",
    "\n",
    "geoms = []\n",
    "tiles = list(mercantile.tiles(west=min_lon, south=min_lat, east=max_lon, north=max_lat, zooms=zoom))\n",
    "for tile in tiles:\n",
    "    tile_geoms = process_tile(tile)\n",
    "    if len(tile_geoms) > 0:\n",
    "        geoms.extend(tile_geoms)\n",
    "geom = cascaded_union(geoms)\n",
    "geom = simplify(geom)\n",
    "\n",
    "\n",
    "def plot_geometry(geom):\n",
    "    bounds = geom.bounds\n",
    "    fig, ax = plt.subplots()\n",
    "    if isinstance(geom, Polygon):\n",
    "        geom = [geom]\n",
    "    for g in geom:\n",
    "        patch = PolygonPatch(g, facecolor='#6699cc', edgecolor='#6699cc',\n",
    "                             alpha=0.5, zorder=2)\n",
    "        ax.add_patch(patch)\n",
    "    ax.set_xlim(bounds[0], bounds[2])\n",
    "    ax.set_ylim(bounds[1], bounds[3])\n",
    "    plt.show()\n",
    "\n",
    "plot_geometry(geom)\n",
    "\n",
    "print(json.dumps(mapping(geom), indent=4))\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}