GeoPandas og PostGIS

---

Referanser som gir bakgunn for å forstå eksemplene nedenfor:

Utvalgte referanser

• geopandas.GeoDataFrame.to_postgis
• geopandas.read_postgis
• pandas.read_sql
• geopandas.GeoSeries.contains
• PostGIS: ST_Contains
• How do I select a subset of a DataFrame
• How do I filter specific rows from a DataFrame?
• How do I select specific columns from a DataFrame?

GeoPandas generelle referanser

• GeoSeries
• GeoDataFrame
• Input/Output
• Tools - inkludert Spatial Join
• API Reference

PostGIS generelle referanser

• 3.3. Creating spatial databases
• Chapter 4. PostGIS Usage
• 4.2. Spatial Queries
• Chapter 5. PostGIS Reference
• PostGIS hovedside
• Workshop: Introduction to PostGIS

NTNU 05.03.2021 - Sverre Stikbakke

---

GeoDataFrame: teiger

import pandas as pd
import geopandas as gpd
import fiona as fi
from shapely.geometry import Polygon, Point

# Python dictionary (en innebygd datastruktur i Python)

data = {'smaragd': {
    'localid': 446124344,
    'objtype': 'Teig',
    'SHAPE_Length': 306.207095792018,
    'SHAPE_Area': 5747.21023925478,
    'geometry': Polygon([
        [591518.2, 6740464.0499],
        [591442.74, 6740401.1999],
        [591399.69, 6740452.8699],
        [591460.55, 6740518.9299],
        [591464.35, 6740525.3999],
        [591471.27, 6740517.5099],
        [591516.35, 6740466.1499],
        [591518.2, 6740464.0499],
        ]),
    }, 'beryll': {
    'localid': 49174530,
    'objtype': 'Teig',
    'SHAPE_Length': 306.207095792018,
    'SHAPE_Area': 5747.21023925478,
    'geometry': Polygon([
        [591535.31, 6740542.8499],
        [591546.79, 6740525.9499],
        [591565.3179, 6740498.2083],
        [591527.4193, 6740467.2872],
        [591516.35, 6740466.1499],
        [591471.27, 6740517.5099],
        [591464.35, 6740525.3999],
        [591494.5, 6740576.6499],
        [591514.4, 6740554.4999],
        [591525.4, 6740550.7899],
        [591535.31, 6740542.8499],
        ]),
    }}

teiger = gpd.GeoDataFrame.from_dict(data, orient='index')
teiger.set_crs('EPSG:25832')

	localid	objtype	SHAPE_Length	SHAPE_Area	geometry
smaragd	446124344	Teig	306.207096	5747.210239	POLYGON ((591518.200 6740464.050, 591442.740 6...
beryll	49174530	Teig	306.207096	5747.210239	POLYGON ((591535.310 6740542.850, 591546.790 6...

teiger.plot(figsize=(8,8))

<AxesSubplot:>

GeoDataFrame: bygninger

# skjuler warnings i cellen nedenfor
import warnings
warnings.filterwarnings('ignore')

# fila som importeres her har Geopackage-format. Den er eksportert fra QGIS.
bygninger = gpd.read_file('bygninger_ved_ntnu.gpkg')

bygninger.plot(figsize=(8,8))

<AxesSubplot:>

Postgis

Postgis-operasjoner med bruk av sqlalcemy, geoalchemy2 og psycopg2. Se Installasjon og oppstart av Geopandas i Windows 10 for info om installasjon av nødvendige Python-pakker.

Brukernavn og passord til databasen, samt vertsnavn (host), portnummer og databasenavn legges i en egen fil, mypostgisdb.py, slik:

postgis_connection_string = "postgresql://brukernavn:passord@host:5432/databasenavn"

from sqlalchemy import create_engine
from mypostgisdb import postgis_connection_string

db = create_engine(postgis_connection_string)

# Disse linjene laster opp teiger og bygninger (GeoDataFrames) til databasen og oppretter tabeller med angitte navn.
# teiger.to_postgis('teiger', db) 
# bygninger.to_postgis('bygninger', db)

SQL-spørringer som ikke returnerer geometri: bruk pd.read_sql()

Enkle SQL-spørringer mot databasen gjøres med Pandas-funksjonen read_sql(). Spørringene vil returnere en Pandas DataFrame.

sql = 'SELECT * from teiger'
pandas_df_teiger = pd.read_sql(sql, db)
pandas_df_teiger

	localid	objtype	SHAPE_Length	SHAPE_Area	geometry
0	446124344	Teig	306.207096	5747.210239	01030000000100000008000000666666663C0D2241C58F...
1	49174530	Teig	306.207096	5747.210239	0103000000010000000B000000EC51B89E5E0D2241F8C2...

# Geometri-kolonnen ble med her, men datatypen for enkelt-objektene er _str_. Den kan ikke brukes videre til geometri-operasjoner.
type(pandas_df_teiger.loc[0,'geometry'])

str

Vi kan bruke PostGIS-funksjoner i databasen med Pandas read_sql() dersom de ikke returnerer geometri-data:

sql = 'SELECT localid, ST_Area(geometry) FROM teiger'
pandas_df_area = pd.read_sql(sql, db)
pandas_df_area

	localid	st_area
0	446124344	7166.961750
1	49174530	5747.210239

SQL-spørringer som returnerer geometri: bruk gpd.read_postgis()

GeoPandas-funksjonen read_postgis() krever at SQL-spørringen skal returnere en geometri-kolonne.
Det vil si at SELECT-delen av SQL-spørringen må inkludere en geometri-kolonne, enten eksplisitt eller implisitt som i SELECT *.
Merk også at kolonnenavn som inneholder store bokstaver må settes i anførselstegn.

Dette fungerer:

sql = 'SELECT * FROM teiger'

Dette fungerer ikke - gir feilmelding:

sql = 'SELECT localid, "SHAPE_Length" FROM teiger'

Dette fungerer:

sql = 'SELECT localid, "SHAPE_Length", geometry FROM teiger'

# Disse variablene brukes i spørringer med gpd.read_postgis()
geom_col = 'geometry'
crs = 'EPSG:25832'

sql = 'SELECT * FROM teiger'
geopandas_df_teiger = gpd.read_postgis(sql, db, geom_col, crs)
geopandas_df_teiger

	localid	objtype	SHAPE_Length	SHAPE_Area	geometry
0	446124344	Teig	306.207096	5747.210239	POLYGON ((591518.200 6740464.050, 591442.740 6...
1	49174530	Teig	306.207096	5747.210239	POLYGON ((591535.310 6740542.850, 591546.790 6...

# Geometri-objektene har nå datatypen _shapely.geometry.polygon.Polygon_
type(geopandas_df_teiger.loc[0,'geometry'])

shapely.geometry.polygon.Polygon

sql = 'SELECT localid, "SHAPE_Length", geometry FROM teiger'
geopandas_df_teiger_length = gpd.read_postgis(sql, db, geom_col, crs)
geopandas_df_teiger_length

	localid	SHAPE_Length	geometry
0	446124344	306.207096	POLYGON ((591518.200 6740464.050, 591442.740 6...
1	49174530	306.207096	POLYGON ((591535.310 6740542.850, 591546.790 6...

Sammenligning: Egenhofer-topologi-spørring med PostGIS-funksjoner vs. GeoPandas-funksjoner

Med spørringen vil vi finne de bygningene som ligger innenfor teig-polygonene for Beryll- og Smaragd-byggene.

• PostGIS-varianten forutsetter at dataene er lastet opp til tabeller i en PostGIS-database.
• PostGIS-spørringen returnerer en komplett GeoDataFrame med de kolonnene som er spesifiert i spørringen.
• GeoPandas-varianten returnerer en Pandas Series med True- eller False-verdier. Denne legges direkte inn i en eksisterende GeoDataFrame for bygninger.
• Generell Pandas-funksjonalitet brukes for å hente ut radene med True-verdier, og kolonnene objectid og geometry.

Se Egenhofer-topologi for nærmere beskrivelse av GeoPandas-varianten.

sql = 'SELECT objectid, bygninger.geometry FROM bygninger, teiger WHERE ST_CONTAINS(teiger.geometry, bygninger.geometry)'
utvalg_bygninger_postgis = gpd.read_postgis(sql, db, geom_col, crs)
utvalg_bygninger_postgis

	objectid	geometry
0	3140	MULTIPOLYGON Z (((591464.400 6740433.940 184.9...
1	29039	MULTIPOLYGON Z (((591457.090 6740476.920 191.2...
2	9906	MULTIPOLYGON Z (((591482.390 6740523.390 198.8...
3	34964	MULTIPOLYGON Z (((591516.280 6740536.880 190.1...

utvalg_bygninger_postgis.plot()

<AxesSubplot:>

# oppdaterer Geopandas GeoDataFrame bygninger med en ny kolonne, Beryll_Smaragd, med True/False, som angir om bygningene ligger innenfor noen av teig-polygonene 
bygninger['Beryll_Smaragd'] = bygninger.geometry.apply(lambda bygning: teiger.geometry.contains(bygning).any())
utvalg_bygninger_geopandas = bygninger[bygninger.Beryll_Smaragd == True][["objectid", "geometry"]]
utvalg_bygninger_geopandas

	objectid	geometry
91	3140	MULTIPOLYGON Z (((591464.400 6740433.940 184.9...
283	9906	MULTIPOLYGON Z (((591482.390 6740523.390 198.8...
834	29039	MULTIPOLYGON Z (((591457.090 6740476.920 191.2...
1048	34964	MULTIPOLYGON Z (((591516.280 6740536.880 190.1...

utvalg_bygninger_geopandas.plot()

<AxesSubplot:>

utvalg_bygninger_geopandas.loc[91][['geometry']]

geometry    (POLYGON Z ((591464.3997999998 6740433.9399 18...
Name: 91, dtype: object