One of the most popular extensions to PostgreSQL is PostGIS, which adds support for storing geospatial geometries, as well as functionality for reasoning about and performing operations on those geometries.
This is a demo showing how to assemble ibis expressions for a PostGIS-enabled database. We will be using a database that has been loaded with an Open Street Map extract for Southern California. This extract can be found here, and loaded into PostGIS using a tool like ogr2ogr.
In [1]:
# Launch the postgis container.
# This may take a bit of time if it needs to download the image.
!docker run -d -p 5432:5432 --name postgis-db -e POSTGRES_PASSWORD=supersecret mdillon/postgis:9.6-alpine
Next, we download our OSM extract (about 400 MB):
In [2]:
!wget https://download.geofabrik.de/north-america/us/california/socal-latest.osm.pbf
Finally, we load it into the database using ogr2ogr (this may take some time):
In [3]:
!ogr2ogr -f PostgreSQL PG:"dbname='postgres' user='postgres' password='supersecret' port=5432 host='localhost'" -lco OVERWRITE=yes --config PG_USE_COPY YES socal-latest.osm.pbf
In [4]:
import os
import geopandas
import ibis
%matplotlib inline
In [5]:
client = ibis.postgres.connect(
url='postgres://postgres:supersecret@localhost:5432/postgres'
)
Let's look at the tables available in the database:
In [6]:
client.list_tables()
Out[6]:
As you can see, this Open Street Map extract stores its data according to the geometry type. Let's grab references to the polygon and line tables:
In [7]:
polygons = client.table('multipolygons')
lines = client.table('lines')
In [8]:
cities = polygons[polygons.admin_level == '8']
cities = cities[
cities.name.name('city_name'),
cities.wkb_geometry.name('city_geometry')
]
We can assemble a specific query for the city of Los Angeles, and execute it to get the geometry of the city. This will be useful later when reasoning about other geospatial relationships in the LA area:
In [9]:
los_angeles = cities[cities.city_name == 'Los Angeles']
la_city = los_angeles.execute()
la_city_geom = la_city.iloc[0].city_geometry
la_city_geom
Out[9]:
Let's also extract freeways from the lines table,
which are indicated by the value 'motorway' in the highway column:
In [10]:
highways = lines[(lines.highway == 'motorway')]
highways = highways[
highways.name.name('highway_name'),
highways.wkb_geometry.name('highway_geometry')
]
Let's test a spatial join by selecting all the highways that intersect the city of Los Angeles, or one if its neighbors.
We begin by assembling an expression for Los Angeles and its neighbors. We consider a city to be a neighbor if it has any point of intersection (by this critereon we also get Los Angeles itself).
We can pass in the city geometry that we selected above when making our query by marking it as a literal value in ibis:
In [11]:
la_neighbors_expr = cities[
cities.city_geometry.intersects(
ibis.literal(la_city_geom, type='multipolygon;4326:geometry')
)
]
In [12]:
la_neighbors = la_neighbors_expr.execute().dropna()
la_neighbors
Out[12]:
Now we join the neighbors expression with the freeways expression, on the condition that the highways intersect any of the city geometries:
In [13]:
la_highways_expr = highways.inner_join(
la_neighbors_expr,
highways.highway_geometry.intersects(la_neighbors_expr.city_geometry)
).materialize()
In [14]:
la_highways = la_highways_expr.execute()
la_highways.plot()
Out[14]:
In [15]:
ocean = geopandas.read_file('https://www.naturalearthdata.com/http//www.naturalearthdata.com/download/10m/physical/ne_10m_ocean.zip')
land = geopandas.read_file('https://www.naturalearthdata.com/http//www.naturalearthdata.com/download/10m/physical/ne_10m_land.zip')
In [16]:
ax = la_neighbors.dropna().plot(figsize=(16,16), cmap='rainbow', alpha=0.9)
ax.set_autoscale_on(False)
ax.set_axis_off()
land.plot(ax=ax, color='tan', alpha=0.4)
ax = ocean.plot(ax=ax, color='navy')
la_highways.plot(ax=ax, color='maroon')
Out[16]: