Рисование геометрий с помощью GeoPandas

28 Apr 2020, 12:02

Ведение

На днях у меня возникла потребность отобразить несколько геометрий на карте, для того чтобы посмотреть на сколько они перекрывают друг друга.

Варианты решения

Первое что мне пришло в голову - это скачать QGIS, это геоинформационная система с открытым исходным кодом, которая бы выполнила мою задачу, но она вести около 1 Gb места, и мне показалось что для такой простой задачи это большой оверхэд.

Следующее решение было - сделать html файлик и использовать javascript библиотеку openlayers для отображения моих геометрий. Итого у меня получилось следующее:

<!doctype html>
<html lang="en">
  <head>
    <link rel="stylesheet" href="https://cdn.jsdelivr.net/gh/openlayers/openlayers.github.io@master/en/v6.3.1/css/ol.css" type="text/css">
    <style>
      .map {
        height: 400px;
        width: 100%;
      }
    </style>
    <script src="https://cdn.jsdelivr.net/gh/openlayers/openlayers.github.io@master/en/v6.3.1/build/ol.js"></script>
    <title>OpenLayers example</title>
  </head>
  <body>
    <h2>My Map</h2>
    <div id="map" class="map"></div>
    <script type="text/javascript">

    var geojsonObject = {"type": "FeatureCollection",
      "features": [
        {"type":"Polygon","coordinates":...},
        {"type":"LineString","coordinates":...}
      ]
    }

      var map = new ol.Map({
        target: 'map',
        layers: [
          new ol.layer.Tile({
            source: new ol.source.OSM()
          }),
          new ol.layer.Vector({
            source: new ol.source.Vector({
                features: (new ol.format.GeoJSON()).readFeatures(geojsonObject)
            })
          })
        ],
        view: new ol.View({
          center: [2186457.60376491,7493014.69347035],
          zoom: 17
        })
      });
    </script>
  </body>
</html>

Но тут меня смущало, что геометрии находятся в БД и для того чтобы их отобразить, надо ручками перенести их в этот файл.

Затем я посоветовался со своим другом и он подсказал мне библиотеку про которую я совсем забыл - geopandas. Эта библиотека дает возможность работать расширить стандартный DataFrame из pandas‘a, функциями для работы с геоданными.

Работа с Geopandas

Первое что нужно было сделать это получить нужные геометрии из БД и загрузить их в DataFrame:

import pandas as pd
import contextily as ctx
from shapely import wkt
import geopandas
import psycopg2

db = psycopg2.connect("dbname=vts user=postgres password=Xae4aep9 host=10.127.32.41 port=5432")
cur = db.cursor()

cur.execute("""select id, st_astext(ST_Transform(geometry, 3857))
             from test
             where id in (4050079, 130558)""")

df = pd.DataFrame(cur.fetchall(), columns =['id', 'geom'])

Чтобы для простого запроса не тащить всю SQL Alchemy, через которую pandas работает с БД, я создал DataFrame на основе списка кортежей (результат sql запроса).

Также нужно отметить что геометрию я сначала перевел в проекцию EPSG:3785, а потом преобразовал в формат WKT, для дальнейшей загрузки в GeoPandas.

После этого создает GeoDataFrame:

df['geom'] = df['geom'].apply(wkt.loads)
gdf = geopandas.GeoDataFrame(df, geometry='geom', crs=3857)

Теперь можно отобразить геометрию и посмотреть что получилось:

gdf.plot(figsize=(10, 10), alpha=0.5, edgecolor='k')

В итоге получим что-то похожее на:

Теперь нужно подложить подложку карты (в моем случае это OpenStreetMap) чтобы можно было понять что за объекты:

ax = gdf.plot(figsize=(10, 10), alpha=0.5, edgecolor='k')
ctx.add_basemap(ax, zoom=17, source=ctx.providers.OpenStreetMap.Mapnik)

Итоговая картинка будет выглядеть так:

Для того, чтобы отобразить несколько слоев на карте, нужно использовать plt.subplots примерно так:

f, ax = plt.subplots(figsize=(20, 20))
gdf1.plot(ax=ax)
gdf2.plot(ax=ax)
ctx.add_basemap(ax)

Заключение

В итоге я написал утилиту которая по заданному id объекта выгружает такую картинку, и все вседено в одно простое действие.