Geopandas and Cartopy

Why not pandas?¶

Let’s reload again our example dataset of conventional power plants in Europe as a pd.DataFrame.

This dataset includes coordinates (latitude and longitude), which allows us to plot the location and capacity of all power plants in a scatter plot:

However, this graphs misses some geographic reference point, we’d normally expect for a map like shorelines, country borders and so on.

Why geopandas?¶

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.

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:

We can convert the latitude and longitude values given in the dataset to formal geometries (to be exact, a shapely.Point object) 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, the so-called WGS84 system.

Now, the resulting gdf looks like this:

With the additional geometry columns, it is now even easier to plot the geographic data:

We can also start up an interactive map to explore the geodata in more detail:

Projections with cartopy¶

Cartopy is a Python package designed for geospatial data processing and has an interface for handling projections in matplotlib maps.

Why is it needed?

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.

First, we need to import the relevant parts of the cartopy package:

Let’s draw a first map with cartopy outlining the global coastlines in the so-called plate carrée projection (equirectangular projection):

We can combine the functionality of cartopy with geopandas plots and even add geographical features:

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:

Reprojecting a GeoDataFrame¶

In geopandas, we can use the function .to_crs() to convert a GeoDataFrame with a given projection to another desired target coordinate reference system, using a cartopy CRS object:

Reading and Writing Files¶

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:

It is good practice to set an index. You can use .set_index() for that:

We can also check out the geometries in the dataset with .geometry:

With .crs we can check in which coordinate reference system the data is given:

With .total_bounds we can check the bounding box of the data in the original coordinate reference system:

Let’s now filter the GeoDataFrame by NUTS-1 level...

... and explore what kind of geometries we have in the dataset ...

To write a GeoDataFrame back to file use GeoDataFrame.to_file(). The file format is inferred from the file ending.

Buffers and Areas¶

Many geospatial analyses require the calculation of areas or buffers around geometries. However, these calculations depend on the coordinate reference system used.

The first thing we need to do in order to calculate area or buffers is to reproject the GeoDataFrame to an equal-area projection (here: EPSG:3035 which is valid only within Europe; a globally valid alternative is the Mollweide projection EPSG:54009)

The area can be accessed via .area and is given in m² (after projection). Let’s convert to km²:

We can also build a buffer of 1km around each geometry using .buffer():

It is also possible to calculate the centroid of each geometry, which is the geometric center of a shape. It can be useful for labelling or constructing a point-based network representation of polygonal data:

To calculate the distance between two points (or geometries more generally), we can use the .distance() method. For instance, using the nuts1.representative_point() method, we can find a central point that represents each geometry but other than the centroid is guaranteed to be within the geometry. However, note that the distance is only meaningful if the geometries are in an appropriate projected coordinate reference system (e.g. an equidistant projection):

The .distance() function uses the haversine formula to calculate the great-circle distance between two points on the Earth’s surface, given their latitude and longitude coordinates.

Dissolve¶

The .dissolve() function is used to aggregate geometries within a GeoDataFrame based on a specified attribute or set of attributes. Practically, it is very similar to the pandas.DataFrame.groupby() function, but instead of aggregating numerical values, it merges geometries.

For instance, we can dissolve the NUTS regions by country code to get country-level geometries:

Spatial Joins¶

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, overlapping or contained within one another). You can read more about the specific options here.

Source: https://pythongis.org

To perform a spatial join, you need to provide the following information:

• The two geopandas.GeoDataFrame objects you want to join (e.g., left_df and right_df).

• The type of spatial relationship to test (e.g., intersects, contains, within). This is specified using the predicate parameter.

• The type of join to perform (e.g., inner, left, right). This is specified using the how parameter. When using left, the resulting GeoDataFrame will retain all records from the left_df, while matching records from the right_df will be added where the spatial relationship holds true. With inner, only records with matches in both GeoDataFrames will be retained. With right, all records from the right_df will be retained, with matching records from the left_df added where applicable.

The .sjoin() function will then go through the geometries in both GeoDataFrames, evaluate the specified spatial relationship, and join the matching records (behind the scenes, it uses spatial indexing to speed up the process).

For example, if the predicate 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.

Source: https://pythongis.org

The result of the spatial join is a new GeoDataFrame that retains the geometries from the left_df and includes attributes from both GeoDataFrames for the matching records.

When no match is found, the resulting attributes from the right_df will be set to NaN (e.g. when some power plants are located outside of NUTS-1 regions).

Let’s do an example, starting with first reprojecting the gdf object of power plant to the same CRS as nuts1 polygons and reducing it to coal plants:

Then, let’s have a look at both datasets at once. We want to find out which points (i.e. the power plants) lie within which shape (i.e. the NUTS1 regions).

We can now apply the .sjoin() function. By default, .sjoin() looks for intersections and keeps the geometries of the left GeoDataFrame.

If we look at this new GeoDataFrame, we now have additional columns from the NUTS1 data:

We can now use these new columns for .groupby() operations on the capacities (as well as sorting, converting values to a suitable unit and renaming the index from NUTS identifiers to NUTS region names):

With our intersection and left join, some NUTS1 regions disappear in cap because they do not have any power plants. We can fix this by reindexing cap to the full set of NUTS1 regions:

Additionally, some power plants might be located outside of NUTS1 regions (e.g. the ones in Ukraine). Hence, the total capacity in cap is lower than the total capacity in gdf.

Finally, we can plot the total hard coal generation capacity per NUTS1 region on a map. For regions without any coal power plants we can use the missing_kwds parameter to set a specific appearance (here red hatching).

This concludes the geopandas and cartopy tutorial!

Exercises¶

Take the following code snippet as a starting point. It loads the NUTS regions of Europe, the power plant dataset, and a shapefile for the Danish Natura2000 natural protection areas.

import cartopy
import pandas as pd
import geopandas as gpd
import cartopy.crs as ccrs
import matplotlib.pyplot as plt

url = "https://tubcloud.tu-berlin.de/s/RHZJrN8Dnfn26nr/download/NUTS_RG_10M_2021_4326.geojson"
nuts = gpd.read_file(url)

fn = "https://raw.githubusercontent.com/PyPSA/powerplantmatching/master/powerplants.csv"
df = pd.read_csv(fn, index_col=0)
geometry = gpd.points_from_xy(df["lon"], df["lat"])
ppl = gpd.GeoDataFrame(df, geometry=geometry, crs=4326)

url = "https://tubcloud.tu-berlin.de/s/mEpdmgBtmMbyjAr/download/Natura2000_end2021-DK.gpkg"
natura = gpd.read_file(url)

Task 1: Identify the coordinate reference system of the natura GeoDataFrame.

Task 2: Plot the natura GeoDataFrame on a map without transforming its CRS. Use cartopy for setting the projection of the figure and add coastlines and borders.

Task 3: Identify the name of the largest protected area in the natura GeoDataFrame.

Task 4: What is the total protection area in square kilometers.

Task 5: The natura GeoDataFrame has a column SITETYPE that indicates the type of protected area. Calculate the total area for each type of protected area (again in square kilometers).

Task 6: By how much (in percent) would the total area of protected areas increase if a buffer of 1 km around each protected area were also protected?

Task 7: List the power plants that are located within protected areas. How many power plants are located within protected areas? Use the .sjoin() function. Check the result by plotting these power plants on top of the protected areas.

Task 8 (advanced): What fraction of the natural protection area is located offshore? Use set operations with the .overlay() function and the NUTS regions GeoDataFrame.