%matplotlib inline
import pandas as pd
import geopandas
pd.options.display.max_rows = 10
countries = geopandas.read_file("zip://./data/ne_110m_admin_0_countries.zip")
cities = geopandas.read_file("zip://./data/ne_110m_populated_places.zip")
rivers = geopandas.read_file("zip://./data/ne_50m_rivers_lake_centerlines.zip")
An important aspect of geospatial data is that we can look at spatial relationships: how two spatial objects relate to each other (whether they overlap, intersect, contain, .. one another).
The topological, set-theoretic relationships in GIS are typically based on the DE-9IM model. See https://en.wikipedia.org/wiki/Spatial_relation for more information.
(Image by Krauss, CC BY-SA 3.0)
Let's first create some small toy spatial objects:
A polygon (note: we use .squeeze()
here to to extract the scalar geometry object from the GeoSeries of length 1):
belgium = countries.loc[countries['name'] == 'Belgium', 'geometry'].squeeze()
Two points:
paris = cities.loc[cities['name'] == 'Paris', 'geometry'].squeeze()
brussels = cities.loc[cities['name'] == 'Brussels', 'geometry'].squeeze()
And a linestring:
from shapely.geometry import LineString
line = LineString([paris, brussels])
Let's visualize those 4 geometry objects together (I only put them in a GeoSeries to easily display them together with the geopandas .plot()
method):
geopandas.GeoSeries([belgium, paris, brussels, line]).plot(cmap='tab10')
You can recognize the abstract shape of Belgium.
Brussels, the capital of Belgium, is thus located within Belgium. This is a spatial relationship, and we can test this using the individual shapely geometry objects as follow:
brussels.within(belgium)
And using the reverse, Belgium contains Brussels:
belgium.contains(brussels)
On the other hand, Paris is not located in Belgium:
belgium.contains(paris)
paris.within(belgium)
The straight line we draw from Paris to Brussels is not fully located within Belgium, but it does intersect with it:
belgium.contains(line)
line.intersects(belgium)
The same methods that are available on individual shapely
geometries as we have seen above, are also available as methods on GeoSeries
/ GeoDataFrame
objects.
For example, if we call the contains
method on the world dataset with the paris
point, it will do this spatial check for each country in the world
dataframe:
countries.contains(paris)
Because the above gives us a boolean result, we can use that to filter the dataframe:
countries[countries.contains(paris)]
And indeed, France is the only country in the world in which Paris is located.
Another example, extracting the linestring of the Amazon river in South America, we can query through which countries the river flows:
amazon = rivers[rivers['name'] == 'Amazonas'].geometry.squeeze()
countries[countries.crosses(amazon)] # or .intersects
Overview of the different functions to check spatial relationships (spatial predicate functions):
See https://shapely.readthedocs.io/en/stable/manual.html#predicates-and-relationships for an overview of those methods.
See https://en.wikipedia.org/wiki/DE-9IM for all details on the semantics of those operations.Next to the spatial predicates that return boolean values, Shapely and GeoPandas aslo provide analysis methods that return new geometric objects.
See https://shapely.readthedocs.io/en/stable/manual.html#spatial-analysis-methods for more details.
For example, using the toy data from above, let's construct a buffer around Brussels (which returns a Polygon):
geopandas.GeoSeries([belgium, brussels.buffer(1)]).plot(alpha=0.5, cmap='tab10')
and now take the intersection, union or difference of those two polygons:
brussels.buffer(1).intersection(belgium)
brussels.buffer(1).union(belgium)
brussels.buffer(1).difference(belgium)
Another useful method is the unary_union
attribute, which converts the set of geometry objects in a GeoDataFrame into a single geometry object by taking the union of all those geometries.
For example, we can construct a single object for the Africa continent:
africa_countries = countries[countries['continent'] == 'Africa']
africa = africa_countries.unary_union
africa
print(str(africa)[:1000])
GeoPandas (and Shapely for the individual objects) provides a whole lot of basic methods to analyse the geospatial data (distance, length, centroid, boundary, convex_hull, simplify, transform, ....), much more than the few that we can touch in this tutorial.