[PYTHON] [Study memo] Build GeoJSON vector tile server with Fast API

Motivation, overview

In GeoDjango, for example, here In the form, it seems that there is a function that can acquire information from the DB table (using PostGIS) and operate it as a vector tile server.

I wondered if I could do something similar with Python's micro web framework (here I use FastAPI), so I did a lot of research and made a note for myself.

Ideally, it should be PostGIS, or if it is NoSQL, MongoDB's Geospatial Query We built a system using around .mongodb.com/manual/geospatial-queries/), and the distribution format is not only GeoJSON but also more practical binary vector tile format (Reference. I would like to support com / vector-tiles / specification /)) etc., but it is difficult to support all of them suddenly, so here

--Tile server endpoint generation using micro web framework --Receive an HTTP request (GET) in the <url> /{z}/{x}/{y}.geojson format and return a response (GeoJSON or 404) in the appropriate format

The main focus is on. On the other hand

--Connection / linkage with DB (PostgreSQL, MongoDB, etc.) --Delivery in binary format such as mvt (pbf) --Performance evaluation after incorporating the above two functions ――What kind of data query is likely to make the machine load and response speed realistic?

This is not possible at this point, and it will be an issue for the future. (In other words, ** the situation is still far from practical use ** ... If possible, we will investigate and add articles in the future)

Practice

demo

Peek 2020-08-30 01-28.gif

Execution environment

I remember that the dependencies and builds of geospatial libraries were unexpectedly troublesome, so here I am building the environment using conda.

The conda virtual environment is, for example

conda_packages.yml


name: fastapi_geojsontileserver
channels:
  - conda-forge
  - defaults
dependencies:
  # core
  - python==3.7.*
  # FastAPI
  - fastapi
  - uvicorn
  - aiofiles
  # for handling GeoSpatial data
  - pandas>=1.1
  - geopandas>=0.8
  - gdal==3.0.*
  - shapely

Prepare a YAML file like

#Create a virtual environment from a YAML file
conda env create -f conda_packages.yml

#Activate the virtual environment
conda activate <Virtual environment name>

To do.

In addition, edit the virtual environment name (name: part) and the library to be installed (dependencies: part) of the YAML file as needed.

Directory structure / files

For example:

.
├── app
│   ├── __init__.py
│   ├── client
│   │   └── index.html
│   └── main.py
└── data
    └── test.geojson

--ʻApp / main.pyimplements tile server processing --In addition,init.pyis an empty file --For the purpose of checking the operation of the tile server, put HTML files for visualization under ʻapp / client / --Here, I mainly use Leaflet --Put appropriate GeoJSON format data in data / test.geojson as test data --Originally, you should prepare a DB etc., but this time we will simply read GeoJSON with GeoPandas and reproduce the query operation etc. in a pseudo manner.

Examples of ʻapp / main.py and ʻapp / client / index.html are detailed below.

app/main.py

main.py


"""
app main

GeoJSON VectorTileLayer Test
"""

import pathlib
import json
import math
from typing import Optional

from fastapi import (
    FastAPI,
    HTTPException,
)
from fastapi.staticfiles import StaticFiles
from fastapi.responses import RedirectResponse, HTMLResponse
import geopandas as gpd
import shapely.geometry

# const
PATH_STATIC = str(pathlib.Path(__file__).resolve().parent / "client")

EXT_DATA_PATH = "./data/"   # TMP

# test data
gdf_testdata = gpd.read_file(EXT_DATA_PATH + "test.geojson")


def create_app():
    """create app"""
    _app = FastAPI()

    # static
    _app.mount(
        "/client",
        StaticFiles(directory=PATH_STATIC, html=True),
        name="client",
    )

    return _app


app = create_app()


@app.get('/', response_class=HTMLResponse)
async def site_root():
    """root"""
    return RedirectResponse("/client")


@app.get("/tiles/test/{z}/{x}/{y}.geojson")
async def test_tile_geojson(
    z: int,
    x: int,
    y: int,
    limit_zmin: Optional[int] = 8,
) -> dict:
    """
    return GeoJSON Tile
    """

    if limit_zmin is not None:
        if z < limit_zmin:
            raise HTTPException(status_code=404, detail="Over limit of zoom")

    # test data
    gdf = gdf_testdata.copy()

    bbox_polygon = tile_bbox_polygon(z, x, y)

    # filtering
    intersections = gdf.geometry.intersection(bbox_polygon)
    gs_filtered = intersections[~intersections.is_empty]    # geoseries
    gdf_filtered = gpd.GeoDataFrame(
        gdf.loc[gs_filtered.index, :].drop(columns=['geometry']),
        geometry=gs_filtered,
    )

    # NO DATA
    if len(gs_filtered) == 0:
        raise HTTPException(status_code=404, detail="No Data")

    # return geojson
    return json.loads(
        gdf_filtered.to_json()
    )


def tile_coord(
    zoom: int,
    xtile: int,
    ytile: int,
) -> (float, float):
    """
    This returns the NW-corner of the square. Use the function
    with xtile+1 and/or ytile+1 to get the other corners.
    With xtile+0.5 & ytile+0.5 it will return the center of the tile.
    http://wiki.openstreetmap.org/wiki/Slippy_map_tilenames#Tile_numbers_to_lon..2Flat._2
    """
    n = 2.0 ** zoom
    lon_deg = xtile / n * 360.0 - 180.0
    lat_rad = math.atan(math.sinh(math.pi * (1 - 2 * ytile / n)))
    lat_deg = math.degrees(lat_rad)
    return (lon_deg, lat_deg)


def tile_bbox_polygon(
    zoom: int,
    xtile: int,
    ytile: int,
) -> shapely.geometry.Polygon:
    """
    create bbox for Tile by using shapely.geometry
    """
    z = zoom
    x = xtile
    y = ytile

    # get bbox
    nw = tile_coord(z, x, y)
    se = tile_coord(z, x+1, y+1)
    bbox = shapely.geometry.Polygon(
        [
            nw, (se[0], nw[1]),
            se, (nw[0], se[1]), nw
        ]
    )
    return bbox

--Test data (test.geojson) is read by specifying the path with geopandas. --Since both test data and geopandas are provisionally prepared and used, the details are passed. --Function: def test_tile_geojson is the main body of the process, and{z} / {x} / {y}(zoom, position) is Path Parameter Receive as -params /) and in the form of Query Parameter (? = <*** after URL if any other required parameters are available You can receive it in the form> ) (Here, the minimum zoom value limit_zmin is set for testing) --Convert the received {z} / {x} / {y} (each integer value) to a rectangular range of latitude and longitude (bbox), and cut out from the test data (test.geojson) in the range of bbox ( intersection) Returns in the form of geojson. If there is no data, 404 is thrown (Reference). --The conversion of each integer value of {z} / {x} / {y} to latitude / longitude coordinates is performed by the helper functions tile_coord and tile_bbox_polygon. --The processing is mostly based on the processing of django-geojson's TiledGeoJSONLayerView ([source] ](Https://github.com/makinacorpus/django-geojson/blob/master/djgeojson/views.py)) --For tile_coord, most of the processing described in openstreetmap.org wiki is as it is. --HTML files for visualization are mounted and distributed as Static Files in Fast API. --Because a CORS error will occur if you do not deliver in the same domain as the tile server part

app/client/index.html

index.html


<!DOCTYPE html>
<html>

<head>
  <meta charset="UTF-8">
  <meta name="viewport" content="initial-scale=1.0, maximum-scale=1.0" />
  <title>Leaflet GeoJSON TileLayer(Polygon) Test</title>

  <link rel="stylesheet" href="https://unpkg.com/[email protected]/dist/leaflet.css" />
  <script src="https://unpkg.com/[email protected]/dist/leaflet.js"></script>
  <script src="https://unpkg.com/[email protected]/leaflet-hash.js"></script>

  <style>
    #map {
      position: absolute;
      top: 0;
      left: 0;
      bottom: 0;
      right: 0;
    }

    .leaflet-control-container::after {
      content: url(https://maps.gsi.go.jp/image/map/crosshairs.png);
      z-index: 1000;
      display: block;
      position: absolute;
      top: 50%;
      left: 50%;
      transform: translate(-50%, -50%);
    }

  </style>
</head>

<body>
  <div id="map"></div>
  <script>
    // Initalize map
    const map = L.map("map", L.extend({
      minZoom: 5,
      zoom: 14,
      maxZoom: 22,
      center: [35.5, 139.5],
    }, L.Hash.parseHash(location.hash)));
    map.zoomControl.setPosition("bottomright");
    L.hash(map);


    // GeoJSON VectorTileLayer
    const tile_geojson_sample = Object.assign(new L.GridLayer({
        attribution: "hoge",
        minZoom: 4,
        maxZoom: 22,
    }), {
        createTile: function(coords) {
            const template = "http://localhost:8000/tiles/test/{z}/{x}/{y}.geojson?limit_zmin=7";
            const div = document.createElement('div');
            div.group = L.layerGroup();
            fetch(L.Util.template(template, coords)).then(a => a.ok ? a.json() : null).then(geojson => {
                if (!div.group) return;
                if (!this._map) return;
                if (!geojson) return;
                div.group.addLayer(L.geoJSON(geojson, {
                    style: () => {
                        return {}
                    }
                }).bindPopup("test"));
                div.group.addTo(this._map);
            });
                return div;
            }
        }).on("tileunload", function(e) {
        if (e.tile.group) {
            if (this._map) this._map.removeLayer(e.tile.group);
            delete e.tile.group;
        }
    });

    // basemap layers
    const osm = L.tileLayer('http://tile.openstreetmap.jp/{z}/{x}/{y}.png', {
        attribution: "<a href='http://osm.org/copyright' target='_blank'>OpenStreetMap</a> contributors",
        // minZoom: 10,
        maxNativeZoom: 18,
        maxZoom: 22,
    });

    const gsi_std = L.tileLayer(
'https://cyberjapandata.gsi.go.jp/xyz/std/{z}/{x}/{y}.png',
        {
            attribution: "<a href='http://portal.cyberjapan.jp/help/termsofuse.html' target='_blank'>Geographical Survey Institute tile (standard map)</a>",
            maxNativeZoom: 18,
            maxZoom: 22,
            opacity:1
        });

    const gsi_pale = L.tileLayer(
        'http://cyberjapandata.gsi.go.jp/xyz/pale/{z}/{x}/{y}.png',
        {
            attribution: "<a href='http://portal.cyberjapan.jp/help/termsofuse.html' target='_blank'>Geographical Survey Institute tile (light color map)</a>",
            maxNativeZoom: 18,
            maxZoom: 22,
        });

    const gsi_ort = L.tileLayer(
'https://cyberjapandata.gsi.go.jp/xyz/ort/{z}/{x}/{y}.jpg',
        {
            attribution: "<a href='http://portal.cyberjapan.jp/help/termsofuse.html' target='_blank'>Geospatial Information Authority of Japan tile (Ortho)</a>",
            maxNativeZoom: 17,
            maxZoom: 22,
            opacity:0.9
        });

    const gsi_blank = L.tileLayer(
'https://cyberjapandata.gsi.go.jp/xyz/blank/{z}/{x}/{y}.png',
        {
            attribution: "<a href='http://portal.cyberjapan.jp/help/termsofuse.html' target='_blank'>Geographical Survey Institute tile (blank map)</a>",
            maxNativeZoom: 14,
            maxZoom: 22,
            opacity:1,
        });


    L.control.scale({
        imperial: false,
        metric: true,
    }).addTo(map);

    const baseLayers ={
        "Geographical Survey Institute tile (standard map)": gsi_std,
        "Geographical Survey Institute tile (light color map)": gsi_pale,
        "Geospatial Information Authority of Japan tile (Ortho)": gsi_ort,
        "Geographical Survey Institute tile (blank map)": gsi_blank,
        'osm': osm.addTo(map),
    };

    const overlays = {"GeoJSON TileLayer(sample)": tile_geojson_sample};

    L.control.layers(baseLayers, overlays, {position:'topright',collapsed:true}).addTo(map);
    const hash = L.hash(map);

  </script>
</body>

</html>

In the const tile_geojson_sample part, specify the URL of the tile server you created and run + the query string and read it. The main thing is to check if the tile server is working properly, and it is not the essential part, so further details are omitted.

To create the above HTML,

I referred to.

Launch and execute

To start the server, for example, type the following command:

uvicorn app.main:app

The API of the created tile server can be confirmed as Swagger (Open API) from, for example, http: // localhost: 8000. (Automatic document generation + You can easily check the operation on the spot)

image.png

If you are launching locally, you can check the HTML for visualization created above by opening http: // localhost: 8000 / client etc .:

image.png

image.png

image.png

The above execution example uses city boundary data (Polygon, MultiPolygon), but it can be executed in almost the same way by using LineString or Point as the data. By the way, the part that is strangely missing is the problem of the contents of the test data. (When performing intersection processing with geopandas, data selection etc. is performed in advance such as removing ʻinvalid` ones)

If you look at the console log of FastAPI (uvicorn) at runtime,

INFO:     127.0.0.1:54910 - "GET /tiles/test/9/451/202.geojson?limit_zmin=7 HTTP/1.1" 200 OK
INFO:     127.0.0.1:54902 - "GET /tiles/test/9/456/202.geojson?limit_zmin=7 HTTP/1.1" 404 Not Found
INFO:     127.0.0.1:54908 - "GET /tiles/test/9/450/199.geojson?limit_zmin=7 HTTP/1.1" 200 OK
INFO:     127.0.0.1:54910 - "GET /tiles/test/9/457/199.geojson?limit_zmin=7 HTTP/1.1" 404 Not Found
INFO:     127.0.0.1:54896 - "GET /tiles/test/9/457/200.geojson?limit_zmin=7 HTTP/1.1" 404 Not Found
INFO:     127.0.0.1:54904 - "GET /tiles/test/9/450/201.geojson?limit_zmin=7 HTTP/1.1" 200 OK
INFO:     127.0.0.1:54906 - "GET /tiles/test/9/457/201.geojson?limit_zmin=7 HTTP/1.1" 404 Not Found
INFO:     127.0.0.1:54908 - "GET /tiles/test/9/457/202.geojson?limit_zmin=7 HTTP/1.1" 404 Not Found
INFO:     127.0.0.1:54904 - "GET /tiles/test/9/451/198.geojson?limit_zmin=7 HTTP/1.1" 200 OK
INFO:     127.0.0.1:54902 - "GET /tiles/test/9/450/202.geojson?limit_zmin=7 HTTP/1.1" 200 OK
INFO:     127.0.0.1:54896 - "GET /tiles/test/9/452/198.geojson?limit_zmin=7 HTTP/1.1" 200 OK
INFO:     127.0.0.1:54910 - "GET /tiles/test/9/450/198.geojson?limit_zmin=7 HTTP/1.1" 200 OK
INFO:     127.0.0.1:54906 - "GET /tiles/test/9/453/198.geojson?limit_zmin=7 HTTP/1.1" 200 OK
INFO:     127.0.0.1:54910 - "GET /tiles/test/9/454/198.geojson?limit_zmin=7 HTTP/1.1" 200 OK
INFO:     127.0.0.1:54902 - "GET /tiles/test/9/449/199.geojson?limit_zmin=7 HTTP/1.1" 404 Not Found
INFO:     127.0.0.1:54904 - "GET /tiles/test/9/449/200.geojson?limit_zmin=7 HTTP/1.1" 200 OK
INFO:     127.0.0.1:54908 - "GET /tiles/test/9/449/198.geojson?limit_zmin=7 HTTP/1.1" 404 Not Found
INFO:     127.0.0.1:54910 - "GET /tiles/test/9/448/200.geojson?limit_zmin=7 HTTP/1.1" 404 Not Found
INFO:     127.0.0.1:54906 - "GET /tiles/test/9/449/201.geojson?limit_zmin=7 HTTP/1.1" 200 OK
INFO:     127.0.0.1:54896 - "GET /tiles/test/9/448/199.geojson?limit_zmin=7 HTTP/1.1" 404 Not Found
INFO:     127.0.0.1:54906 - "GET /tiles/test/7/112/49.geojson?limit_zmin=7 HTTP/1.1" 200 OK
INFO:     127.0.0.1:54908 - "GET /tiles/test/9/448/198.geojson?limit_zmin=7 HTTP/1.1" 404 Not Found
INFO:     127.0.0.1:54902 - "GET /tiles/test/9/455/198.geojson?limit_zmin=7 HTTP/1.1" 200 OK
INFO:     127.0.0.1:54904 - "GET /tiles/test/9/448/201.geojson?limit_zmin=7 HTTP/1.1" 200 OK
INFO:     127.0.0.1:54906 - "GET /tiles/test/7/111/49.geojson?limit_zmin=7 HTTP/1.1" 404 Not Found
INFO:     127.0.0.1:54902 - "GET /tiles/test/7/114/49.geojson?limit_zmin=7 HTTP/1.1" 200 OK
INFO:     127.0.0.1:54896 - "GET /tiles/test/7/113/49.geojson?limit_zmin=7 HTTP/1.1" 200 OK
INFO:     127.0.0.1:54910 - "GET /tiles/test/7/113/48.geojson?limit_zmin=7 HTTP/1.1" 200 OK
INFO:     127.0.0.1:54902 - "GET /tiles/test/7/112/51.geojson?limit_zmin=7 HTTP/1.1" 200 OK

It feels like.

Summary

Most of the data and processing used are provisional, but I was able to build a tile server that works like that.

As mentioned at the beginning, there are actually many issues, but for the time being, it was a good study of map tiles, vector tiles, and Python's micro web framework in general.

reference

Previous research

Geographical Survey Institute Vector Tile Provision Experiment

Overall quite helpful

-Tile specifications -Tile coordinate confirmation page --Since you can see which area is which {z} / {x} / {y} on the map, it is also convenient for API debugging etc.

Vector tile Format information

GeoJSON

protobuf

The extension is .mvt or .pbf. This time it is not implemented, but as a subsequent issue

Other

Different from "Dynamic Vector Tile Server" but helpful

static vector tile

geobuf

Dynamic tile servers can have load and response time issues, so in some cases it's quite possible that this is a more realistic approach than tiled.

Recommended Posts

[Study memo] Build GeoJSON vector tile server with Fast API
Build a Fast API environment with docker-compose
[Fast API + Firebase] Build an API server for Bearer authentication
Build a speed of light web API server with Falcon
Persist Flask API server with forever