{ "cells": [ { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "# Examples of spatial querying" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "[![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/DOV-Vlaanderen/pydov/master?filepath=docs%2Fnotebooks%2Fspatial_querying.ipynb)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "In a first section, different examples will be shown on how to use existing webservices to setup spatial queries:\n", "\n", "- Filter data within a certain community using its geographic borders\n", "- Filter data within a geographic boundary of a feature in layer \"Bekken\"\n", "- Filter data within a geographic boundary of a feature in layer \"Hydrogeologische homogene zones\"\n", "\n", "In other situations, a spatial query is need to be defined by a file (e.g. an ESRI shapefile) containing information on the particular area of interest. In a second section of this notebook, the usage of file based queries is shown. \n", "\n", "__Note:__ To run the file based examples, the installation of the `fiona` package and/or the `geopandas` package is required. See the installation instructions." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "%matplotlib inline\n", "from io import BytesIO\n", "import matplotlib.pyplot as plt\n", "\n", "from owslib.etree import etree\n", "from owslib.fes import PropertyIsEqualTo\n", "from owslib.wfs import WebFeatureService\n", "\n", "from pydov.search.boring import BoringSearch\n", "from pydov.search.grondwaterfilter import GrondwaterFilterSearch\n", "\n", "from pydov.util.location import (\n", " GmlFilter,\n", " Within,\n", " GeometryFilter, \n", " GeopandasFilter\n", ")\n", "\n", "from pydov.util.owsutil import get_namespace\n", "\n", "# import the necessary modules (not included in the requirements of pydov!)\n", "import folium\n", "from folium.plugins import MarkerCluster\n", "from pyproj import Transformer\n", "import fiona\n", "from fiona.crs import from_epsg\n", "import geopandas\n", "\n", "# convert the coordinates to lat/lon for folium\n", "def convert_latlon(x1, y1):\n", " transformer = Transformer.from_crs(\"epsg:31370\", \"epsg:4326\", always_xy=True)\n", " x2,y2 = transformer.transform(x1, y1)\n", " return x2, y2\n", "\n", "# convert to bytes if necessary\n", "def to_bytes(data):\n", " if isinstance(data, bytes):\n", " return data\n", " elif isinstance(data, str):\n", " return data.encode('utf8')" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## Spatial querying using web services" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "### Filter data within a certain community using its geographic borders" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "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.\n", "\n", "We have to start by defining the third party WFS service that contains the featuretype (layer) that contains the geometry we're interested in:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "gemeentegrenzen = WebFeatureService(\n", " 'https://geo.api.vlaanderen.be/VRBG/wfs',\n", " version='1.1.0')" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "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.\n", "\n", "We can get the schema (i.e. all the available fields) of a layer to find a field to query on:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{'UIDN': 'decimal',\n", " 'OIDN': 'decimal',\n", " 'TERRID': 'decimal',\n", " 'NISCODE': 'string',\n", " 'NAAM': 'string',\n", " 'DATPUBLBS': 'date',\n", " 'NUMAC': 'string'}" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "gemeentegrenzen.get_schema('VRBG:Refgem')['properties']" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "We can now build a search query to get the feature (including geometry) of the community with the `NAAM` of `Lievegem`:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "naam_filter = PropertyIsEqualTo(propertyname='NAAM', literal='Lievegem')" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "By default, the feature is returned as a GML feature, which stands for `Geography Markup Language` and is a XML grammar used to express geographical features.\n", "We change the `outputFormat` to GML3.2, since pydov's GmlFilter expects GML version 3.2." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "gemeente_gml = gemeentegrenzen.getfeature(\n", " typename='VRBG:Refgem',\n", " filter=etree.tostring(naam_filter.toXML()).decode(\"utf8\"),\n", " outputFormat='text/xml; subtype=gml/3.2').read()\n" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Next, we'll use the GML feature inside a `GmlFilter` query to find all boreholes that are located `Within` the community borders:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[000/001] .\n" ] } ], "source": [ "bs = BoringSearch()\n", "df = bs.search(\n", " location=GmlFilter(gemeente_gml, Within),\n", " return_fields=('pkey_boring', 'x', 'y'))" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "The `df` DataFrame now contains the `pkey_boring` together with the x and y coordinates:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
pkey_boringxy
0https://www.dov.vlaanderen.be/data/boring/1894...99605.0198147.0
1https://www.dov.vlaanderen.be/data/boring/1894...99209.0198139.0
2https://www.dov.vlaanderen.be/data/boring/1894...99103.0198002.0
3https://www.dov.vlaanderen.be/data/boring/1894...98629.0198057.0
4https://www.dov.vlaanderen.be/data/boring/1894...98518.0197907.0
\n", "
" ], "text/plain": [ " pkey_boring x y\n", "0 https://www.dov.vlaanderen.be/data/boring/1894... 99605.0 198147.0\n", "1 https://www.dov.vlaanderen.be/data/boring/1894... 99209.0 198139.0\n", "2 https://www.dov.vlaanderen.be/data/boring/1894... 99103.0 198002.0\n", "3 https://www.dov.vlaanderen.be/data/boring/1894... 98629.0 198057.0\n", "4 https://www.dov.vlaanderen.be/data/boring/1894... 98518.0 197907.0" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "df.head()" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "We can show these result on the map as well:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "df['lon'], df['lat'] = zip(*map(convert_latlon, df['x'], df['y']))\n", "\n", "# convert to list\n", "loclist = df[['lat', 'lon']].values.tolist()" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
Make this Notebook Trusted to load map: File -> Trust Notebook
" ], "text/plain": [ "" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# initialize the Folium map on the centre of the selected locations, play with the zoom until ok\n", "fmap = folium.Map(location=[df['lat'].mean(), df['lon'].mean()], zoom_start=10)\n", "marker_cluster = MarkerCluster().add_to(fmap)\n", "for loc in range(0, len(loclist)):\n", " folium.Marker(loclist[loc], popup=df['pkey_boring'][loc]).add_to(marker_cluster)\n", "fmap" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "### Filter data within a geographic boundary of a feature in layer \"Bekken\"" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "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.\n", "\n", "We have to start by defining the third party WFS service that contains the featuretype (layer) that contains the geometry we're interested in:" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "bekkens = WebFeatureService(\n", " 'https://geo.api.vlaanderen.be/Watersystemen/wfs',\n", " version='1.1.0'\n", ")" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "We can list all the available featuretypes (layers) in the service if we are unsure which one to use:" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "['Watersystemen:WsBekken',\n", " 'Watersystemen:WsDeelbek',\n", " 'Watersystemen:WsStrmgeb']" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "list(bekkens.contents)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "We can also get the schema (i.e. all the available fields) of a layer to find a field to query on:" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{'OIDN': 'decimal',\n", " 'UIDN': 'decimal',\n", " 'BEKNR': 'short',\n", " 'BEKNAAM': 'string',\n", " 'STRMGEB': 'string'}" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "bekkens.get_schema('Watersystemen:WsBekken')['properties']" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "We can get all the distinct values of a specific field:" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{'Bekken van de Brugse Polders',\n", " 'Bekken van de Gentse Kanalen',\n", " 'Benedenscheldebekken',\n", " 'Bovenscheldebekken',\n", " 'Demerbekken',\n", " 'Denderbekken',\n", " 'Dijle- en Zennebekken',\n", " 'IJzerbekken',\n", " 'Leiebekken',\n", " 'Maasbekken',\n", " 'Netebekken'}" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "namespace = get_namespace(bekkens, 'Watersystemen:WsBekken')\n", "\n", "tree = etree.fromstring(to_bytes(bekkens.getfeature('Watersystemen:WsBekken', propertyname='BEKNAAM').read()))\n", "set((i.text for i in tree.findall('.//{%s}BEKNAAM' % namespace)))" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "We can now build a search query to get the feature (including geometry) of the waterbody with the `BEKNAAM` of `Demerbekken`:" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [], "source": [ "naam_filter = PropertyIsEqualTo(propertyname='BEKNAAM', literal='Demerbekken')" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [], "source": [ "bekken_poly = bekkens.getfeature(\n", " typename='Watersystemen:WsBekken',\n", " filter=etree.tostring(naam_filter.toXML()).decode(\"utf8\"),\n", " outputFormat='text/xml; subtype=gml/3.2').read()" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "And use this GML feature inside a `GmlFilter` query to find all groundwater screens that are located `Within` the waterbody:" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[000/001] .\n" ] } ], "source": [ "filter_search = GrondwaterFilterSearch()\n", "df = filter_search.search(\n", " max_features = 100, # Note - we limit the results here to 100 for the demo case\n", " location=GmlFilter(bekken_poly, Within),\n", " return_fields=('pkey_filter', 'x', 'y')\n", ")" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
pkey_filterxy
0https://www.dov.vlaanderen.be/data/filter/2006...185347.00186221.00
1https://www.dov.vlaanderen.be/data/filter/2003...195990.07173681.03
2https://www.dov.vlaanderen.be/data/filter/1900...186630.00184140.00
3https://www.dov.vlaanderen.be/data/filter/1900...219095.00181256.00
4https://www.dov.vlaanderen.be/data/filter/1991...179856.00184474.00
\n", "
" ], "text/plain": [ " pkey_filter x y\n", "0 https://www.dov.vlaanderen.be/data/filter/2006... 185347.00 186221.00\n", "1 https://www.dov.vlaanderen.be/data/filter/2003... 195990.07 173681.03\n", "2 https://www.dov.vlaanderen.be/data/filter/1900... 186630.00 184140.00\n", "3 https://www.dov.vlaanderen.be/data/filter/1900... 219095.00 181256.00\n", "4 https://www.dov.vlaanderen.be/data/filter/1991... 179856.00 184474.00" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "df.head()" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "We can show the result on the map:" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [], "source": [ "df['lon'], df['lat'] = zip(*map(convert_latlon, df['x'], df['y']))\n", "# convert to list\n", "loclist = df[['lat', 'lon']].values.tolist()" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
Make this Notebook Trusted to load map: File -> Trust Notebook
" ], "text/plain": [ "" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# initialize the Folium map on the centre of the selected locations, play with the zoom until ok\n", "fmap = folium.Map(location=[df['lat'].mean(), df['lon'].mean()], zoom_start=10)\n", "marker_cluster = MarkerCluster().add_to(fmap)\n", "for loc in range(0, len(loclist)):\n", " folium.Marker(loclist[loc], popup=df['pkey_filter'][loc]).add_to(marker_cluster)\n", "fmap" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "### Filter data within a geographic boundary of a feature in layer \"HHZ\"" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "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.\n", "\n", "We can query the HHZ layer from DOV using pydov:" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [], "source": [ "from pydov.search.generic import WfsSearch\n", "\n", "hhz = WfsSearch('gw_varia:hhz')" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Next we can get the field of the layer to find a field to query on:" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", " \n", "
\n", "

id - None

  • type: string
  • notnull: False
  • query: True
  • cost: 1
  • multivalue: False
\n", "
\n", " \n", " \n", "
\n", "

hhz_naam - None

  • type: string
  • notnull: False
  • query: True
  • cost: 1
  • multivalue: False
\n", "
\n", " \n", " \n", "
\n", "

opp_m2 - None

  • type: integer
  • notnull: False
  • query: True
  • cost: 1
  • multivalue: False
\n", "
\n", " \n", " \n", "
\n", "

hhz - None

  • type: string
  • notnull: False
  • query: True
  • cost: 1
  • multivalue: False
\n", "
\n", " \n", " \n", "
\n", "

geom - None

  • type: geometry
  • notnull: False
  • query: False
  • cost: 1
  • multivalue: False
\n", "
\n", "
" ], "text/plain": [ "{'id': {'name': 'id', 'definition': None, 'type': 'string', 'multivalue': False, 'notnull': False, 'query': True, 'cost': 1}, 'hhz_naam': {'name': 'hhz_naam', 'definition': None, 'type': 'string', 'multivalue': False, 'notnull': False, 'query': True, 'cost': 1}, 'opp_m2': {'name': 'opp_m2', 'definition': None, 'type': 'integer', 'multivalue': False, 'notnull': False, 'query': True, 'cost': 1}, 'hhz': {'name': 'hhz', 'definition': None, 'type': 'string', 'multivalue': False, 'notnull': False, 'query': True, 'cost': 1}, 'geom': {'name': 'geom', 'definition': None, 'type': 'geometry', 'multivalue': False, 'notnull': False, 'query': False, 'cost': 1}}" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "hhz.get_fields()" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "We can now build a search query to get the feature where `hhz_naam` equals `Formatie van Mol`:" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [], "source": [ "from owslib.fes2 import PropertyIsEqualTo as PropertyIsEqualTo2\n", "from pydov.search.fields import GeometryReturnField\n", "\n", "naam_filter = PropertyIsEqualTo2(\n", " propertyname='hhz_naam', literal='Formatie van Mol')" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[000/001] .\n" ] }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
idhhz_naamgeom
02Formatie van MolMULTIPOLYGON (((203474.0966 223036.6511, 20424...
\n", "
" ], "text/plain": [ " id hhz_naam geom\n", "0 2 Formatie van Mol MULTIPOLYGON (((203474.0966 223036.6511, 20424..." ] }, "execution_count": 23, "metadata": {}, "output_type": "execute_result" } ], "source": [ "hhz_poly = hhz.search(query=naam_filter, return_fields=('id', 'hhz_naam', GeometryReturnField('geom', epsg=31370)))\n", "hhz_poly" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "And turn this result into a GeoPandas GeoDataFrame:" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
idhhz_naamgeom
02Formatie van MolMULTIPOLYGON (((203474.097 223036.651, 204245....
\n", "
" ], "text/plain": [ " id hhz_naam geom\n", "0 2 Formatie van Mol MULTIPOLYGON (((203474.097 223036.651, 204245...." ] }, "execution_count": 24, "metadata": {}, "output_type": "execute_result" } ], "source": [ "hhz_poly = geopandas.GeoDataFrame(hhz_poly, geometry='geom', crs='EPSG:31370')\n", "hhz_poly" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "And subsequently use this geometry in a search for groundwater screens:" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[000/001] .\n" ] }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
pkey_filtergeom
0https://www.dov.vlaanderen.be/data/filter/1984...POINT Z (5.1758 51.2172 0)
1https://www.dov.vlaanderen.be/data/filter/1989...POINT Z (5.1796 51.2162 0)
2https://www.dov.vlaanderen.be/data/filter/2005...POINT Z (5.1231 51.3088 0)
3https://www.dov.vlaanderen.be/data/filter/1999...POINT Z (5.0182 51.2388 0)
4https://www.dov.vlaanderen.be/data/filter/2012...POINT Z (5.0096 51.2578 0)
.........
95https://www.dov.vlaanderen.be/data/filter/2003...POINT Z (5.2206 51.2186 0)
96https://www.dov.vlaanderen.be/data/filter/1978...POINT Z (5.0297 51.2747 0)
97https://www.dov.vlaanderen.be/data/filter/2022...POINT Z (5.0816 51.2411 0)
98https://www.dov.vlaanderen.be/data/filter/1984...POINT Z (5.1781 51.2167 0)
99https://www.dov.vlaanderen.be/data/filter/1900...POINT Z (5.13 51.2298 0)
\n", "

100 rows × 2 columns

\n", "
" ], "text/plain": [ " pkey_filter \\\n", "0 https://www.dov.vlaanderen.be/data/filter/1984... \n", "1 https://www.dov.vlaanderen.be/data/filter/1989... \n", "2 https://www.dov.vlaanderen.be/data/filter/2005... \n", "3 https://www.dov.vlaanderen.be/data/filter/1999... \n", "4 https://www.dov.vlaanderen.be/data/filter/2012... \n", ".. ... \n", "95 https://www.dov.vlaanderen.be/data/filter/2003... \n", "96 https://www.dov.vlaanderen.be/data/filter/1978... \n", "97 https://www.dov.vlaanderen.be/data/filter/2022... \n", "98 https://www.dov.vlaanderen.be/data/filter/1984... \n", "99 https://www.dov.vlaanderen.be/data/filter/1900... \n", "\n", " geom \n", "0 POINT Z (5.1758 51.2172 0) \n", "1 POINT Z (5.1796 51.2162 0) \n", "2 POINT Z (5.1231 51.3088 0) \n", "3 POINT Z (5.0182 51.2388 0) \n", "4 POINT Z (5.0096 51.2578 0) \n", ".. ... \n", "95 POINT Z (5.2206 51.2186 0) \n", "96 POINT Z (5.0297 51.2747 0) \n", "97 POINT Z (5.0816 51.2411 0) \n", "98 POINT Z (5.1781 51.2167 0) \n", "99 POINT Z (5.13 51.2298 0) \n", "\n", "[100 rows x 2 columns]" ] }, "execution_count": 25, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from pydov.search.fields import GeometryReturnField\n", "\n", "filter_search = GrondwaterFilterSearch()\n", "df = filter_search.search(\n", " max_features = 100, # Note - we limit the results here to 100 for the demo case\n", " location=GeopandasFilter(hhz_poly, Within),\n", " return_fields=('pkey_filter', GeometryReturnField('geom', 4326))\n", ")\n", "df" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "And plot the results on a map:" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [], "source": [ "loclist = [[point.xy[1][0], point.xy[0][0]] for point in df.geom]" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
Make this Notebook Trusted to load map: File -> Trust Notebook
" ], "text/plain": [ "" ] }, "execution_count": 27, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# initialize the Folium map on the centre of the selected locations, play with the zoom until ok\n", "fmap = folium.Map(zoom_start=12)\n", "marker_cluster = MarkerCluster().add_to(fmap)\n", "\n", "# add filters\n", "for loc in range(0, len(loclist)):\n", " folium.Marker(loclist[loc], popup=df['pkey_filter'][loc]).add_to(marker_cluster)\n", "\n", "# add HHZ polygon\n", "hhz_poly_series = geopandas.GeoSeries(hhz_poly[\"geom\"]).simplify(tolerance=0.001)\n", "hhz_poly_json = hhz_poly_series.to_json()\n", "hhz_poly_json = folium.GeoJson(data=hhz_poly_series, style_function=lambda x: {\"fillColor\": \"orange\"})\n", "hhz_poly_json.add_to(fmap)\n", "\n", "fmap.fit_bounds(fmap.get_bounds())\n", "fmap" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## File based spatial querying" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "### Using a vector file" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Intead of using a webservice, we can also use a __file__ as input to define a spatial query. The pydov package provides the `GeometryFilter` to convert a vector file (e.g. ESRI shape file) into a GML file which can be used in the same way as the previous examples. In this case, we want all the boringen within the polygons defined in the vector file:" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "We have an example shape file prepared here, but this could be any other shape file as well:" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [], "source": [ "shapefile = \"../../tests/data/util/location/polygon_multiple_31370.shp\"" ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[000/001] .\n" ] } ], "source": [ "df = bs.search(\n", " location=GeometryFilter(shapefile, Within),\n", " return_fields=('pkey_boring', 'x', 'y'))" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "[Fiona](https://pypi.org/project/Fiona/) is used to do the data transformation, which is why this is an additional dependency when using the `GeometryFilter` functionality." ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "The resulting data points can be represented on a map:" ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [], "source": [ "df['lon'], df['lat'] = zip(*map(convert_latlon, df['x'], df['y']))\n", "# convert to list\n", "loclist = df[['lat', 'lon']].values.tolist()" ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
Make this Notebook Trusted to load map: File -> Trust Notebook
" ], "text/plain": [ "" ] }, "execution_count": 31, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# initialize the Folium map on the centre of the selected locations, play with the zoom until ok\n", "fmap = folium.Map(location=[df['lat'].mean(), df['lon'].mean()], zoom_start=10)\n", "marker_cluster = MarkerCluster().add_to(fmap)\n", "for loc in range(0, len(loclist)):\n", " folium.Marker(loclist[loc], popup=df['pkey_boring'][loc]).add_to(marker_cluster)\n", "fmap" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "### Interact with GeoPandas" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "In the previous example, the vector file as available on disk was used to create the spatial query. However, when working on geospatial vector data in Python the usage of [GeoPandas](https://geopandas.readthedocs.io/en/latest/) is convenient as it combines the data table functionalities of the [Pandas](https://pandas.pydata.org/docs/) package with spatial functionalities. GeoPandas also relies on fiona to read vector files, so using GeoPandas DataFrames can be used to define spatial filters as well." ] }, { "cell_type": "code", "execution_count": 32, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
gml_idgeometryname
0polygon_multiple_31370.0POLYGON ((108636.15 194960.844, 109195.574 195...site 1
1polygon_multiple_31370.1POLYGON ((107485.786 196741.544, 108297.344 19...site 2
\n", "
" ], "text/plain": [ " gml_id \\\n", "0 polygon_multiple_31370.0 \n", "1 polygon_multiple_31370.1 \n", "\n", " geometry name \n", "0 POLYGON ((108636.15 194960.844, 109195.574 195... site 1 \n", "1 POLYGON ((107485.786 196741.544, 108297.344 19... site 2 " ] }, "execution_count": 32, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import geopandas as gpd\n", "\n", "shapefile = \"../../tests/data/util/location/polygon_multiple_31370.shp\"\n", "\n", "# User reads geometry-file from disk (shapefile, geojson,...)\n", "gdf = gpd.read_file(shapefile)\n", "gdf[\"name\"] = [\"site 1\", \"site 2\"] # add dummy names to each of the polygones in the vector file for the demo case\n", "gdf" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Using the GeoDataFrame directly to define a spatial query:" ] }, { "cell_type": "code", "execution_count": 33, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[000/001] .\n" ] } ], "source": [ "my_filter = GeopandasFilter(gdf, Within)\n", "\n", "bs = BoringSearch()\n", "df = bs.search(\n", " location=my_filter,\n", " return_fields=('pkey_boring', 'x', 'y'))" ] }, { "cell_type": "code", "execution_count": 34, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
pkey_boringxy
0https://www.dov.vlaanderen.be/data/boring/2018...108025.00196593.00
1https://www.dov.vlaanderen.be/data/boring/2019...107947.29196640.52
2https://www.dov.vlaanderen.be/data/boring/2020...107991.00196706.00
3https://www.dov.vlaanderen.be/data/boring/2022...107842.53196371.08
4https://www.dov.vlaanderen.be/data/boring/2023...108144.95196771.09
5https://www.dov.vlaanderen.be/data/boring/2022...108131.66196779.43
6https://www.dov.vlaanderen.be/data/boring/2025...107924.96196787.73
7https://www.dov.vlaanderen.be/data/boring/2025...107862.05196777.40
8https://www.dov.vlaanderen.be/data/boring/2025...107863.09196776.79
9https://www.dov.vlaanderen.be/data/boring/2025...107896.57196754.28
10https://www.dov.vlaanderen.be/data/boring/2025...107895.42196755.16
11https://www.dov.vlaanderen.be/data/boring/2025...107925.14196786.74
12https://www.dov.vlaanderen.be/data/boring/2026...108932.00194350.80
13https://www.dov.vlaanderen.be/data/boring/1893...108900.00194425.00
14https://www.dov.vlaanderen.be/data/boring/1894...107618.00196709.00
15https://www.dov.vlaanderen.be/data/boring/1894...107791.00196516.00
16https://www.dov.vlaanderen.be/data/boring/1895...109050.00194990.00
17https://www.dov.vlaanderen.be/data/boring/1895...108742.00194936.00
\n", "
" ], "text/plain": [ " pkey_boring x y\n", "0 https://www.dov.vlaanderen.be/data/boring/2018... 108025.00 196593.00\n", "1 https://www.dov.vlaanderen.be/data/boring/2019... 107947.29 196640.52\n", "2 https://www.dov.vlaanderen.be/data/boring/2020... 107991.00 196706.00\n", "3 https://www.dov.vlaanderen.be/data/boring/2022... 107842.53 196371.08\n", "4 https://www.dov.vlaanderen.be/data/boring/2023... 108144.95 196771.09\n", "5 https://www.dov.vlaanderen.be/data/boring/2022... 108131.66 196779.43\n", "6 https://www.dov.vlaanderen.be/data/boring/2025... 107924.96 196787.73\n", "7 https://www.dov.vlaanderen.be/data/boring/2025... 107862.05 196777.40\n", "8 https://www.dov.vlaanderen.be/data/boring/2025... 107863.09 196776.79\n", "9 https://www.dov.vlaanderen.be/data/boring/2025... 107896.57 196754.28\n", "10 https://www.dov.vlaanderen.be/data/boring/2025... 107895.42 196755.16\n", "11 https://www.dov.vlaanderen.be/data/boring/2025... 107925.14 196786.74\n", "12 https://www.dov.vlaanderen.be/data/boring/2026... 108932.00 194350.80\n", "13 https://www.dov.vlaanderen.be/data/boring/1893... 108900.00 194425.00\n", "14 https://www.dov.vlaanderen.be/data/boring/1894... 107618.00 196709.00\n", "15 https://www.dov.vlaanderen.be/data/boring/1894... 107791.00 196516.00\n", "16 https://www.dov.vlaanderen.be/data/boring/1895... 109050.00 194990.00\n", "17 https://www.dov.vlaanderen.be/data/boring/1895... 108742.00 194936.00" ] }, "execution_count": 34, "metadata": {}, "output_type": "execute_result" } ], "source": [ "df" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "A user can do all sort of spatial and non-spatial operations in Geopandas. For example, one is only interested in the Boreholes of _site 1_ and selects this polygon only:" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Using the `GeopandasFilter`, we can use this to create a spatial query:" ] }, { "cell_type": "code", "execution_count": 35, "metadata": {}, "outputs": [], "source": [ "gdf_subselection = gdf[gdf[\"name\"] == \"site 1\"]" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "As long as the result is a GeoDataFrame with a crs defined, it can be used as input to define a spatial filter:" ] }, { "cell_type": "code", "execution_count": 36, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[000/001] .\n" ] } ], "source": [ "my_filter = GeopandasFilter(gdf_subselection, Within)\n", "\n", "bs = BoringSearch()\n", "df = bs.search(\n", " location=my_filter,\n", " return_fields=('pkey_boring', 'x', 'y'))" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "The returned set of boreholes is now limited to Site 1:" ] }, { "cell_type": "code", "execution_count": 37, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
pkey_boringxy
0https://www.dov.vlaanderen.be/data/boring/2026...108932.0194350.8
1https://www.dov.vlaanderen.be/data/boring/1893...108900.0194425.0
2https://www.dov.vlaanderen.be/data/boring/1895...109050.0194990.0
3https://www.dov.vlaanderen.be/data/boring/1895...108742.0194936.0
\n", "
" ], "text/plain": [ " pkey_boring x y\n", "0 https://www.dov.vlaanderen.be/data/boring/2026... 108932.0 194350.8\n", "1 https://www.dov.vlaanderen.be/data/boring/1893... 108900.0 194425.0\n", "2 https://www.dov.vlaanderen.be/data/boring/1895... 109050.0 194990.0\n", "3 https://www.dov.vlaanderen.be/data/boring/1895... 108742.0 194936.0" ] }, "execution_count": 37, "metadata": {}, "output_type": "execute_result" } ], "source": [ "df" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "__Note:__ To select a single row (i.e. feature) of a GeoDataFrame, make sure the result is still a valid GeoDataFrame:" ] }, { "cell_type": "code", "execution_count": 38, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
gml_idgeometryname
0polygon_multiple_31370.0POLYGON ((108636.15 194960.844, 109195.574 195...site 1
\n", "
" ], "text/plain": [ " gml_id \\\n", "0 polygon_multiple_31370.0 \n", "\n", " geometry name \n", "0 POLYGON ((108636.15 194960.844, 109195.574 195... site 1 " ] }, "execution_count": 38, "metadata": {}, "output_type": "execute_result" } ], "source": [ "gdf.iloc[[0]]" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "instead of:" ] }, { "cell_type": "code", "execution_count": 39, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "gml_id polygon_multiple_31370.0\n", "geometry POLYGON ((108636.150020818 194960.844295764, 1...\n", "name site 1\n", "Name: 0, dtype: object" ] }, "execution_count": 39, "metadata": {}, "output_type": "execute_result" } ], "source": [ "gdf.iloc[0]" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "The latter won't be a valid input to define a spatial query." ] } ], "metadata": { "kernelspec": { "display_name": ".venv (3.13.5)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.13.5" } }, "nbformat": 4, "nbformat_minor": 4 }