Examples of search queries based on external WFS geometries

Use cases

Filter data within a certain community using its geographic borders
Filter data within a geographic boundary of a feature in layer "Bekken"
Filter data within a geographic boundary of a feature in layer "Hydrogeologische homogene zones"



In [1]:

    
%matplotlib inline
import matplotlib.pyplot as plt

from owslib.etree import etree
from owslib.fes import PropertyIsEqualTo
from owslib.wfs import WebFeatureService

from pydov.search.boring import BoringSearch
from pydov.search.grondwaterfilter import GrondwaterFilterSearch

from pydov.util.location import (
    GmlFilter,
    Within,
)

from pydov.util.owsutil import get_namespace

# 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

# convert to bytes if necessary
def to_bytes(data):
    if isinstance(data, bytes):
        return data
    elif isinstance(data, str):
        return data.encode('utf8')

Filter data within a certain community using its geographic borders

This example will use a third party WFS service to retrieve to geometry of (the boundary of) a community and subsequently use this boundary as a spatial filter to query boreholes from the DOV service.

We have to start by defining the third party WFS service that contains the featuretype (layer) that contains the geometry we're interested in:



In [2]:

    
gemeentegrenzen = WebFeatureService(
    'https://geoservices.informatievlaanderen.be/overdrachtdiensten/VRBG/wfs',
    version='1.1.0')

In this case we know in advance which featuretype from the WFS service we need (VRBG:RefGem) but are unsure which field (attribute) we can use to get the community we want.

We can get the schema (i.e. all the available fields) of a layer to find a field to query on:



In [3]:

    
gemeentegrenzen.get_schema('VRBG:Refgem')['properties']









    Out[3]:





{'UIDN': 'decimal',
 'OIDN': 'decimal',
 'TERRID': 'decimal',
 'NISCODE': 'string',
 'NAAM': 'string',
 'DATPUBLBS': 'date',
 'NUMAC': 'string'}

We can now build a search query to get the feature (including geometry) of the community with the NAAM of Lievegem:



In [4]:

    
naam_filter = PropertyIsEqualTo(propertyname='NAAM', literal='Lievegem')



In [5]:

    
gemeente_poly = gemeentegrenzen.getfeature(
    typename='VRBG:Refgem',
    filter=etree.tostring(naam_filter.toXML()).decode("utf8")).read()

And use this GML feature inside a GmlFilter query to find all boreholes that are located Within the community borders:



In [6]:

    
bs = BoringSearch()
df = bs.search(
    location=GmlFilter(gemeente_poly, Within),
    return_fields=('pkey_boring', 'x', 'y'))

Finally we can show the result on the map:



In [7]:

    
df['lon'], df['lat'] = zip(*map(convert_latlon, df['x'], df['y']))
# convert to list
loclist = df[['lat', 'lon']].values.tolist()



In [8]:

    
# 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=10)
marker_cluster = MarkerCluster().add_to(fmap)
for loc in range(0, len(loclist)):
    folium.Marker(loclist[loc], popup=df['pkey_boring'][loc]).add_to(marker_cluster)
fmap









    Out[8]:




Make this Notebook Trusted to load map: File -> Trust Notebook

Filter data within a geographic boundary of a feature in layer "Bekken"

This example will use a third party WFS service to retrieve to geometry of a waterbody and subsequently use this boundary as a spatial filter to query groundwater screens from the DOV service.

We have to start by defining the third party WFS service that contains the featuretype (layer) that contains the geometry we're interested in:



In [9]:

    
bekkens = WebFeatureService(
    'https://geoservices.informatievlaanderen.be/overdrachtdiensten/VHAZones/wfs',
    version='1.1.0'
)

We can list all the available featuretypes (layers) in the service if we are unsure which one to use:



In [10]:

    
list(bekkens.contents)









    Out[10]:





['VHAZones:Bekken', 'VHAZones:Deelbekken', 'VHAZones:Vhazone']

We can also get the schema (i.e. all the available fields) of a layer to find a field to query on:



In [11]:

    
bekkens.get_schema('VHAZones:Bekken')['properties']









    Out[11]:





{'UIDN': 'decimal',
 'OIDN': 'decimal',
 'BEKNR': 'short',
 'BEKNAAM': 'string',
 'STRMGEB': 'string'}

We can get all the distinct values of a specific field:



In [12]:

    
namespace = get_namespace(bekkens, 'VHAZones:Bekken')

tree = etree.fromstring(to_bytes(bekkens.getfeature('VHAZones:Bekken', propertyname='BEKNAAM').read()))
set((i.text for i in tree.findall('.//{%s}BEKNAAM' % namespace)))









    Out[12]:





{'Bekken Brugse polders',
 'Bekken Gentse kanalen',
 'Beneden-Scheldebekken',
 'Boven-Scheldebekken',
 'Demerbekken',
 'Denderbekken',
 'Dijlebekken',
 'Ijzerbekken',
 'Leiebekken',
 'Maasbekken',
 'Netebekken'}

We can now build a search query to get the feature (including geometry) of the waterbody with the BEKNAAM of Bekken Brugse polders:



In [13]:

    
naam_filter = PropertyIsEqualTo(propertyname='BEKNAAM', literal='Bekken Brugse polders')



In [14]:

    
bekken_poly = bekkens.getfeature(
    typename='VHAZones:Bekken',
    filter=etree.tostring(naam_filter.toXML()).decode("utf8")).read()

And use this GML feature inside a GmlFilter query to find all groundwater screens that are located Within the waterbody:



In [15]:

    
filter_search = GrondwaterFilterSearch()
df = filter_search.search(
    max_features = 100,
    location=GmlFilter(bekken_poly, Within),
    return_fields=('pkey_filter', 'x', 'y')
)

Finally we can show the result on the map:



In [16]:

    
df['lon'], df['lat'] = zip(*map(convert_latlon, df['x'], df['y']))
# convert to list
loclist = df[['lat', 'lon']].values.tolist()



In [17]:

    
# 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=10)
marker_cluster = MarkerCluster().add_to(fmap)
for loc in range(0, len(loclist)):
    folium.Marker(loclist[loc], popup=df['pkey_filter'][loc]).add_to(marker_cluster)
fmap









    Out[17]:




Make this Notebook Trusted to load map: File -> Trust Notebook

Filter data within a geographic boundary of a feature in layer "HHZ"

This example will use the DOV WFS service to retrieve to geometry of a HHZ (hydrogeological homogenous area) and subsequently use this boundary as a spatial filter to query groundwater screens from the DOV service.

We have to start by defining the DOV WFS service:



In [18]:

    
hhz = WebFeatureService(
    'https://www.dov.vlaanderen.be/geoserver/wfs',
    version='1.1.0'
)

Next we can get the schema (i.e. all the available fields) of a layer to find a field to query on:



In [19]:

    
hhz.get_schema('gw_varia:hhz')['properties']









    Out[19]:





{'id': 'string', 'hhz_naam': 'string', 'opp_m2': 'long', 'hhz': 'string'}

And get all the distinct values that occur for the hhz_naam field:



In [20]:

    
namespace = get_namespace(hhz, 'gw_varia:hhz')

tree = etree.fromstring(to_bytes(hhz.getfeature('gw_varia:hhz', propertyname='hhz_naam').read()))
set((i.text for i in tree.findall('.//{%s}hhz_naam' % namespace)))









    Out[20]:





{'Complex van de Kempen',
 'Duingebieden',
 'Dun quartair dek boven Ieperiaan klei',
 'Dun quartair dek boven Paniseliaan klei',
 'Dun quartair dek boven de Bartoon klei',
 'Dun quartair dek boven de Rupel klei',
 'Formatie van Berchem',
 'Formatie van Bolderberg',
 'Formatie van Brasschaat (+Merksplas)',
 'Formatie van Diest',
 'Formatie van Diest in de heuvelstreken',
 'Formatie van Eigenbilzen',
 'Formatie van Kasterlee',
 'Formatie van Kattendijk',
 'Formatie van Lillo en Poederlee',
 'Formatie van Mol',
 'Heersiaan',
 'Hoogterras-afzettingen',
 'Krijtafzettingen',
 'Landeniaan',
 'Ledo-Paniseliaan',
 'Ledo-Paniseliaan in de heuvelstreken',
 'Maas-Rijn-afzettingen',
 'Onder-Oligoceen (Tongeren+Bilzen)',
 'Paleoceen',
 'Polders (verzilte gebieden)',
 'Vlaamse Vallei (+bijrivieren en kustvlakte)',
 'Zanden van Brussel',
 'Zanden van Brussel in de heuvelstreken',
 'Zanden van Egem',
 'Zanden van Egem in de heuvelstreken',
 'Zanden van Mons-en-Pévèle',
 'Zanden van Mons-en-Pévèle in de heuvelstreken'}

We can now build a search query to get the feature where hhz_naam equals Formatie van Brasschaat (+Merksplas):



In [21]:

    
naam_filter = PropertyIsEqualTo(propertyname='hhz_naam', literal='Formatie van Brasschaat (+Merksplas)')



In [22]:

    
hhz_poly = hhz.getfeature(
    typename='gw_varia:hhz',
    filter=etree.tostring(naam_filter.toXML()).decode("utf8")).read()

And subsequently use this geometry in a search for groundwater screens:



In [23]:

    
filter_search = GrondwaterFilterSearch()
df = filter_search.search(
    max_features = 100,
    location=GmlFilter(hhz_poly, Within),
    return_fields=('pkey_filter', 'x', 'y')
)

And plot the results on a map:



In [24]:

    
df['lon'], df['lat'] = zip(*map(convert_latlon, df['x'], df['y']))
# convert to list
loclist = df[['lat', 'lon']].values.tolist()



In [25]:

    
# 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=10)
marker_cluster = MarkerCluster().add_to(fmap)
for loc in range(0, len(loclist)):
    folium.Marker(loclist[loc], popup=df['pkey_filter'][loc]).add_to(marker_cluster)
fmap









    Out[25]:




Make this Notebook Trusted to load map: File -> Trust Notebook



In [ ]: