Introduction to geopandas
and cartopy
#
Basic Setup#
Again, we will be using pandas
and matplotlib
in this tutorial.
import pandas as pd
import matplotlib.pyplot as plt
Note
If you have not yet set up Python on your computer, you can execute this tutorial in your browser via Google Colab. Click on the rocket in the top right corner and launch “Colab”. If that doesn’t work download the .ipynb
file and import it in Google Colab.
Then install the following packages by executing the following command in a Jupyter cell at the top of the notebook.
!pip install pandas geopandas matplotlib cartopy mapclassify
Why do we need something other than pandas
?#
Let’s reload again our example dataset of conventional power plants in Europe as a pd.DataFrame
.
fn = "https://raw.githubusercontent.com/PyPSA/powerplantmatching/master/powerplants.csv"
ppl = pd.read_csv(fn, index_col=0)
This dataset includes coordinates (latitude and longitude), which allows us to plot the location and capacity of all power plants in a scatter plot:
ppl.plot.scatter("lon", "lat", s=ppl.Capacity / 1e3)
<Axes: xlabel='lon', ylabel='lat'>
However, this graphs misses some geographic reference point, we’d normally expect for a map like shorelines, country borders etc.
Geopandas - a Pandas extension for geospatial data#
Geopandas extends pandas
by adding support for geospatial data.
The core data structure in GeoPandas is the geopandas.GeoDataFrame
, a subclass of pandas.DataFrame
, that can store geometry columns and perform spatial operations.
Note
Documentation for this package is available at https://geopandas.org/en/stable/.
Typical geometries are points, lines, and polygons. They come from another library called shapely
, which helps you create, analyze, and manipulate two-dimensional shapes and their properties, such as points, lines, and polygons.
First, we need to import the geopandas
package. The conventional alias is gpd
:
import geopandas as gpd
We can convert the latitude and longitude values given in the dataset to formal geometries (to be exact: shapely.Point
objects but we won’t go into detail regarding this) using the gpd.points_from_xy()
function, and use this to gpd.GeoDataFrame
. We should also specify a so-called coordinate reference system (CRS). The code ‘4326’ means latitude and longitude values.
geometry = gpd.points_from_xy(ppl["lon"], ppl["lat"])
gdf = gpd.GeoDataFrame(ppl, geometry=geometry, crs=4326)
Now, the gdf
looks like this:
gdf.head(3)
Name | Fueltype | Technology | Set | Country | Capacity | Efficiency | DateIn | DateRetrofit | DateOut | lat | lon | Duration | Volume_Mm3 | DamHeight_m | StorageCapacity_MWh | EIC | projectID | geometry | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
id | |||||||||||||||||||
0 | Borssele | Hard Coal | Steam Turbine | PP | Netherlands | 485.0 | NaN | 1973.0 | NaN | 2034.0 | 51.4332 | 3.7160 | NaN | 0.0 | 0.0 | 0.0 | {'49W000000000054X'} | {'BEYONDCOAL': {'BEYOND-NL-2'}, 'ENTSOE': {'49... | POINT (3.71600 51.43320) |
1 | Sainte Croix | Hydro | Reservoir | Store | France | 132.3 | NaN | 1974.0 | NaN | NaN | 43.7375 | 6.1343 | NaN | 0.0 | 0.0 | 0.0 | {'17W100P100P0297X'} | {'ENTSOE': {'17W100P100P0297X'}, 'GEM': {'G601... | POINT (6.13430 43.73750) |
2 | Pied De Borne | Hydro | Reservoir | Store | France | 109.4 | NaN | 1965.0 | NaN | NaN | 44.4788 | 3.9858 | NaN | 0.0 | 0.0 | 0.0 | {'17W100P100P0289W'} | {'ENTSOE': {'17W100P100P0289W'}, 'GEM': {'G601... | POINT (3.98580 44.47880) |
With the additional geometry
columns, it is now even easier to plot the geographic data:
gdf.plot(
column="Fueltype",
markersize=gdf.Capacity / 1e2,
)
<Axes: >
We can also start up an interactive map to explore the geodata in more detail:
gdf.explore(column="Fueltype")
Map Projections with Cartopy#
Cartopy is a Python package designed for geospatial data processing and has exposed an interface to enable easy map creation using matplotlib
.
The Earth is a globe, but we present maps usually on two-dimensional surfaces. Hence, we typically need to project data points onto flat surfaces (e.g. screens, paper). However, we will always loose some information in doing so.
A map projection is:
a systematic transformation of the latitudes and longitudes of locations from the surface of a sphere or an ellipsoid into locations on a plane. Wikipedia: Map projection.
Different projections preserve different metric properties. As a result, converting geodata from one projection to another is a common exercise in geographic data science.
conformal projections preserve angles/directions (e.g. Mercator projection)
equal-area projections preserve area measure (e.g. Mollweide)
equidistant projections preserve distances between points (e.g. Plate carrée)
compromise projections seek to strike a balance between distortions (e.g. Robinson)
If you like the “Orange-as-Earth” analogy for projections, checkout this numberphile video by Hannah Fry.
Note
Documentation for this package is available at https://scitools.org.uk/cartopy/docs/latest/.
First, we need to import the relevant parts of the cartopy
package:
import cartopy
import cartopy.crs as ccrs
Let’s draw a first map with cartopy
outlining the global coastlines in the so-called plate carrée projection (equirectangular projection):
ax = plt.axes(projection=ccrs.PlateCarree())
ax.coastlines()
<cartopy.mpl.feature_artist.FeatureArtist at 0x7f8358c8b990>
A list of the available projections can be found on the Cartopy projection list page.
ax = plt.axes(projection=ccrs.Mollweide())
ax.stock_img()
<matplotlib.image.AxesImage at 0x7f8358cc8a50>
We can combine the functionality of cartopy
with geopandas
plots:
fig = plt.figure(figsize=(7, 7))
ax = plt.axes(projection=ccrs.PlateCarree())
gdf.plot(
ax=ax,
column="Fueltype",
markersize=gdf.Capacity / 1e2,
)
<GeoAxes: >
We can add further geographic features to this map for better orientation.
For instance, we can add the coastlines…
ax.coastlines()
fig
… country borders …
ax.add_feature(cartopy.feature.BORDERS, color="grey", linewidth=0.5)
fig
… colour in the ocean in blue …
ax.add_feature(cartopy.feature.OCEAN, color="azure")
fig
…and color in the land area in yellow …
ax.add_feature(cartopy.feature.LAND, color="cornsilk")
fig
Geopandas will automatically calculate sensible bounds for the plot given the geographic data.
But we can also manually zoom in or out by setting the spatial extent with the .set_extent()
method:
ax.set_extent([5, 16, 47, 55])
fig
Reprojecting a GeoDataFrame
#
In geopandas
, we can use the function .to_crs()
to convert a GeoDataFrame
to a desired coordinate reference system. In this particular case, we use the proj4_init
string of an initialised cartopy
projection to reproject our power plant GeoDataFrame
.
A
proj4_init
string is a text-based representation of a coordinate reference system (CRS) that defines the parameters for transforming geographical coordinates between different spatial reference systems, used by the PROJ library. It will look similar to this: “+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs”.
fig = plt.figure(figsize=(7, 7))
crs = ccrs.AlbersEqualArea()
ax = plt.axes(projection=crs)
gdf.to_crs(crs.proj4_init).plot(
ax=ax,
column="Fueltype",
markersize=gdf.Capacity / 1e2,
)
ax.coastlines()
<cartopy.mpl.feature_artist.FeatureArtist at 0x7f8355b93110>
Reading and Writing Files with geopandas
#
In the following example, we’ll load a dataset containing the NUTS regions:
Nomenclature of Territorial Units for Statistics or NUTS (French: Nomenclature des unités territoriales statistiques) is a geocode standard for referencing the subdivisions of countries for statistical purposes.
Our ultimate goal for this part of the tutorial is to map the power plant capacities to the NUTS-1 region they belong to.
Common filetypes for vector-based geospatial datasets are GeoPackage (.gpkg
), GeoJSON (.geojson
), File Geodatabase (.gdb
), or Shapefiles (.shp
).
In geopandas
we can use the gpd.read_file()
function to read such files. So let’s start:
url = "https://tubcloud.tu-berlin.de/s/RHZJrN8Dnfn26nr/download/NUTS_RG_10M_2021_4326.geojson"
nuts = gpd.read_file(url)
nuts.head(3)
id | NUTS_ID | LEVL_CODE | CNTR_CODE | NAME_LATN | NUTS_NAME | MOUNT_TYPE | URBN_TYPE | COAST_TYPE | FID | geometry | |
---|---|---|---|---|---|---|---|---|---|---|---|
0 | BG423 | BG423 | 3 | BG | Pazardzhik | Пазарджик | 3.0 | 2 | 3 | BG423 | POLYGON ((24.42101 42.55306, 24.41032 42.46950... |
1 | BG424 | BG424 | 3 | BG | Smolyan | Смолян | 3.0 | 3 | 3 | BG424 | POLYGON ((25.07422 41.79348, 25.05851 41.75177... |
2 | BG425 | BG425 | 3 | BG | Kardzhali | Кърджали | 3.0 | 3 | 3 | BG425 | POLYGON ((25.94863 41.32034, 25.90644 41.30757... |
It is good practice to set an index. You can use .set_index()
for that:
nuts = nuts.set_index("id")
We can also check out the geometries in the dataset with .geometry
:
nuts.geometry
id
BG423 POLYGON ((24.42101 42.55306, 24.41032 42.46950...
BG424 POLYGON ((25.07422 41.79348, 25.05851 41.75177...
BG425 POLYGON ((25.94863 41.32034, 25.90644 41.30757...
CH011 MULTIPOLYGON (((6.86623 46.90929, 6.89621 46.9...
CH012 POLYGON ((8.47767 46.52760, 8.39953 46.48872, ...
...
LV POLYGON ((27.35158 57.51824, 27.54521 57.53444...
ME POLYGON ((20.06394 43.00682, 20.32958 42.91149...
MK POLYGON ((22.36021 42.31116, 22.51041 42.15516...
SK0 POLYGON ((19.88393 49.20418, 19.96275 49.23031...
IT MULTIPOLYGON (((12.47792 46.67984, 12.69064 46...
Name: geometry, Length: 2010, dtype: geometry
With .crs
we can check in which coordinate reference system the data is given:
nuts.crs
<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World.
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich
nuts.total_bounds
array([-63.08825, -21.38917, 55.83616, 80.76427])
Let’s filter by NUTS-1 level…
nuts1 = nuts.query("LEVL_CODE == 1")
… and explore what kind of geometries we have in the dataset …
nuts1.explore()
To write a GeoDataFrame back to file use GeoDataFrame.to_file()
. The file format is inferred from the file ending.
nuts1.to_file("tmp.geojson")
Calculating the areas and buffers#
The first thing we need to do to calculate area or buffers is to reproject the GeoDataFrame
to an equal-area projection (here: EPSG:3035 which is valid only within Europe; global alternative is the Mollweide projection EPSG:54009):
nuts1 = nuts1.to_crs(3035)
The area can be accessed via .area
and is given in m² (after projection). Let’s convert to km²:
area = nuts1.area / 1e6
area
id
AT1 23545.286205
AT2 25894.953057
EL4 17388.679384
EE0 45315.713593
EL3 3799.676547
...
PL7 29846.398582
PL8 63217.536546
PL9 35563.812826
RO2 72545.290405
SK0 49008.115415
Length: 125, dtype: float64
nuts1.explore(column=area, vmax=1e5)
We can also build a buffer of 1km around each geometry using .buffer()
:
nuts1.buffer(1000).explore()
Joining spatial datasets#
Multiple GeoDataFrames
can be combined via spatial joins.
Observations from two datasets are combined with the .sjoin()
function based on their spatial relationship to one another (e.g. whether they are intersecting or overlapping). You can read more about the specific options here.
To perform a spatial join, you need to provide the following information:
The two GeoDataFrames you want to join (e.g.,
left_df
andright_df
).The type of spatial relationship to test (e.g.,
intersects
,contains
,within
). This is specified using theop
parameter.The type of join to perform (e.g.,
inner
,left
,right
). This is specified using thehow
parameter.
The .sjoin()
function will then iterate through the geometries in both GeoDataFrames, evaluate the specified spatial relationship, and join the matching records.
For example, if the op
parameter is set to ‘intersects’, the function will check if the geometries of each record in left_df
intersect with the geometries of any records in right_df
. If a match is found, the attributes of the corresponding records from both GeoDataFrames will be combined into a new record in the output GeoDataFrame.
By performing spatial joins, you can efficiently combine and analyze geospatial data based on their spatial relationships, without the need for explicit coordinate-based calculations.
Let’s first reproject the gdf
object to the same CRS as nuts1
:
gdf = gdf.to_crs(3035)
Then, let’s have a look at both datasets at once. We want to find out which points (representing power plants) lie within which shape (representing NUTS regions).
fig = plt.figure(figsize=(7, 7))
ax = plt.axes(projection=ccrs.epsg(3035))
nuts1.plot(ax=ax, edgecolor="black", facecolor="lightgrey")
gdf.to_crs(3035).plot(
ax=ax, column="Fueltype", markersize=gdf.Capacity / 20, legend=True
)
ax.set_extent([5, 19, 47, 55])
We can now apply the .sjoin
function to look for which power plants lie within which NUTS1 region. By default, sjoin
looks for intersections and keeps the geometries of the left GeoDataFrame
.
joined = gdf.sjoin(nuts1)
If we look at this new GeoDataFrame
, we now have additional columns from the NUTS1 data:
joined.head(3)
Name | Fueltype | Technology | Set | Country | Capacity | Efficiency | DateIn | DateRetrofit | DateOut | ... | index_right | NUTS_ID | LEVL_CODE | CNTR_CODE | NAME_LATN | NUTS_NAME | MOUNT_TYPE | URBN_TYPE | COAST_TYPE | FID | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
id | |||||||||||||||||||||
1 | Sainte Croix | Hydro | Reservoir | Store | France | 132.3 | NaN | 1974.0 | NaN | NaN | ... | FRL | FRL | 1 | FR | Provence-Alpes-Côte d’Azur | Provence-Alpes-Côte d’Azur | 0.0 | 0 | 0 | FRL |
2 | Pied De Borne | Hydro | Reservoir | Store | France | 109.4 | NaN | 1965.0 | NaN | NaN | ... | FRJ | FRJ | 1 | FR | Occitanie | Occitanie | 0.0 | 0 | 0 | FRJ |
3 | Pouget | Hydro | Reservoir | Store | France | 446.9 | NaN | 1951.0 | NaN | NaN | ... | FRJ | FRJ | 1 | FR | Occitanie | Occitanie | 0.0 | 0 | 0 | FRJ |
3 rows × 29 columns
We can now use these new columns to group the capacities (and convert to a suitable unit):
cap = joined.groupby("NUTS_ID").Capacity.sum() / 1000 # GW
Let’s quickly check if all NUTS1 regions have power plants:
nuts1.index.difference(cap.index)
Index(['CY0', 'ES7', 'FI2', 'FRY', 'IS0', 'LI0', 'MK0', 'MT0', 'PT2', 'TR1',
'TR2', 'TR3', 'TR4', 'TR5', 'TR6', 'TR7', 'TR8', 'TR9', 'TRA', 'TRB',
'TRC'],
dtype='object')
This is not the case. Then it is good practice to reindex the series to include all NUTS1 regions, even if this leads to some NaN values.
cap = cap.reindex(nuts1.index)
cap
id
AT1 4.463100
AT2 4.896500
EL4 0.427967
EE0 2.155000
EL3 1.343142
...
PL7 7.077524
PL8 1.805551
PL9 6.033288
RO2 3.514566
SK0 5.713420
Name: Capacity, Length: 125, dtype: float64
Finally, we can plot the total generation capacity per NUTS1 region on a map.
nuts1.plot(figsize=(7, 7), column=cap, legend=True)
<Axes: >
This concludes the geopandas
and cartopy
tutorial.
Exercises#
Task 1: Recreate the figure above (i.e. generation capacity per NUTS1 region)
using 3 different cartopy projections of your choice.
such that the capacities are normalised by the area of the NUTS1 region (unit: MW/km²).
such that it only shows the hard coal power plant capacities.
such that it only shows NUTS1 shapes for Germany (and Denmark).
but with NUTS2 regions instead of NUTS1 regions.
(advanced) such that it shows the capacity-weighted average age of the power plants instead of the capacities. Where are the oldest generators? Can you think of a reason why?
Which German NUTS2 region has the highest conventional generation capacity? Can you think of a reason why?
Three projections of choice:
Show code cell content
for crs in [ccrs.EqualEarth(), ccrs.AlbersEqualArea(), ccrs.Orthographic()]:
fig = plt.figure()
ax = plt.axes(projection=crs)
nuts1.to_crs(crs.proj4_init).plot(ax=ax, column=cap, legend=True)
Normalised by area in MW/km²:
Show code cell content
nuts1.plot(figsize=(7, 7), column=cap / (nuts1.area / 1e9), legend=True)
<Axes: >
Only hard coal generators:
Show code cell content
hard_coal_cap = (
joined.query("Fueltype == 'Hard Coal'")
.groupby("NUTS_ID")
.Capacity.sum()
.reindex(nuts1.index)
.div(1e3)
)
nuts1.plot(
figsize=(7, 7),
column=hard_coal_cap,
legend=True,
missing_kwds=dict(color="lightgrey"),
)
plt.ylim(1e6, 6e6)
plt.xlim(2e6, 6e6);
Only Germany:
Show code cell content
subregion = nuts1.query("CNTR_CODE == 'DE'")
subregion.plot(column=cap.reindex(subregion.index), legend=True)
<Axes: >
Only Germany and Denmark:
Show code cell content
countries = ["DE", "DK"]
subregion = nuts1.query("CNTR_CODE in @countries") # alternative A
subregion = nuts1.loc[nuts1.CNTR_CODE.isin(countries)] # alternative B
subregion.plot(column=cap.reindex(subregion.index), legend=True)
<Axes: >
In NUTS2 rather than NUTS1:
Show code cell content
nuts2 = nuts.query("LEVL_CODE == 2").to_crs(3035)
joined2 = gdf.sjoin(nuts2)
cap2 = joined2.groupby("NUTS_ID").Capacity.sum().reindex(nuts2.index).div(1000) # GW
nuts2.plot(figsize=(7, 7), column=cap2, legend=True)
<Axes: >
NUTS2 region with highest capacity:
Show code cell content
cap2.sort_values(ascending=False).head(1) # alternative A
id
FRK2 22.271118
Name: Capacity, dtype: float64
Show code cell content
cap2.idxmax() # alternative B
'FRK2'
Capacity-weighted average age:
Show code cell content
def calculate_weighted_age(x):
"""
Formula:
\bar{x} = \frac{ \sum\limits_{i=1}^n w_i x_i}{\sum\limits_{i=1}^n w_i}
"""
numerator = (x.Capacity * x.DateIn).sum()
# only build sum of capacities where a date is given
denominator = x.Capacity.where(~x.DateIn.isna()).sum()
return numerator / denominator
age = joined.groupby("NUTS_ID").apply(calculate_weighted_age).reindex(nuts1.index)
nuts1.plot(figsize=(7, 7), column=age, legend=True)
<Axes: >
Task 2: Load the following .gpkg
file containing the Danish Natura2000 natural protection areas as GeoDataFrame
:
https://tubcloud.tu-berlin.de/s/mEpdmgBtmMbyjAr/download/Natura2000_end2021-DK.gpkg
How many individual natural protection areas are there?
What is the name of the largest natural protection area in Denmark?
What is the total protection area in square kilometers (make sure you don’t double-count overlapping areas by running the
.dissolve()
function).(advanced) Use set operations with the
.overlay()
function and the NUTS1 regionsGeoDataFrame
to identify what share of Danish natural protection areas is at sea.
Show code cell content
url = "https://tubcloud.tu-berlin.de/s/mEpdmgBtmMbyjAr/download/Natura2000_end2021-DK.gpkg"
Show code cell content
natura = gpd.read_file(url)
Show code cell content
natura.plot(alpha=0.5)
<Axes: >
Show code cell content
natura.crs
<Projected CRS: EPSG:3035>
Name: ETRS89-extended / LAEA Europe
Axis Info [cartesian]:
- Y[north]: Northing (metre)
- X[east]: Easting (metre)
Area of Use:
- name: Europe - European Union (EU) countries and candidates. Europe - onshore and offshore: Albania; Andorra; Austria; Belgium; Bosnia and Herzegovina; Bulgaria; Croatia; Cyprus; Czechia; Denmark; Estonia; Faroe Islands; Finland; France; Germany; Gibraltar; Greece; Hungary; Iceland; Ireland; Italy; Kosovo; Latvia; Liechtenstein; Lithuania; Luxembourg; Malta; Monaco; Montenegro; Netherlands; North Macedonia; Norway including Svalbard and Jan Mayen; Poland; Portugal including Madeira and Azores; Romania; San Marino; Serbia; Slovakia; Slovenia; Spain including Canary Islands; Sweden; Switzerland; Türkiye (Turkey); United Kingdom (UK) including Channel Islands and Isle of Man; Vatican City State.
- bounds: (-35.58, 24.6, 44.83, 84.73)
Coordinate Operation:
- name: Europe Equal Area 2001
- method: Lambert Azimuthal Equal Area
Datum: European Terrestrial Reference System 1989 ensemble
- Ellipsoid: GRS 1980
- Prime Meridian: Greenwich
How many:
Show code cell content
len(natura)
350
Name of largest protection area:
Show code cell content
i = natura.area.idxmax()
i
309
Show code cell content
natura.loc[i]
SITECODE DK00FX112
SITENAME Skagens Gren og Skagerak
MS DK
SITETYPE B
INSPIRE_ID Dk.nst.ps.SAC1
geometry POLYGON ((4323183.452199999 3832152.5764000034...
Name: 309, dtype: object
Total area in km²:
Show code cell content
total = natura.dissolve()
Show code cell content
total.plot(alpha=0.5)
<Axes: >
Show code cell content
total.area / 1e6 # sqkm
0 22646.375701
dtype: float64
Share of natural protection area offshore:
Show code cell content
onshore_natura = gpd.overlay(nuts1.dissolve(), total)
Show code cell content
onshore_natura.explore()
Show code cell content
1 - onshore_natura.area[0] / total.area[0]
0.8133701227617232