Tech Blog 7: Plot a map with Eurostat NUTS2 regions in 7 easy steps

During our continued work on the C3S Energy project (a component of Copernicus Climate Change), we’ve been using Eurostat NUTS areas (nomenclature of territorial units for statistics) to define our European regions, both at country (NUTS0) and regional (NUTS2) levels.

Whilst there are plenty of resources available online for identifying the NUTS0 areas, there isn’t much out there detailing the NUTS2 regions and where they appear on a map.

After some research and experimentation at the behest of one of our colleagues, we were able to generate the resource here as a high resolution png file. Read on to find out how we did it.

How to show NUTS2 regions on a map

In this tutorial we will go though a workflow to produce a map with labelled regions. This should work for any geographic shapefiles.

Step one

Using Jupyter Notebooks, import the following packages:

# import packages
import geopandas as gpd
import matplotlib.pyplot as plt
from descartes import PolygonPatch
import shapefile
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from matplotlib.offsetbox import AnchoredText

Step two

Load up your NUTS shapefile (obtain from eurostat)

# read the NUTS shapefile and extract the polygons for a individual countries
nuts=gpd.read_file('/path/to/your/shapefiles/NUTS_RG_10M_2016_4326_LEVL_2.shp')

Step three

The next piece of code defines a central point within each shape for labelling use

# use representitive_point method to obtain point from within each shapefile
# creates new column 'coords' containing the representitive points
nuts['coords'] = nuts['geometry'].apply(lambda x: x.representative_point().coords[:])
nuts['coords'] = [coords[0] for coords in nuts['coords']]

Step four

Using matplotlib, plot the map.

sf=shapefile.Reader("/Users/user/Documents/C3S_Energy_Docs/NUTS2/ref-nuts-2016-10m.shp/NUTS_RG_10M_2016_4326_LEVL_2.shp")
def main():
    fig = plt.figure(figsize=(80,40)) 
    
    ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())
    
    for poly in sf.shapes():
        poly_geo=poly.__geo_interface__
        ax.add_patch(PolygonPatch(poly_geo, fc='#ffffff', ec='#000000', alpha=0.5, zorder=2 ))
    ax.axis('scaled')


if __name__ == '__main__':
    main()

And that’s it, we’ve plotted the NUTS2 regions!

But as you can see, it still needs some tweaking to make the regions clearer for the user.

Step five

Let’s first change the map x,y limits to focus on the EU domain.

sf=shapefile.Reader("/Users/user/Documents/C3S_Energy_Docs/NUTS2/ref-nuts-2016-10m.shp/NUTS_RG_10M_2016_4326_LEVL_2.shp")
def main():
    fig = plt.figure(figsize=(80,40)) 
    
    ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())
    
    for poly in sf.shapes():
        poly_geo=poly.__geo_interface__
        ax.add_patch(PolygonPatch(poly_geo, fc='#ffffff', ec='#000000', alpha=0.5, zorder=2 ))
    ax.axis('scaled')
    
    # adjust plot domain to focus on EU region
    plt.xlim(-25, 47)
    plt.ylim(30, 73)   

if __name__ == '__main__':
    main()

Step six

Then Matplotlib can call the open source natural earth map resources for use as map elements (such as land/sea areas etc).

sf=shapefile.Reader("/Users/user/Documents/C3S_Energy_Docs/NUTS2/ref-nuts-2016-10m.shp/NUTS_RG_10M_2016_4326_LEVL_2.shp")
def main():
    fig = plt.figure(figsize=(80,40)) 
    
    ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())
    #hi res polygons for land, sea etc
    ax.coastlines(resolution='10m', color='black', linewidth=1)
    ax.natural_earth_shp(name='land', resolution='10m', category='physical')
    ax.natural_earth_shp(name='ocean', resolution='50m', category='physical',facecolor='lightblue')
    
    for poly in sf.shapes():
        poly_geo=poly.__geo_interface__
        ax.add_patch(PolygonPatch(poly_geo, fc='#ffffff', ec='#000000', alpha=0.5, zorder=2 ))
    ax.axis('scaled')
    
    # adjust plot domain to focus on EU region
    plt.xlim(-25, 47)
    plt.ylim(30, 73)   

if __name__ == '__main__':
    main()

Step seven

Finally, looping through the central points in the shapefiles as defined earlier, gives us the labelling for each NUTS2 area.

sf=shapefile.Reader("/path/to/your/shapefiles/ref-nuts-2016-10m.shp/NUTS_RG_10M_2016_4326_LEVL_2.shp")
def main():
    fig = plt.figure(figsize=(80,40)) 
    
    ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())
    
    #hi res polygons for land, sea etc
    ax.coastlines(resolution='10m', color='black', linewidth=1)
    ax.natural_earth_shp(name='land', resolution='10m', category='physical')
    ax.natural_earth_shp(name='ocean', resolution='50m', category='physical',facecolor='lightblue')

    
    for poly in sf.shapes():
        poly_geo=poly.__geo_interface__
        ax.add_patch(PolygonPatch(poly_geo, fc='#ffffff', ec='#000000', alpha=0.5, zorder=2 ))
    ax.axis('scaled')

    # annotate plot with NUTS_ID 
    for idx, row in nuts.iterrows():
        plt.annotate(s=row['NUTS_ID'], xy=row['coords'], horizontalalignment='center')

    # adjust plot domain to focus on EU region
    plt.xlim(-25, 47)
    plt.ylim(30, 73)     

    # Add a text annotation for the license information to the
    # the bottom right corner.
    text = AnchoredText(r'$\mathcircled{{c}}$ Luke Sanger - WEMC (2020)',loc=4, prop={'size': 12}, frameon=True)
    ax.add_artist(text)
    plt.show()

if __name__ == '__main__':
    main()

There you have it.

We know that the London regions are quite densely clustered at this level and as a result, the labels are overlapping, so it’s not ideal. But this does demonstrate one way you can visualize NUTS2 Eurostat regions in the EU domain.

Luke Sanger (WEMC Data Engineer)

Recent Posts