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:

[pastacode lang=”python” manual=”%23%20import%20packages%0Aimport%20geopandas%20as%20gpd%0Aimport%20matplotlib.pyplot%20as%20plt%0Afrom%20descartes%20import%20PolygonPatch%0Aimport%20shapefile%0Aimport%20cartopy.crs%20as%20ccrs%0Aimport%20cartopy.feature%20as%20cfeature%0Afrom%20matplotlib.offsetbox%20import%20AnchoredText” message=”” highlight=”” provider=”manual”/]

Step two

Load up your NUTS shapefile (obtain from eurostat)

[pastacode lang=”python” manual=”%23%20read%20the%20NUTS%20shapefile%20and%20extract%20the%20polygons%20for%20a%20individual%20countries%0Anuts%3Dgpd.read_file(‘%2Fpath%2Fto%2Fyour%2Fshapefiles%2FNUTS_RG_10M_2016_4326_LEVL_2.shp’)” message=”” highlight=”” provider=”manual”/]

Step three

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

[pastacode lang=”python” manual=”%23%20use%20representitive_point%20method%20to%20obtain%20point%20from%20within%20each%20shapefile%0A%23%20creates%20new%20column%20’coords’%20containing%20the%20representitive%20points%0Anuts%5B’coords’%5D%20%3D%20nuts%5B’geometry’%5D.apply(lambda%20x%3A%20x.representative_point().coords%5B%3A%5D)%0Anuts%5B’coords’%5D%20%3D%20%5Bcoords%5B0%5D%20for%20coords%20in%20nuts%5B’coords’%5D%5D” message=”” highlight=”” provider=”manual”/]

Step four

Using matplotlib, plot the map.

[pastacode lang=”python” manual=”sf%3Dshapefile.Reader(%22%2FUsers%2Fuser%2FDocuments%2FC3S_Energy_Docs%2FNUTS2%2Fref-nuts-2016-10m.shp%2FNUTS_RG_10M_2016_4326_LEVL_2.shp%22)%0Adef%20main()%3A%0A%20%20%20%20fig%20%3D%20plt.figure(figsize%3D(80%2C40))%20%0A%20%20%20%20%0A%20%20%20%20ax%20%3D%20fig.add_subplot(1%2C%201%2C%201%2C%20projection%3Dccrs.PlateCarree())%0A%20%20%20%20%0A%20%20%20%20for%20poly%20in%20sf.shapes()%3A%0A%20%20%20%20%20%20%20%20poly_geo%3Dpoly.__geo_interface__%0A%20%20%20%20%20%20%20%20ax.add_patch(PolygonPatch(poly_geo%2C%20fc%3D’%23ffffff’%2C%20ec%3D’%23000000’%2C%20alpha%3D0.5%2C%20zorder%3D2%20))%0A%20%20%20%20ax.axis(‘scaled’)%0A%0A%0Aif%20__name__%20%3D%3D%20’__main__’%3A%0A%20%20%20%20main()” message=”” highlight=”” provider=”manual”/]

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.

[pastacode lang=”python” manual=”sf%3Dshapefile.Reader(%22%2FUsers%2Fuser%2FDocuments%2FC3S_Energy_Docs%2FNUTS2%2Fref-nuts-2016-10m.shp%2FNUTS_RG_10M_2016_4326_LEVL_2.shp%22)%0Adef%20main()%3A%0A%20%20%20%20fig%20%3D%20plt.figure(figsize%3D(80%2C40))%20%0A%20%20%20%20%0A%20%20%20%20ax%20%3D%20fig.add_subplot(1%2C%201%2C%201%2C%20projection%3Dccrs.PlateCarree())%0A%20%20%20%20%0A%20%20%20%20for%20poly%20in%20sf.shapes()%3A%0A%20%20%20%20%20%20%20%20poly_geo%3Dpoly.__geo_interface__%0A%20%20%20%20%20%20%20%20ax.add_patch(PolygonPatch(poly_geo%2C%20fc%3D’%23ffffff’%2C%20ec%3D’%23000000’%2C%20alpha%3D0.5%2C%20zorder%3D2%20))%0A%20%20%20%20ax.axis(‘scaled’)%0A%20%20%20%20%0A%20%20%20%20%23%20adjust%20plot%20domain%20to%20focus%20on%20EU%20region%0A%20%20%20%20plt.xlim(-25%2C%2047)%0A%20%20%20%20plt.ylim(30%2C%2073)%20%20%20%0A%0Aif%20__name__%20%3D%3D%20’__main__’%3A%0A%20%20%20%20main()” message=”” highlight=”” provider=”manual”/]

Step six

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

[pastacode lang=”python” manual=”sf%3Dshapefile.Reader(%22%2FUsers%2Fuser%2FDocuments%2FC3S_Energy_Docs%2FNUTS2%2Fref-nuts-2016-10m.shp%2FNUTS_RG_10M_2016_4326_LEVL_2.shp%22)%0Adef%20main()%3A%0A%20%20%20%20fig%20%3D%20plt.figure(figsize%3D(80%2C40))%20%0A%20%20%20%20%0A%20%20%20%20ax%20%3D%20fig.add_subplot(1%2C%201%2C%201%2C%20projection%3Dccrs.PlateCarree())%0A%20%20%20%20%23hi%20res%20polygons%20for%20land%2C%20sea%20etc%0A%20%20%20%20ax.coastlines(resolution%3D’10m’%2C%20color%3D’black’%2C%20linewidth%3D1)%0A%20%20%20%20ax.natural_earth_shp(name%3D’land’%2C%20resolution%3D’10m’%2C%20category%3D’physical’)%0A%20%20%20%20ax.natural_earth_shp(name%3D’ocean’%2C%20resolution%3D’50m’%2C%20category%3D’physical’%2Cfacecolor%3D’lightblue’)%0A%20%20%20%20%0A%20%20%20%20for%20poly%20in%20sf.shapes()%3A%0A%20%20%20%20%20%20%20%20poly_geo%3Dpoly.__geo_interface__%0A%20%20%20%20%20%20%20%20ax.add_patch(PolygonPatch(poly_geo%2C%20fc%3D’%23ffffff’%2C%20ec%3D’%23000000’%2C%20alpha%3D0.5%2C%20zorder%3D2%20))%0A%20%20%20%20ax.axis(‘scaled’)%0A%20%20%20%20%0A%20%20%20%20%23%20adjust%20plot%20domain%20to%20focus%20on%20EU%20region%0A%20%20%20%20plt.xlim(-25%2C%2047)%0A%20%20%20%20plt.ylim(30%2C%2073)%20%20%20%0A%0Aif%20__name__%20%3D%3D%20’__main__’%3A%0A%20%20%20%20main()” message=”” highlight=”” provider=”manual”/]

Step seven

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

[pastacode lang=”python” manual=”sf%3Dshapefile.Reader(%22%2Fpath%2Fto%2Fyour%2Fshapefiles%2Fref-nuts-2016-10m.shp%2FNUTS_RG_10M_2016_4326_LEVL_2.shp%22)%0Adef%20main()%3A%0A%20%20%20%20fig%20%3D%20plt.figure(figsize%3D(80%2C40))%20%0A%20%20%20%20%0A%20%20%20%20ax%20%3D%20fig.add_subplot(1%2C%201%2C%201%2C%20projection%3Dccrs.PlateCarree())%0A%20%20%20%20%0A%20%20%20%20%23hi%20res%20polygons%20for%20land%2C%20sea%20etc%0A%20%20%20%20ax.coastlines(resolution%3D’10m’%2C%20color%3D’black’%2C%20linewidth%3D1)%0A%20%20%20%20ax.natural_earth_shp(name%3D’land’%2C%20resolution%3D’10m’%2C%20category%3D’physical’)%0A%20%20%20%20ax.natural_earth_shp(name%3D’ocean’%2C%20resolution%3D’50m’%2C%20category%3D’physical’%2Cfacecolor%3D’lightblue’)%0A%0A%20%20%20%20%0A%20%20%20%20for%20poly%20in%20sf.shapes()%3A%0A%20%20%20%20%20%20%20%20poly_geo%3Dpoly.__geo_interface__%0A%20%20%20%20%20%20%20%20ax.add_patch(PolygonPatch(poly_geo%2C%20fc%3D’%23ffffff’%2C%20ec%3D’%23000000’%2C%20alpha%3D0.5%2C%20zorder%3D2%20))%0A%20%20%20%20ax.axis(‘scaled’)%0A%0A%20%20%20%20%23%20annotate%20plot%20with%20NUTS_ID%20%0A%20%20%20%20for%20idx%2C%20row%20in%20nuts.iterrows()%3A%0A%20%20%20%20%20%20%20%20plt.annotate(s%3Drow%5B’NUTS_ID’%5D%2C%20xy%3Drow%5B’coords’%5D%2C%20horizontalalignment%3D’center’)%0A%0A%20%20%20%20%23%20adjust%20plot%20domain%20to%20focus%20on%20EU%20region%0A%20%20%20%20plt.xlim(-25%2C%2047)%0A%20%20%20%20plt.ylim(30%2C%2073)%20%20%20%20%20%0A%0A%20%20%20%20%23%20Add%20a%20text%20annotation%20for%20the%20license%20information%20to%20the%0A%20%20%20%20%23%20the%20bottom%20right%20corner.%0A%20%20%20%20text%20%3D%20AnchoredText(r’%24%5Cmathcircled%7B%7Bc%7D%7D%24%20Luke%20Sanger%20-%20WEMC%20(2020)’%2Cloc%3D4%2C%20prop%3D%7B’size’%3A%2012%7D%2C%20frameon%3DTrue)%0A%20%20%20%20ax.add_artist(text)%0A%20%20%20%20plt.show()%0A%0Aif%20__name__%20%3D%3D%20’__main__’%3A%0A%20%20%20%20main()” message=”” highlight=”” provider=”manual”/]

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)