Exploring Groningen Earthquake data in Python

Stefan Crooijmans
4 min readJan 5, 2022


Heat map of earthquake occurrences since January 2000. Polygons include Groningen province (black) and gas/oil fields (grey). Created with folium.


In the mid 1990’s, the Royal Netherlands Meteorological Institute (KNMI) constructed a network of seismometers to measure the occurrences of earthquakes in Groningen. Earthquake data is monitored by the KNMI and can be scraped through API requests. Based on the longitude and latitude of the data we can filter the data within the geometry of a shapefile, in this case the province of Groningen.

from shapely.geometry import Pointprovinces = 'data/provincies/B1_Provinciegrenzen_van_NederlandPolygon.shp'
provinces = gpd.read_file(provinces)
provinces.to_crs(epsg=4326, inplace=True)
df = pd.read_csv('output_df_earthquakes.csv')
df['within'] = ""
gdf = GeoDataFrame(df, geometry=geometry)
gdf.to_crs(epsg=4326, inplace = True)
groningen = provinces[provinces.PROV_NAAM=='Groningen']
within_list = []
for lon,lat in zip(df.lon, df.lat):
pt = Point(lon, lat)
groningen = groningen.geometry
within = pt.within(groningen.geometry.values[0]) # this yields values within the province
df['within'] = within_list
df_within = df.loc[df.within ==True]

Github Repositories

The code used to generate the graphics in this article can be found in this GitHub repositories below. Packages used for this project include selenium, json, pandas, seaborn, folium, and geopandas.

Declining Production of the Groningen Gas Field

The Groningen gas field was discovered in 1957. Since then, hundreds of wells have been drilled to enhance production from the field. Production ramped up to 80–90 billion Nm³ per annum in 1970’s and became responsible for nearly 7% of the GDP for the next decade.

Gas from the Groningen gas field is produced from a reservoir at a depth of approximately 3km. These reservoir contains compressed sand-grains which make up the matrix of the rock. Fluid filled pore space makes up the remainder of the rock. Fluid (gas) depletion causes the sponge to compress which in turn causes subsidence and differential stress along faults. The stress across faults builds and when movement occurs energy is released resulting in the occurrence of an earthquake.

Have Earthquakes decreased in magnitude with decreasing gas production?

Reservoir depletion combined with increasing political pressure due to damage of surrounding houses have caused gas production to steadily decline at an annual rate of 10–15% from 52 billion Nm³ in 2013. Politicians have decided to cease production from the largest gas field in the Netherlands in the coming year(s).

Significant increase in earthquakes in Groningen since the early 2000's. Data highlights a disproportional increase in lower ’magnitude earthquakes since the beginning of the 21st century. See code below.
import seaborn as sns
import numpy as np
import matplotlib.pyplot as plt
def add_mean_line(data, var=None, **kws):
# Derive mean for each group
m = np.mean(data[var])
nums = np.array(data[var])
# get axis
ax = plt.gca()
#add line at group_mean
ax.axvline(m, color='black', lw=2, ls='--')
# Annotate Group Mean
x_pos = 0.5
ax.text(x_pos, .6, f'mean={m:.2f}' '\n' f'count={nums.size:.0f}',
transform = ax.transAxes, fontsize=12)
bins = [0,0.2,0.4,0.6,0.8,1.0,1.2,1.4,1.6,1.8,2,2.2,2.4,2.6,2.8,3.0,3.2,3.4,3.6,3.8]
var ='mag'
g = sns.FacetGrid(df, col='cat_y',col_wrap=2, sharex=False, xlim=(0,3.5))
g.map_dataframe(sns.histplot, x= 'mag', edgecolor="black", bins=bins)
g.map_dataframe(add_mean_line, var =var)
g.set_axis_labels('Magnitude', 'Frequency', fontsize=14)
g.set_titles(col_template = '{col_name}', size=14)
g.savefig('facetgrid.png', dpi=300)

The KNMI data highlights a twofold increase in earthquake occurrences in Groningen in the past decade relative to the beginning of the century. The number of earthquakes annually has increased to 100 per year in Groningen alone. More than 80% of the recorded earthquakes are below a magnitude of 2.0. Although it is possible that an increasing number of earthquakes occur, the KNMI explains that it has increased the density of its stations since 2014–15 and subsequently has become better at identifying lower magnitude earthquakes. The growth in earthquakes appears to have flattened in the past few years and may be linked to the sharp decline in gas production. The number of >3 magnitude earthquakes appears steady, although they do continue to occur. The big questions is whether the decreasing production will pay dividends and lead to a further reduction of both the frequency and magnitude of earthquakes. Such a result may prove significant in affecting future government policies in- and outside the Netherlands.

Spatial Distribution of Earthquakes across the Netherlands

def generate_heat_map(df, df2, date, colorGradient, colorGradient2):#Give a start location
location=[53.2194, 6.7563]
# Filter dataframe by dates
df = df.loc[df.dates<=date]
df2 = df2.loc[df2.dates<=date]
# Assign heat map
heat_map = folium.Map(location=location, tiles = 'CartoDB positron', zoom_start=9)
# save all lat and lon and rescaled magnitudes to list values
latlons = df2[['lat', 'lon']].values.tolist()
HeatMap(latlons, gradient=colorGradient2,min_opacity=.7, blur=1, radius=2).add_to(heat_map)
heat_map = add_provinces(provinces,heat_map)
heat_map = add_fields(fields, heat_map)
dates = df.dates.unique()
start_date = str(dates[0])[:10]
stop_date = str(dates[-1])[:10]
start = datetime.strptime(start_date, "%Y-%m-%d")
stop = datetime.strptime(stop_date, "%Y-%m-%d")
colorGradient = {0.0:'blue',0.25:'cyan',0.5:'lime',0.75:'yellow', 1.0:'red'}
colorGradient2 = {0.0:'grey',1.0:'grey'}
#loop through dates
while start < stop:
generate_heat_map(df, df2, start, colorGradient,colorGradient2)
start = start + relativedelta(months=+1)

The heat map shows the spatial distribution of of earthquakes since 2000 across the Groningen province. The earthquakes coincide with the outlines of the gas fields and appear to spread more towards the outer edges of the field (deeper parts of the gas reservoir) that get depleted with increasing production.



Stefan Crooijmans

Shell Geoscientist passionate about building products.