Skip to content

Instantly share code, notes, and snippets.

@arbakker
Last active Mar 16, 2021
Embed
What would you like to do?
QGIS processing tools for the PDOK locatie server

PDOK Locatieserver QGIS Processing Tool

QGIS processing tools for the PDOK Locatieserver (LS). This gist contains two processing tools:

  • pdok-geocoder.py: for geocoding i.e. convert a text-based description of a location to a geometry; for instance add point geometries to a spreadsheet of addresses
  • pdok-reverse-geocoder.py: for reverse geocoding i.e. convert a geometry in a text-based description of a location; for instance convert point geometries to addresses

Tested with QGIS version 3.16-hannover, will most likely work with all QGIS version 3.X.

Installation steps:

  1. Download pdok-reverse-geocoder.py
  2. Open the QGIS Processing Toolbox and click Add Script to Toolbox... and browse to the downloaded pdok-reverse-geocoder.py file:

  1. Now the PDOK Reverse Geocoder tool should be available in your toolbox:

  1. Start using the PDOK Reverse Geocoder:

⚠️ Be aware that for each feature 1 or more queries are send to the PDOK Locatieserver.

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""pdok-geocoder.py: QGIS Processing tool for geocoding with the PDOK \
Locatieserver. Tested with QGIS version 3.16, but will probably work with any \
3.X version."""
# MIT License
# Copyright (c) 2021 Anton Bakker
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.
__author__ = "Anton Bakker"
__copyright__ = "Copyright 2021, Anton Bakker"
__license__ = "MIT"
__version__ = "1.0.0"
__maintainer__ = "Anton Bakker"
__email__ = "anton.bakker@kadaster.nl"
__date__ = "2021-02-05"
from qgis.PyQt.QtCore import QCoreApplication, QVariant
from qgis.core import (
QgsProject,
QgsProcessing,
QgsField,
QgsGeometry,
QgsPointXY,
QgsWkbTypes,
QgsFeature,
QgsUnitTypes,
QgsFeatureSink,
QgsCoordinateReferenceSystem,
QgsCoordinateTransform,
QgsProcessingAlgorithm,
QgsProcessingException,
QgsProcessingParameterDistance,
QgsProcessingParameterCrs,
QgsProcessingParameterFeatureSource,
QgsProcessingParameterEnum,
QgsProcessingParameterBoolean,
QgsProcessingParameterString,
QgsProcessingParameterNumber,
QgsProcessingParameterFeatureSink,
QgsProcessingParameterField,
)
from qgis import processing
import requests
import json
import re
class PDOKGeocoder(QgsProcessingAlgorithm):
"""
This processing tool queries the PDOK Locatieserver fe geocoder service for each point in the input
layer and adds the first result to the target attribute.
"""
GEOM_TYPE_MAP = {
"weg": QgsWkbTypes.MultiLineString,
"adres": QgsWkbTypes.Point,
"gemeente": QgsWkbTypes.MultiPolygon,
"postcode": QgsWkbTypes.Point,
"woonplaats": QgsWkbTypes.MultiPolygon,
}
USER_AGENT_HEADER = {"User-Agent": "qgis-pdok-processing-tools"}
def tr(self, string):
"""
Returns a translatable string with the self.tr() function.
"""
return QCoreApplication.translate("Processing", string)
def createInstance(self):
# Must return a new copy of your algorithm.
return PDOKGeocoder()
def name(self):
"""
Returns the unique algorithm name.
"""
return "pdok-geocoder"
def displayName(self):
"""
Returns the translated algorithm name.
"""
return self.tr("PDOK Geocoder")
def group(self):
"""
Returns the name of the group this algorithm belongs to.
"""
return self.tr("PDOK Tools")
def groupId(self):
"""
Returns the unique ID of the group this algorithm belongs
to.
"""
return "pdok-tools"
def shortHelpString(self):
"""
Returns a localised short help string for the algorithm.
"""
return self.tr(
'This is processing tool queries the PDOK Locatieserver (LS) geocoder service for each\
feature in the input layer. The geometry returned by the LS, \
based on the target attribute of the feature, will be added to the output layer.\
Layers without geometry such as csv and xslx based layers are also suported. \
Existing attributes will be overwritten in the output layer. To query based on\
postal code and house number, ensure your input data conforms to this format: \
"{postal-code} {house-nr}" (note the space separating the postal code en the house number).\
So for example "1071XX 1". See the LS documentation: \
https://github.com/PDOK/locatieserver/wiki/API-Locatieserver\n\
Parameters:\n\n\
- Input layer: for each feature the LS geocoder service will be queried\n\
- Attribute to geocode: attribute in input layer to query LS with\n\
- Geocode result type, default - "adres"\n\
- Output layer: resulting output layer\n\
- Target CRS: CRS of the resulting output layer\n\
- Retrieve actual geometry (instead of centroid), default - false: will return MultiLineString geometry for weg, and MultiPolygon for gemeente en woonplaats. Not applicable to adres and postcode.\n\
- Add x and Y attribute, default - false: add "x" and "y" attributes to the output layer containing the \
geometry centroid coordinates\n\
- Add "weergavenaam" (display name) attribute, default - false: add "weergavenaam" attribute to the output \
layer, displayname is a field returned by LS.\n\
- Score treshold, optional: objects returned by the LS geocoder each have a score, \
to indicate how well they match with the query. Results with a score lower than the treshold \
are excluded\n\
'
)
def initAlgorithm(self, config=None):
"""
Here we define the inputs and outputs of the algorithm.
"""
self.predicates = [
("adres", self.tr("adres")),
("gemeente", self.tr("gemeente")),
("postcode", self.tr("postcode")),
("weg", self.tr("weg")),
("woonplaats", self.tr("woonplaats")),
]
self.TARGET_CRS = "TARGET_CRS"
self.INPUT = "INPUT" # recommended name for the main input parameter
self.ADD_XY_FIELD = "ADD_XY_FIELD"
self.SRC_FIELD = "SRC_FIELD"
self.RESULT_TYPE = "RESULT_TYPE"
self.SCORE_TRESHOLD = "SCORE_TRESHOLD"
self.OUTPUT = "OUTPUT" # recommended name for the main output parameter
self.ADD_DISPLAY_NAME = "ADD_DISPLAY_NAME"
self.GET_ACTUAL_GEOM = "GET_ACTUAL_GEOM"
self.addParameter(
QgsProcessingParameterFeatureSource(
self.INPUT,
self.tr("Input layer"),
types=[QgsProcessing.TypeFile],
)
)
self.addParameter(
QgsProcessingParameterField(
self.SRC_FIELD,
self.tr("Attribute to geocode"),
None,
"INPUT",
)
)
self.addParameter(
QgsProcessingParameterEnum(
self.RESULT_TYPE,
self.tr("Result type to geocode"),
options=[p[1] for p in self.predicates],
defaultValue=0,
optional=False,
)
)
self.addParameter(
QgsProcessingParameterFeatureSink(self.OUTPUT, self.tr("Output layer"))
)
self.addParameter(
QgsProcessingParameterCrs(
self.TARGET_CRS, self.tr("Target CRS"), "EPSG:4326"
)
)
self.addParameter(
QgsProcessingParameterBoolean(
self.GET_ACTUAL_GEOM,
self.tr("Retrieve actual geometry (instead of centroid)"),
False,
)
)
self.addParameter(
QgsProcessingParameterBoolean(
self.ADD_XY_FIELD, self.tr("Add x and Y attribute"), False
)
)
self.addParameter(
QgsProcessingParameterBoolean(
self.ADD_DISPLAY_NAME,
self.tr('Add "weergavenaam" (display name) attribute'),
False,
)
)
self.addParameter(
QgsProcessingParameterNumber(
self.SCORE_TRESHOLD,
self.tr("Score treshold"),
type=QgsProcessingParameterNumber.Double,
defaultValue=None,
optional=True,
minValue=0,
)
)
def get_geom(self, get_actual_geom, result_type, data, feedback):
"""
Returns a geometry depending on get_actual_geom boolean.
If false: return geom based on "centroide_ll" from the data
If true: retrieve the actual object from the lookup service and
return the geom based on "geometrie_ll" from the lookup response
"""
if not get_actual_geom or result_type in ["adres", "postcode"]:
wkt_point = data["response"]["docs"][0]["centroide_ll"]
return QgsGeometry.fromWkt(wkt_point)
else:
ls_id = wkt_point = data["response"]["docs"][0]["id"]
url = f"https://geodata.nationaalgeoregister.nl/locatieserver/v3/lookup?id={ls_id}&fl=id,geometrie_ll"
feedback.pushInfo(f"INFO: HTTP GET {url}")
response = requests.get(url, headers=PDOKGeocoder.USER_AGENT_HEADER)
if response.status_code != 200:
raise QgsProcessingException(
f"Unexpected response from HTTP GET {url}, response code: {response.status_code}"
)
data = response.json()
if len(data["response"]["docs"]) != 1:
raise QgsProcessingException(
f"Unexpected response body from HTTP GET {url}"
)
wkt_geom = data["response"]["docs"][0]["geometrie_ll"]
return QgsGeometry.fromWkt(wkt_geom)
def processAlgorithm(self, parameters, context, feedback):
try:
# read out parameters
input_layer = self.parameterAsVectorLayer(parameters, self.INPUT, context)
out_crs = parameters[self.TARGET_CRS]
result_type = [
self.predicates[i][0]
for i in self.parameterAsEnums(parameters, self.RESULT_TYPE, context)
][0]
score_treshold = parameters[self.SCORE_TRESHOLD]
add_xy_field = parameters[self.ADD_XY_FIELD]
add_display_name = parameters[self.ADD_DISPLAY_NAME]
src_field = parameters[self.SRC_FIELD]
get_actual_geom = parameters[self.GET_ACTUAL_GEOM]
# start processing
transform = None
fields = input_layer.fields()
field_names = [field.name() for field in fields]
if add_xy_field:
fields.append(QgsField("x", QVariant.Double))
fields.append(QgsField("y", QVariant.Double))
display_name_att_name = "weergavenaam"
if add_display_name:
fields.append(QgsField(display_name_att_name, QVariant.String))
(sink, dest_id) = self.parameterAsSink(
parameters,
self.OUTPUT,
context,
fields,
PDOKGeocoder.GEOM_TYPE_MAP[result_type],
out_crs,
)
if feedback.isCanceled():
return {}
for feature in input_layer.getFeatures():
query = feature.attribute(src_field)
url = f"https://geodata.nationaalgeoregister.nl/locatieserver/v3/free?q={query} and type:{result_type}"
match = re.search("([0-9]{4}[A-Za-z]{2})\s(.*)", query)
if match and len(match.groups()) == 2:
postal_code = match.group(1)
house_nr = match.group(2)
url = f"http://geodata.nationaalgeoregister.nl/locatieserver/free?fq=postcode:{postal_code}&fq=huisnummer~{house_nr}*&q=type:{result_type}"
feedback.pushInfo(f"INFO: HTTP GET {url}")
response = requests.get(url, headers=PDOKGeocoder.USER_AGENT_HEADER)
sc = response.status_code
if response.status_code != 200:
raise QgsProcessingException(
f"Unexpected response from HTTP GET {url}, response code: {response.status_code}"
)
data = response.json()
geom = None
display_name = ""
if len(data["response"]["docs"]) > 0:
if (
score_treshold != None
and data["response"]["docs"][0]["score"] <= score_treshold
):
geom = None
else:
geom = self.get_geom(
get_actual_geom, result_type, data, feedback
)
display_name = data["response"]["docs"][0]["weergavenaam"]
if geom:
attrs = feature.attributes()
new_ft = QgsFeature(fields)
for i in range(len(attrs)):
attr = attrs[i]
field_name = field_names[i]
new_ft.setAttribute(field_name, attr)
in_crs = QgsCoordinateReferenceSystem.fromEpsgId(4326)
if out_crs.authid() != "EPSG:4326":
transform = QgsCoordinateTransform(
in_crs, out_crs, QgsProject.instance()
)
geom.transform(transform)
if add_xy_field:
point_geom = QgsGeometry.asPoint(geom.centroid())
pxy = QgsPointXY(point_geom)
x = pxy.x()
y = pxy.y()
new_ft.setAttribute("x", x)
new_ft.setAttribute("y", y)
if add_display_name:
new_ft.setAttribute(display_name_att_name, display_name)
new_ft.setGeometry(geom)
sink.addFeature(new_ft, QgsFeatureSink.FastInsert)
if feedback.isCanceled():
return {}
results = {}
results[self.OUTPUT] = dest_id
return results
except Exception as e:
raise QgsProcessingException(
f"Unexpected error occured while running PDOKGeocoder: {str(e)}"
)
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""pdok-reverse-geocoder.py: QGIS Processing tool for reverse geocoding with the PDOK \
Locatieserver. Tested with QGIS version 3.16, but will probably work with any \
3.X version.
"""
# MIT License
# Copyright (c) 2021 Anton Bakker
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.
__author__ = "Anton Bakker"
__copyright__ = "Copyright 2021, Anton Bakker"
__license__ = "MIT"
__version__ = "1.0.0"
__maintainer__ = "Anton Bakker"
__email__ = "anton.bakker@kadaster.nl"
__date__ = "2021-02-05"
from qgis.PyQt.QtCore import QCoreApplication, QVariant
from qgis.core import (
QgsProject,
QgsProcessing,
QgsField,
QgsGeometry,
QgsPointXY,
QgsWkbTypes,
QgsFeature,
QgsUnitTypes,
QgsFeatureSink,
QgsCoordinateReferenceSystem,
QgsCoordinateTransform,
QgsProcessingAlgorithm,
QgsProcessingException,
QgsProcessingParameterDistance,
QgsProcessingParameterFeatureSource,
QgsProcessingParameterEnum,
QgsProcessingParameterString,
QgsProcessingParameterFeatureSink,
)
from qgis import processing
import requests
import json
class PDOKReverseGeocoder(QgsProcessingAlgorithm):
"""
This processing tool queries the PDOK Locatieserver reverse geocoder service for each point in the input
layer and adds the first result to the target attribute.
"""
USER_AGENT_HEADER = {"User-Agent": "qgis-pdok-processing-tools"}
def tr(self, string):
"""
Returns a translatable string with the self.tr() function.
"""
return QCoreApplication.translate("Processing", string)
def createInstance(self):
# Must return a new copy of your tool.
return PDOKReverseGeocoder()
def name(self):
"""
Returns the unique tool name.
"""
return "pdok-reverse-geocoder"
def displayName(self):
"""
Returns the translated tool name.
"""
return self.tr("PDOK Reverse Geocoder")
def group(self):
"""
Returns the name of the group this tool belongs to.
"""
return self.tr("PDOK Tools")
def groupId(self):
"""
Returns the unique ID of the group this tool belongs
to.
"""
return "pdok-tools"
def shortHelpString(self):
"""
Returns a localised short help string for the tool.
"""
return self.tr(
'This processing tool queries the PDOK Locatieserver reverse geocoder service for each\
point in the input layer and adds the first result to the target attribute.\n\
See PDOK Locatieserver documentation: https://github.com/PDOK/locatieserver/wiki/API-Reverse-Geocoder\n\
Parameters:\n\n\
- Input point layer (any projection): for each point the PDOK locatieserver reverse geocoder service will be queried\n\
- Result type to query: defaults to "adres"\n\
- Distance treshold, optional: objects returned by the PDOK locatieserver reverse geocoder \
with a distance greater than the threshold will be excluded\n\
- Attribute name, optional: defaults to result type, target attribute name results will be written to'
)
def initAlgorithm(self, config=None):
"""
Here we define the inputs and outputs of the tool.
"""
self.predicates = [
("adres", self.tr("adres")),
("appartementsrecht", self.tr("appartementsrecht")),
("buurt", self.tr("buurt")),
("gemeente", self.tr("gemeente")),
("hectometerpaal", self.tr("hectometerpaal")),
("perceel", self.tr("perceel")),
("postcode", self.tr("postcode")),
("provincie", self.tr("provincie")),
("waterschap", self.tr("waterschap")),
("weg", self.tr("weg")),
("wijk", self.tr("wijk")),
("woonplaats", self.tr("woonplaats")),
]
self.INPUT = "INPUT" # recommended name for the main input parameter
self.ATTRIBUTE_NAME = "ATTRIBUTE_NAME"
self.RESULT_TYPE = "RESULT_TYPE"
self.DISTANCE_TRESHOLD = "DISTANCE_TRESHOLD"
self.OUTPUT = "OUTPUT" # recommended name for the main output parameter
self.addParameter(
QgsProcessingParameterFeatureSource(
self.INPUT,
self.tr("Input point layer"),
types=[QgsProcessing.TypeVectorPoint],
)
)
self.addParameter(
QgsProcessingParameterEnum(
self.RESULT_TYPE,
self.tr("Result type to query"),
options=[p[1] for p in self.predicates],
defaultValue=0,
optional=True,
)
)
dist_param = QgsProcessingParameterDistance(
self.DISTANCE_TRESHOLD,
self.tr("Distance treshold"),
defaultValue=None,
optional=True,
minValue=0,
)
dist_param.setDefaultUnit(QgsUnitTypes.DistanceMeters)
self.addParameter(dist_param)
self.addParameter(
QgsProcessingParameterString(
self.ATTRIBUTE_NAME,
self.tr("Attribute name"),
defaultValue=None,
optional=True,
)
)
self.addParameter(
QgsProcessingParameterFeatureSink(
self.OUTPUT, self.tr("Output point layer")
)
)
def processAlgorithm(self, parameters, context, feedback):
try:
# read out algorithm parameters
input_points = self.parameterAsVectorLayer(parameters, self.INPUT, context)
distance_treshold = parameters[self.DISTANCE_TRESHOLD]
target_att_name = parameters[self.ATTRIBUTE_NAME]
result_type = [
self.predicates[i][0]
for i in self.parameterAsEnums(parameters, self.RESULT_TYPE, context)
][0]
# Initialize output layer
fields = input_points.fields()
field_names = [field.name() for field in fields]
if not target_att_name:
target_att_name = result_type
if target_att_name in fields:
raise QgsProcessingException(
f"Target attribute name {field_name} already exists in input layer. \
Supply different target attribute name."
)
fields.append(QgsField(target_att_name, QVariant.String))
(sink, dest_id) = self.parameterAsSink(
parameters,
self.OUTPUT,
context,
fields,
QgsWkbTypes.Point,
input_points.sourceCrs(),
)
# Setup transformation if required
in_crs = input_points.crs()
out_crs = QgsCoordinateReferenceSystem.fromEpsgId(4326)
transform = None
if in_crs.authid() != "EPSG:4326":
transform = QgsCoordinateTransform(
in_crs, out_crs, QgsProject.instance()
)
if feedback.isCanceled():
return {}
# start processing features
for point in input_points.getFeatures():
geom = point.geometry()
if transform:
geom.transform(transform)
point_geom = QgsGeometry.asPoint(geom)
pxy = QgsPointXY(point_geom)
lon = pxy.x()
lat = pxy.y()
url = f"https://geodata.nationaalgeoregister.nl/locatieserver/v4/revgeo/?lon={lon}&type={result_type}&lat={lat}"
feedback.pushInfo(f"INFO: HTTP GET {url}")
response = requests.get(
url, headers=PDOKReverseGeocoder.USER_AGENT_HEADER
)
sc = response.status_code
if response.status_code != 200:
raise QgsProcessingException(
f"Unexpected response from HTTP GET {url}, response code: {response.status_code}"
)
data = response.json()
result = ""
if len(data["response"]["docs"]) > 0:
if (
distance_treshold != None
and data["response"]["docs"][0]["afstand"] > distance_treshold
):
pass
else:
result = data["response"]["docs"][0]["weergavenaam"]
attrs = point.attributes()
new_ft = QgsFeature(fields)
for i in range(len(attrs)):
attr = attrs[i]
field_name = field_names[i]
new_ft.setAttribute(field_name, attr)
new_ft.setAttribute(target_att_name, result)
new_ft.setGeometry(point.geometry())
sink.addFeature(new_ft, QgsFeatureSink.FastInsert)
if feedback.isCanceled():
return {}
results = {}
results[self.OUTPUT] = dest_id
return results
except Exception as e:
raise QgsProcessingException(
f"Unexpected error occured while running PDOKReverseGeocoder: {str(e)}"
)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment