Geopandas and Geoplot¶

Below is a simple example of using geopandas and geoplot to visualize some spatial data. The data is derived from SMHI Open data and is downloaded and stored in GeoJSON by a Python script using Python bindings and utilities for GeoJSON.

Import libraries and ignore annoying warning message

In [7]:

%matplotlib inline

import matplotlib.pyplot as plt
import geopandas as gpd
import geoplot as gplt
import warnings; warnings.filterwarnings('ignore', 'GeoSeries.isna', UserWarning)

Read the shape file downloaded from Natural Earth. The file used is Admin 1 - States, Provinces. Then filter out Sweden.

In [8]:

world = gpd.read_file("data/ne_10m_admin_1_states_provinces.shp")
swe = world[world['admin'] == 'Sweden']  # Filter out Sweden from the world

Now read the data file on average daily temperatures from SMHI, use the coordinate reference system WGS84 ("epsg=4326") and display some data.

In [9]:

avg_temp = gpd.read_file("https://www.viltstigen.se/metobs/latest/02*") # Probl in EPSG 3006 or 3021
avg_temp = avg_temp.to_crs(epsg=4326)
avg_temp.head()

Out[9]:

	id	key	title	summary	updated	timestamp	height	value	geometry
0	Abisko Aut	02	Lufttemperatur	medelvärde 1 dygn, 1 gång/dygn, kl 00	1593896400000	2020-07-04T21:00:00	392.303	8.5	POINT (18.81640 68.35380)
1	Adelsö A	02	Lufttemperatur	medelvärde 1 dygn, 1 gång/dygn, kl 00	1593896400000	2020-07-04T21:00:00	5.612	16.7	POINT (17.52130 59.35790)
2	Arjeplog A	02	Lufttemperatur	medelvärde 1 dygn, 1 gång/dygn, kl 00	1593896400000	2020-07-04T21:00:00	430.839	7.6	POINT (17.83960 66.05130)
3	Arvidsjaur A	02	Lufttemperatur	medelvärde 1 dygn, 1 gång/dygn, kl 00	1593896400000	2020-07-04T21:00:00	382.450	6.2	POINT (19.26420 65.59410)
4	Arvika A	02	Lufttemperatur	medelvärde 1 dygn, 1 gång/dygn, kl 00	1593896400000	2020-07-04T21:00:00	65.758	12.5	POINT (12.63540 59.67430)

Equally, read some pressure data from the same date and display.

In [10]:

pressure = gpd.read_file("https://www.viltstigen.se/metobs/latest/09*")
pressure = pressure.to_crs(epsg=4326)
pressure.head()

Out[10]:

	id	key	title	summary	updated	timestamp	height	value	geometry
0	Abisko Aut	09	Lufttryck reducerat havsytans nivå	vid havsytans nivå, momentanvärde, 1 gång/tim	1593896400000	2020-07-04T21:00:00	392.303	1000.4	POINT (18.81640 68.35380)
1	Arjeplog A	09	Lufttryck reducerat havsytans nivå	vid havsytans nivå, momentanvärde, 1 gång/tim	1593896400000	2020-07-04T21:00:00	430.839	998.7	POINT (17.83960 66.05130)
2	Arvidsjaur A	09	Lufttryck reducerat havsytans nivå	vid havsytans nivå, momentanvärde, 1 gång/tim	1593896400000	2020-07-04T21:00:00	382.450	997.0	POINT (19.26420 65.59410)
3	Arvika A	09	Lufttryck reducerat havsytans nivå	vid havsytans nivå, momentanvärde, 1 gång/tim	1593896400000	2020-07-04T21:00:00	65.758	1001.9	POINT (12.63540 59.67430)
4	Axstål Mo	09	Lufttryck reducerat havsytans nivå	vid havsytans nivå, momentanvärde, 1 gång/tim	1593896400000	2020-07-04T21:00:00	91.000	1002.5	POINT (14.57220 58.57190)

Now do the plotting. Setup the projection, conic projection Albers. 2 plots will be made, for temperature and pressure, hence 2 axis is created.

Do a voronoi diagram for temperatures and pressure, then plot Sweden (including counties) on top. For the voronoi plts, use the colormap ("cmap") coolwarm, select the value as 'hue="value"`.

In [11]:

proj = gplt.crs.AlbersEqualArea()
fig = plt.figure(figsize=(10,5))
ax1 = plt.subplot(121, projection=proj)
ax2 = plt.subplot(122, projection=proj)
gplt.voronoi(avg_temp, 
             cmap='coolwarm', 
             clip=swe, 
             hue="value", 
             legend=True, 
             edgecolor="None", 
             projection=proj, 
             ax=ax1)
gplt.polyplot(swe, edgecolor="Black", zorder=1, linewidth=0.5, projection=proj, ax=ax1)
gplt.voronoi(pressure, 
             cmap='coolwarm', 
             clip=swe, 
             hue="value", 
             legend=True, 
             edgecolor="None", 
             projection=proj, 
             ax=ax2)
gplt.polyplot(swe, edgecolor="Black", zorder=1, linewidth=0.5, projection=proj, ax=ax2)

Out[11]:

<cartopy.mpl.geoaxes.GeoAxesSubplot at 0x7fd4f4a5cda0>

That's it