In [1]:
%matplotlib inline
import inspect, sys
In [2]:
# check pydov path
import pydov
In [3]:
from pydov.search.grondwaterfilter import GrondwaterFilterSearch
gwfilter = GrondwaterFilterSearch()
A description is provided for the 'GrondwaterFilter' datatype:
In [4]:
print(gwfilter.get_description())
The different fields that are available for objects of the 'GrondwaterFilter' datatype can be requested with the get_fields() method:
In [5]:
fields = gwfilter.get_fields()
# print available fields
for f in fields.values():
print(f['name'])
You can get more information of a field by requesting it from the fields dictionary:
In [6]:
# print information for a certain field
fields['aquifer']
Out[6]:
Optionally, if the values of the field have a specific domain the possible values are listed as values:
In [7]:
# if an attribute can have several values, these are listed under 'values', e.g. for 'putsoort':
fields['putsoort']['values']
Out[7]:
Get data for all the groundwater screens that are geographically located within the bounds of the specified box.
The coordinates are in the Belgian Lambert72 (EPSG:31370) coordinate system and are given in the order of lower left x, lower left y, upper right x, upper right y.
In [8]:
from pydov.util.location import Within, Box
df = gwfilter.search(location=Within(Box(93378, 168009, 94246, 169873)))
df.head()
Out[8]:
Using the pkey attributes one can request the details of the corresponding put or filter in a webbrowser:
In [9]:
for pkey_grondwaterlocatie in set(df.pkey_grondwaterlocatie):
print(pkey_grondwaterlocatie)
for pkey_filter in set(df.pkey_filter):
print(pkey_filter)
Next to querying groundwater screens based on their geographic location within a bounding box, we can also search for groundwater screens matching a specific set of properties. For this we can build a query using a combination of the 'GrondwaterFilter' fields and operators provided by the WFS protocol.
A list of possible operators can be found below:
In [10]:
[i for i,j in inspect.getmembers(sys.modules['owslib.fes'], inspect.isclass) if 'Property' in i]
Out[10]:
In this example we build a query using the PropertyIsEqualTo operator to find all groundwater screens that are within the community (gemeente) of 'Hamme':
In [11]:
from owslib.fes import PropertyIsEqualTo
query = PropertyIsEqualTo(
propertyname='gemeente',
literal='Herstappe')
df = gwfilter.search(query=query)
df.head()
Out[11]:
Once again we can use the pkey_filter as a permanent link to the information of the groundwater screens:
In [12]:
for pkey_filter in set(df.pkey_filter):
print(pkey_filter)
In [13]:
query = PropertyIsEqualTo(propertyname='gemeente',
literal='Gent')
df = gwfilter.search(query=query,
return_fields=('pkey_filter', 'x', 'y', 'meetnet'))
df.head()
Out[13]:
In [14]:
query = PropertyIsEqualTo(propertyname='gemeente',
literal='Boortmeerbeek')
df = gwfilter.search(query=query,
return_fields=('pkey_filter', 'meetnet', 'meetnet_code'))
df.head()
Out[14]:
In [15]:
from owslib.fes import PropertyIsLike
query = PropertyIsLike(propertyname='meetnet',
literal='meetnet 9 %')
df = gwfilter.search(query=query,
location=Within(Box(87676, 163442, 91194, 168043)))
df.head()
Out[15]:
Get all groundwater screens in Hamme that have a value for length_filter and either belong to the primary meetnet of VMM or that have a depth bottom screen less than 3 meter.
In [16]:
from owslib.fes import Or, Not, PropertyIsNull, PropertyIsLessThanOrEqualTo, And
query = And([PropertyIsEqualTo(propertyname='gemeente',
literal='Hamme'),
Not([PropertyIsNull(propertyname='lengte_filter')]),
Or([PropertyIsLike(propertyname='meetnet',
literal='meetnet 1%'),
PropertyIsLessThanOrEqualTo(
propertyname='diepte_onderkant_filter',
literal='3')])])
df_hamme = gwfilter.search(query=query,
return_fields=('pkey_filter', 'x', 'y', 'gw_id', 'filternummer', 'diepte_onderkant_filter'))
df_hamme.head()
Out[16]:
Get some data and keep only the relevant waterlevel measurements:
In [17]:
query = PropertyIsEqualTo(
propertyname='pkey_filter',
literal='https://www.dov.vlaanderen.be/data/filter/1997-000494')
df = gwfilter.search(query=query)
df = df[df.filterstatus == "in rust"]
If the tijdstip field contains data, it can be combined with the datum field to create a date.datetime object. Make sure to fill the empty tijdstip fields with a default timestamp.
In [18]:
import pandas as pd
df.reset_index(inplace=True)
df['tijdstip'] = df.tijdstip.fillna('00:00:00')
df['tijd'] = pd.to_datetime(df.datum.astype(str)+' '+df.tijdstip.astype(str))
df.tijd.head()
Out[18]:
For further analysis and visualisation of the time series data, we can use the data analysis library pandas and visualisation library matplotlib.
In [19]:
import pandas as pd
import matplotlib.pyplot as plt
Query the data of a specific filter using its pkey:
In [20]:
query = PropertyIsEqualTo(
propertyname='pkey_filter',
literal='https://www.dov.vlaanderen.be/data/filter/2003-009883')
df = gwfilter.search(query=query)
df.head()
Out[20]:
The date is still stored as a string type. Transforming to a data type using the available pandas function to_datetime and using these dates as row index:
In [21]:
df['datum'] = pd.to_datetime(df['datum'])
df = df.set_index('datum')
The default plotting functionality of Pandas can be used:
In [28]:
df.head()
Out[28]:
In [37]:
ax = df['peil_mtaw'].plot(style='-', figsize=(12, 5))
ax.set_title('Water heads `Put ZWAP205` in Houthalen-Helchteren');
ax.set_ylabel('head (m TAW)');
ax.set_xlabel('');
Or a combination with matplotlib to have full customization options:
In [25]:
from matplotlib.dates import MonthLocator, YearLocator, DateFormatter
from matplotlib.ticker import MaxNLocator, MultipleLocator
# Get height of ground surface
ground_surface = df["mv_mtaw"][0]
# create a plot with 2 subplots
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 6),
sharex=False, sharey=True)
# Plot entire time series in the upper plot
df['peil_mtaw'].plot(ax=ax1, title='Water heads `Put ZWAP205`')
ax1.xaxis.set_major_locator(YearLocator(5))
ax1.xaxis.set_major_formatter(DateFormatter('%Y'))
# Plot the data for 2011 in the lower plot
df['peil_mtaw']["2011"].plot(ax=ax2, title='Water heads `Put ZWAP205` year 2011')
ax2.xaxis.set_major_locator(MonthLocator(interval=3))
ax2.xaxis.set_major_formatter(DateFormatter('%Y-%m'))
# Adjust configuration of plot
for ax in (ax1, ax2):
ax.set_xlabel('')
ax.set_ylabel('head (mTAW)')
for tick in ax.get_xticklabels():
tick.set_rotation(0)
tick.set_horizontalalignment('center')
# Only draw spine between the y-ticks
ax.spines['left'].set_position(('outward', 10))
# Hide the right and top spines
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
ax.yaxis.set_major_locator(MultipleLocator(0.2))
# Add the ground surface (provided in the data) on the subplots
ax.axhline(ground_surface, color = 'brown')
ax.annotate('Ground surface',
xy=(0.05, 0.68),
xycoords='axes fraction',
xytext=(-25, -15), textcoords='offset points',
fontsize=12, color='brown')
fig.tight_layout(h_pad=5)
The Pandas package provides a the functionality to further analyze and process time series data. Particularly, the resample function can be useful.
For example, calculate the yearly minima and maxima of the time series:
In [38]:
df["peil_mtaw"].resample("A").agg(['min', 'max'])
Out[38]:
or the monthly minima and maxima:
In [40]:
ax = df["peil_mtaw"].resample("M").agg(['min', 'max']).plot()
ax.set_title('Monthly minimum and maximum water heads `Put ZWAP205`');
ax.set_ylabel('head (m TAW)');
ax.set_xlabel('');
Calculate 10 and 90 percentiles of the time series with the quantile function:
In [46]:
fig, ax = plt.subplots(figsize=(12, 4))
ax.plot(df["peil_mtaw"], label='head (mTAW)', linewidth=0.75)
ax.axhline(df["peil_mtaw"].quantile(0.1), color = 'brown', label='p10', linewidth=3)
ax.axhline(df["peil_mtaw"].quantile(0.9), color = 'darkblue', label='p90', linewidth=3)
handles, labels = ax.get_legend_handles_labels()
ax.set_title('Water heads `Put ZWAP205` with 10% and 90% quantile');
ax.set_ylabel('head (m TAW)');
ax.set_xlabel('');
fig.legend(handles, labels)
Out[46]:
A duration exceedance curve provides the percentage of time that the water head is above a given value:
In [47]:
import numpy as np
In [50]:
sorted_heads = np.sort(df["peil_mtaw"])[::-1]
exceedence = np.arange(1.,len(sorted_heads)+1) / len(sorted_heads)
fig, ax = plt.subplots()
ax.plot(exceedence*100, sorted_heads)
ax.set_xlabel("Exceedence [%]");
ax.set_ylabel("Water head (mTAW)");
ax.set_title("Duration exceedance curve `Put ZWAP205`");
Using Folium, we can display the results of our search on a map.
In [51]:
# import the necessary modules (not included in the requirements of pydov!)
import folium
from folium.plugins import MarkerCluster
from pyproj import Transformer
In [52]:
# convert the coordinates to lat/lon for folium
def convert_latlon(x1, y1):
transformer = Transformer.from_crs("epsg:31370", "epsg:4326", always_xy=True)
x2,y2 = transformer.transform(x1, y1)
return x2, y2
df_hamme['lon'], df_hamme['lat'] = zip(*map(convert_latlon, df_hamme['x'], df_hamme['y']))
# convert to list
loclist = df_hamme[['lat', 'lon']].values.tolist()
In [53]:
# initialize the Folium map on the centre of the selected locations, play with the zoom until ok
fmap = folium.Map(location=[df_hamme['lat'].mean(), df_hamme['lon'].mean()], zoom_start=12)
marker_cluster = MarkerCluster().add_to(fmap)
for loc in range(0, len(loclist)):
folium.Marker(loclist[loc], popup=df_hamme['gw_id'][loc]).add_to(marker_cluster)
fmap
Out[53]:
In [54]:
gwfilter = GrondwaterFilterSearch()
In [55]:
query = PropertyIsLike(
propertyname='grondwaterlichaam',
literal='BLKS_1000%')
df = gwfilter.search(max_features=10,
query=query,
return_fields=("pkey_filter", "gw_id", "filternummer", "x", "y", "grondwaterlichaam_code"))
df
Out[55]: