%matplotlib inline
import inspect, sys
# check pydov path
import pydov
from pydov.search.interpretaties import InformeleStratigrafieSearch
itp = InformeleStratigrafieSearch()
A description is provided for the 'Informele stratigrafie' datatype:
itp.get_description()
'Een informele stratigrafie van een boring is, net als een formele stratigrafie, een lithostratigrafische interpretatie van een (gecodeerde) lithologische beschrijving per diepte-interval. Het verschil met de formele stratigrafie is dat de informele stratigrafie niet moet voldoen aan een algemeen aanvaarde standaard. De informele stratigrafie is een plaats in DOV om oude interpretaties of interpretaties die niet voldoen aan de huidige stratigrafische indeling toch weer te geven. De informele stratigrafie kan ook een gedetailleerde beschrijving bevatten van een lithostratigrafische eenheid.'
The different fields that are available for objects of the 'Informele stratigrafie' datatype can be requested with the get_fields() method:
fields = itp.get_fields()
# print available fields
for f in fields.values():
print(f['name'])
pkey_interpretatie Type_proef Proefnummer Proeffiche x y Z_mTAW diepte_tot_m gemeente Auteurs Datum Opdrachten betrouwbaarheid_interpretatie Geldig_van Geldig_tot diepte_laag_van diepte_laag_tot beschrijving pkey_boring pkey_sondering
You can get more information of a field by requesting it from the fields dictionary:
fields['Datum']
{'name': 'Datum', 'definition': 'De datum waarop de informele stratigrafie werd beschreven.', 'type': 'date', 'notnull': False, 'query': True, 'cost': 1}
Get data for all the 'informele stratigrafie' interpretations 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.
from pydov.util.location import Within, Box
df = itp.search(location=Within(Box(153145, 206930, 153150, 206935)))
df.head()
[000/001] .
pkey_interpretatie | pkey_boring | pkey_sondering | betrouwbaarheid_interpretatie | x | y | diepte_laag_van | diepte_laag_tot | beschrijving | |
---|---|---|---|---|---|---|---|---|---|
0 | https://www.dov.vlaanderen.be/data/interpretat... | https://www.dov.vlaanderen.be/data/boring/1974... | NaN | goed | 153147.0 | 206931.0 | 0.0 | 3.00 | Quartair |
1 | https://www.dov.vlaanderen.be/data/interpretat... | https://www.dov.vlaanderen.be/data/boring/1974... | NaN | goed | 153147.0 | 206931.0 | 3.0 | 14.05 | Tertiair - Mioceen: Antwerpiaan |
The dataframe contains one 'informele stratigrafie' interpretation where three layers ('laag') were identified. The available data are flattened to represent unique attributes per row of the dataframe.
Using the pkey_interpretatie field one can request the details of this interpretation in a webbrowser:
for pkey_interpretatie in set(df.pkey_interpretatie):
print(pkey_interpretatie)
https://www.dov.vlaanderen.be/data/interpretatie/1975-213484
Next to querying interpretations based on their geographic location within a bounding box, we can also search for interpretations matching a specific set of properties. For this we can build a query using a combination of the 'InformeleStratigrafie' fields and operators provided by the WFS protocol.
A list of possible operators can be found below:
[i for i,j in inspect.getmembers(sys.modules['owslib.fes'], inspect.isclass) if 'Property' in i]
['PropertyIsBetween', 'PropertyIsEqualTo', 'PropertyIsGreaterThan', 'PropertyIsGreaterThanOrEqualTo', 'PropertyIsLessThan', 'PropertyIsLessThanOrEqualTo', 'PropertyIsLike', 'PropertyIsNotEqualTo', 'PropertyIsNull', 'SortProperty']
In this example we build a query using the PropertyIsEqualTo operator to find all interpretations that are within the community (gemeente) of 'Herstappe':
from owslib.fes import PropertyIsEqualTo
query = PropertyIsEqualTo(propertyname='gemeente',
literal='Herstappe')
df = itp.search(query=query)
df.head()
[000/003] ...
pkey_interpretatie | pkey_boring | pkey_sondering | betrouwbaarheid_interpretatie | x | y | diepte_laag_van | diepte_laag_tot | beschrijving | |
---|---|---|---|---|---|---|---|---|---|
0 | https://www.dov.vlaanderen.be/data/interpretat... | https://www.dov.vlaanderen.be/data/boring/2016... | NaN | goed | 224687.4 | 158191.0 | 0.00 | 3.75 | q2 |
1 | https://www.dov.vlaanderen.be/data/interpretat... | https://www.dov.vlaanderen.be/data/boring/2016... | NaN | goed | 224687.4 | 158191.0 | 3.75 | 6.00 | q1c |
2 | https://www.dov.vlaanderen.be/data/interpretat... | https://www.dov.vlaanderen.be/data/boring/1993... | NaN | goed | 224820.0 | 157794.0 | 0.00 | 50.00 | ? |
3 | https://www.dov.vlaanderen.be/data/interpretat... | https://www.dov.vlaanderen.be/data/boring/1993... | NaN | goed | 224820.0 | 157794.0 | 5.00 | 50.00 | Formatie van Gulpen |
Once again we can use the pkey_interpretatie as a permanent link to the information of these interpretations:
for pkey_interpretatie in set(df.pkey_interpretatie):
print(pkey_interpretatie)
https://www.dov.vlaanderen.be/data/interpretatie/2010-106477 https://www.dov.vlaanderen.be/data/interpretatie/2016-271517 https://www.dov.vlaanderen.be/data/interpretatie/1993-097697
We can combine a query on attributes with a query on geographic location to get the interpretations within a bounding box that have specific properties.
The following example requests the interpretations of boreholes only, 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.)
from owslib.fes import PropertyIsEqualTo
query = PropertyIsEqualTo(
propertyname='Type_proef',
literal='Boring')
df = itp.search(
location=Within(Box(205000, 205000, 206000, 206000)),
query=query
)
df.head()
[000/003] ...
pkey_interpretatie | pkey_boring | pkey_sondering | betrouwbaarheid_interpretatie | x | y | diepte_laag_van | diepte_laag_tot | beschrijving | |
---|---|---|---|---|---|---|---|---|---|
0 | https://www.dov.vlaanderen.be/data/interpretat... | https://www.dov.vlaanderen.be/data/boring/1975... | NaN | goed | 205264.2 | 205607.5 | 0.0 | 8.0 | Kwartair |
1 | https://www.dov.vlaanderen.be/data/interpretat... | https://www.dov.vlaanderen.be/data/boring/1896... | NaN | goed | 205970.0 | 205270.5 | 0.0 | 2.3 | q4 |
2 | https://www.dov.vlaanderen.be/data/interpretat... | https://www.dov.vlaanderen.be/data/boring/1896... | NaN | goed | 205952.9 | 205993.0 | 0.0 | 1.0 | R |
3 | https://www.dov.vlaanderen.be/data/interpretat... | https://www.dov.vlaanderen.be/data/boring/1896... | NaN | goed | 205952.9 | 205993.0 | 1.0 | 2.3 | q4 |
We can look at one of the interpretations in a webbrowser using its pkey_interpretatie:
for pkey_interpretatie in set(df.pkey_interpretatie):
print(pkey_interpretatie)
https://www.dov.vlaanderen.be/data/interpretatie/1896-044150 https://www.dov.vlaanderen.be/data/interpretatie/1975-044576 https://www.dov.vlaanderen.be/data/interpretatie/1896-044149
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 'informele stratigrafie' interpretations in the city of Ghent and return their date:
query = PropertyIsEqualTo(propertyname='gemeente',
literal='Gent')
df = itp.search(query=query,
return_fields=('Datum',))
df.head()
Datum | |
---|---|
0 | NaN |
1 | NaN |
2 | NaN |
3 | 1981-02-25 |
4 | NaN |
df.describe()
Datum | |
---|---|
count | 642 |
unique | 248 |
top | 1996-11-25 |
freq | 41 |
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 interpretations as illustrated below.
For example, make a selection of the interpretations in municipality the of Antwerp, before 1/1/1900:
from owslib.fes import And, PropertyIsEqualTo, PropertyIsLessThan
query = And([PropertyIsEqualTo(propertyname='gemeente',
literal='Antwerpen'),
PropertyIsLessThan(propertyname='Datum',
literal='1900-01-01')]
)
df = itp.search(query=query,
return_fields=('pkey_interpretatie', 'Datum'))
df.head()
pkey_interpretatie | Datum | |
---|---|---|
0 | https://www.dov.vlaanderen.be/data/interpretat... | 1895-05-06 |
1 | https://www.dov.vlaanderen.be/data/interpretat... | 1895-05-10 |
2 | https://www.dov.vlaanderen.be/data/interpretat... | 1879-08-09 |
3 | https://www.dov.vlaanderen.be/data/interpretat... | 1879-08-18 |
4 | https://www.dov.vlaanderen.be/data/interpretat... | 1879-03-05 |
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:
query = PropertyIsEqualTo(
propertyname='gemeente',
literal='Herstappe')
df = itp.search(query=query,
return_fields=('pkey_interpretatie', 'pkey_boring', 'pkey_sondering',
'x', 'y', 'Z_mTAW', 'gemeente', 'Auteurs', 'Proefnummer'))
df.head()
pkey_interpretatie | pkey_boring | pkey_sondering | x | y | Z_mTAW | gemeente | Auteurs | Proefnummer | |
---|---|---|---|---|---|---|---|---|---|
0 | https://www.dov.vlaanderen.be/data/interpretat... | https://www.dov.vlaanderen.be/data/boring/2016... | NaN | 224687.4 | 158191.0 | 124.0 | Herstappe | Rutôt, A. - Belgische Geologische Dienst (BGD) | kb33d106e-B236 |
1 | https://www.dov.vlaanderen.be/data/interpretat... | https://www.dov.vlaanderen.be/data/boring/1993... | NaN | 224820.0 | 157794.0 | 132.0 | Herstappe | boormeester - bedrijf-dienst onbekend | kb41d120e-B1027 |
2 | https://www.dov.vlaanderen.be/data/interpretat... | https://www.dov.vlaanderen.be/data/boring/1993... | NaN | 224820.0 | 157794.0 | 132.0 | Herstappe | Lagrou, David - VLAKO | kb41d120e-B1027 |
Using Folium, we can display the results of our search on a map.
# import the necessary modules (not included in the requirements of pydov!)
import folium
from folium.plugins import MarkerCluster
from pyproj import Transformer
# 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()
# 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=12)
marker_cluster = MarkerCluster().add_to(fmap)
for loc in range(0, len(loclist)):
folium.Marker(loclist[loc], popup=df['Proefnummer'][loc]).add_to(marker_cluster)
fmap