In [1]:
%matplotlib inline
import inspect, sys
In [3]:
import pydov
In [4]:
from pydov.search.sondering import SonderingSearch
sondering = SonderingSearch()
A description is provided for the 'Sondering' datatype:
In [5]:
sondering.get_description()
Out[5]:
The different fields that are available for objects of the 'Sondering' datatype can be requested with the get_fields() method:
In [6]:
fields = sondering.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 [7]:
fields['diepte_sondering_tot']
Out[7]:
Optionally, if the values of the field have a specific domain the possible values are listed as values:
In [8]:
fields['conus']['values']
Out[8]:
Get data for all the CPT measurements 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 [9]:
from pydov.util.location import Within, Box
df = sondering.search(location=Within(Box(152999, 206930, 153050, 207935)))
df.head()
Out[9]:
The dataframe contains one CPT measurement where multiple measurement points. The available data are flattened to represent unique attributes per row of the dataframe.
Using the pkey_sondering field one can request the details of this borehole in a webbrowser:
In [10]:
for pkey_sondering in set(df.pkey_sondering):
print(pkey_sondering)
Next to querying CPT based on their geographic location within a bounding box, we can also search for CPT measurements matching a specific set of properties. For this we can build a query using a combination of the 'Sondering' fields and operators provided by the WFS protocol.
A list of possible operators can be found below:
In [11]:
[i for i,j in inspect.getmembers(sys.modules['owslib.fes'], inspect.isclass) if 'Property' in i]
Out[11]:
In this example we build a query using the PropertyIsEqualTo operator to find all CPT measuremetns that are within the community (gemeente) of 'Herstappe':
In [12]:
from owslib.fes import PropertyIsEqualTo
query = PropertyIsEqualTo(propertyname='gemeente',
literal='Elsene')
df = sondering.search(query=query)
df.head()
Out[12]:
Once again we can use the pkey_sondering as a permanent link to the information of these CPT measurements:
In [13]:
for pkey_sondering in set(df.pkey_sondering):
print(pkey_sondering)
We can combine a query on attributes with a query on geographic location to get the CPT measurements within a bounding box that have specific properties.
The following example requests the CPT measurements with a depth greater than or equal to 2000 meters within the given bounding box.
(Note that the datatype of the literal parameter should be a string, regardless of the datatype of this field in the output dataframe.)
In [14]:
from owslib.fes import PropertyIsGreaterThanOrEqualTo
query = PropertyIsGreaterThanOrEqualTo(
propertyname='diepte_sondering_tot',
literal='20')
df = sondering.search(
location=Within(Box(200000, 211000, 205000, 214000)),
query=query
)
df.head()
Out[14]:
We can look at one of the CPT measurements in a webbrowser using its pkey_sondering:
In [15]:
for pkey_sondering in set(df.pkey_sondering):
print(pkey_sondering)
We can limit the columns in the output dataframe by specifying the return_fields parameter in our search.
In this example we query all the CPT measurements in the city of Ghent and return their depth:
In [16]:
query = PropertyIsEqualTo(propertyname='gemeente',
literal='Gent')
df = sondering.search(query=query,
return_fields=('diepte_sondering_tot',))
df.head()
Out[16]:
In [17]:
df.describe()
Out[17]:
In [20]:
ax = df.boxplot()
ax.set_title('Distribution depth CPT measurements in Ghent');
ax.set_ylabel("depth (m)")
Out[20]:
To keep the output dataframe size acceptable, not all available WFS fields are included in the standard output. However, one can use this information to select CPT measurements as illustrated below.
For example, make a selection of the CPT measurements in municipality the of Antwerp, using a conustype 'U':
In [21]:
from owslib.fes import And
query = And([PropertyIsEqualTo(propertyname='gemeente',
literal='Antwerpen'),
PropertyIsEqualTo(propertyname='conus',
literal='U')]
)
df = sondering.search(query=query,
return_fields=('pkey_sondering', 'sondeernummer', 'x', 'y', 'diepte_sondering_tot', 'datum_aanvang'))
df.head()
Out[21]:
As denoted in the previous example, not all available fields are available in the default output frame to keep its size limited. However, you can request any available field by including it in the return_fields parameter of the search:
In [22]:
query = And([PropertyIsEqualTo(propertyname='gemeente', literal='Gent'),
PropertyIsEqualTo(propertyname='conus', literal='U')])
df = sondering.search(query=query,
return_fields=('pkey_sondering', 'sondeernummer', 'diepte_sondering_tot',
'conus', 'x', 'y'))
df.head()
Out[22]:
In [23]:
df
Out[23]:
The data for the reporting of resistivity plots with the online application, see for example this report, is also accessible with the pydov package. Querying the data for this specific sondering:
In [24]:
query = PropertyIsEqualTo(propertyname='pkey_sondering',
literal='https://www.dov.vlaanderen.be/data/sondering/1993-001275')
df_sond = sondering.search(query=query)
df_sond.head()
Out[24]:
We have the depth (z) available, together with the measured values for each depth of the variables (in dutch):
qc: Opgemeten waarde van de conusweerstand, uitgedrukt in MPa.Qt: Opgemeten waarde van de totale weerstand, uitgedrukt in kN.fs: Opgemeten waarde van de plaatelijke kleefweerstand uitgedrukt in kPa.u: Opgemeten waarde van de porienwaterspanning, uitgedrukt in kPa.i: Opgemeten waarde van de inclinatie, uitgedrukt in graden.To recreate the resistivity plot, we also need the resistivity number (wrijvingsgetal rf), see DOV documentation.
Notice: $f_s$ is provide in kPa and $q_c$ in MPa.
Adding rf to the dataframe:
In [25]:
df_sond["rf"] = df_sond["fs"]/df_sond["qc"]/10
Recreate the resistivity plot:
In [26]:
import matplotlib.pyplot as plt
In [27]:
def make_patch_spines_invisible(ax):
ax.set_frame_on(True)
ax.patch.set_visible(False)
for sp in ax.spines.values():
sp.set_visible(False)
In [36]:
fig, ax0 = plt.subplots(figsize=(8, 12))
# Prepare the individual axis
ax_qc = ax0.twiny()
ax_fs = ax0.twiny()
ax_u = ax0.twiny()
ax_rf = ax0.twiny()
for i, ax in enumerate([ax_qc, ax_fs, ax_u]):
ax.spines["top"].set_position(("axes", 1+0.05*(i+1)))
make_patch_spines_invisible(ax)
ax.spines["top"].set_visible(True)
# Plot the data on the axis
df_sond.plot(x="rf", y="z", label="rf", ax=ax_rf, color='purple', legend=False)
df_sond.plot(x="qc", y="z", label="qc (MPa)", ax=ax_qc, color='black', legend=False)
df_sond.plot(x="fs", y="z", label="fs (kPa)", ax=ax_fs, color='green', legend=False)
df_sond.plot(x="u", y="z", label="u (kPa)", ax=ax_u, color='red',
legend=False, xlim=(-100, 300)) # ! 300 is hardocded here for the example
# styling and configuration
ax_rf.xaxis.label.set_color('purple')
ax_fs.xaxis.label.set_color('green')
ax_u.xaxis.label.set_color('red')
ax0.axes.set_visible(False)
ax_qc.axes.yaxis.set_visible(False)
ax_fs.axes.yaxis.set_visible(False)
for i, ax in enumerate([ax_rf, ax_qc, ax_fs, ax_u, ax0]):
ax.spines["right"].set_visible(False)
ax.spines["bottom"].set_visible(False)
ax.xaxis.label.set_fontsize(15)
ax.xaxis.set_label_coords(-0.05, 1+0.05*i)
ax.spines['left'].set_position(('outward', 10))
ax.spines['left'].set_bounds(0, 30)
ax_rf.set_xlim(0, 46)
ax_u.set_title("Resistivity plot CPT measurement GEO-93/023-SII-E", fontsize=12)
ax0.invert_yaxis()
ax_rf.invert_xaxis()
ax_u.set_ylabel("Depth(m)", fontsize=12)
fig.legend(loc='lower center', ncol=4)
fig.tight_layout()
Using Folium, we can display the results of our search on a map.
In [37]:
# import the necessary modules (not included in the requirements of pydov!)
import folium
from folium.plugins import MarkerCluster
from pyproj import Transformer
In [38]:
# 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['lon'], df['lat'] = zip(*map(convert_latlon, df['x'], df['y']))
# convert to list
loclist = df[['lat', 'lon']].values.tolist()
In [39]:
# initialize the Folium map on the centre of the selected locations, play with the zoom until ok
fmap = folium.Map(location=[df['lat'].mean(), df['lon'].mean()], zoom_start=11)
marker_cluster = MarkerCluster().add_to(fmap)
for loc in range(0, len(loclist)):
folium.Marker(loclist[loc], popup=df['sondeernummer'][loc]).add_to(marker_cluster)
fmap
Out[39]: