Coverage for brodata / webservices.py: 65%
159 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-03-20 14:37 +0000
« prev ^ index » next coverage.py v7.13.5, created at 2026-03-20 14:37 +0000
1import logging
2import json
4import geopandas as gpd
5import numpy as np
6import requests
7from requests.exceptions import HTTPError
8from shapely.geometry import MultiPolygon, Point, Polygon
9from .util import tqdm
11logger = logging.getLogger(__name__)
14def get_gdf(kind, extent=None, config=None, index="DINO_NR", **kwargs):
15 """Retrieve a GeoDataFrame for a configured BRO/GDN service.
17 Parameters
18 ----------
19 kind : str
20 Name of the service as defined in the configuration.
21 extent : list, tuple, np.array, shapely.geometry.Polygon or shapely.geometry.MultiPolygon, optional
22 Spatial filter. Use [xmin, xmax, ymin, ymax] for an envelope query or
23 a polygon geometry for an intersects query.
24 config : dict, optional
25 Service configuration dictionary. If omitted, the default configuration
26 from get_configuration() is used.
27 index : str, optional
28 Column name to use as index when present in the returned data.
29 **kwargs
30 Extra keyword arguments passed to requests.get via arcrest.
32 Returns
33 -------
34 geopandas.GeoDataFrame
35 GeoDataFrame with optional post-filtering and index assignment.
36 """
37 if config is None:
38 config = get_configuration()
39 if kind not in config:
40 raise (ValueError(f"Unknown kind: {kind}"))
41 url = config[kind]["mapserver"]
42 layer = 0
43 if "layer" in config[kind]:
44 layer = config[kind]["layer"]
45 gdf = arcrest(url, layer, extent=extent, **kwargs)
46 if "greater_than_0" in config[kind]:
47 greater_than_0 = config[kind]["greater_than_0"]
48 if greater_than_0 in gdf.columns:
49 gdf = gdf[gdf[greater_than_0] > 0]
50 if not gdf.empty and index in gdf.columns:
51 gdf = gdf.set_index(index)
52 return gdf
55def arcrest(
56 url,
57 layer,
58 extent=None,
59 sr=28992,
60 f="geojson",
61 max_record_count=None,
62 timeout=120,
63 **kwargs,
64):
65 """Download data from an arcgis rest FeatureServer.
67 Parameters
68 ----------
69 url : str
70 arcrest url.
71 layer : str
72 layer
73 extent : list, tuple, np.array, shapely.geometry.Polygon or shapely.geometry.MultiPolygon
74 spatial filter as envelope [xmin, xmax, ymin, ymax] or polygon geometry.
75 sr : int, optional
76 coördinate reference system. The default is 28992 (RD).
77 f : str, optional
78 output format. Default is geojson
79 max_record_count : int, optional
80 maximum number of records for request.
81 timeout : int, optional
82 timeout time of request. Default is 120.
84 """
85 params = {
86 "f": f,
87 "outFields": "*",
88 "outSR": sr,
89 "where": "1=1",
90 }
91 if extent is not None:
92 params["spatialRel"] = "esriSpatialRelIntersects"
93 if isinstance(extent, (Polygon, MultiPolygon)):
94 params["geometry"] = json.dumps(_as_arcgis_polygon(extent, sr=sr))
95 params["geometryType"] = "esriGeometryPolygon"
96 else:
97 xmin, xmax, ymin, ymax = extent
98 params["geometry"] = f"{xmin},{ymin},{xmax},{ymax}"
99 params["geometryType"] = "esriGeometryEnvelope"
100 params["inSR"] = sr
101 props = _get_data(url, {"f": "json"}, timeout=timeout, **kwargs)
102 if max_record_count is None:
103 max_record_count = props["maxRecordCount"]
104 else:
105 max_record_count = min(max_record_count, props["maxRecordCount"])
107 params["returnIdsOnly"] = True
108 url_query = f"{url}/{layer}/query"
109 props = _get_data(url_query, params, timeout=timeout, **kwargs)
110 params.pop("returnIdsOnly")
111 if "objectIds" in props:
112 object_ids = props["objectIds"]
113 object_id_field_name = props["objectIdFieldName"]
114 else:
115 object_ids = props["properties"]["objectIds"]
116 object_id_field_name = props["properties"]["objectIdFieldName"]
117 if object_ids is not None and len(object_ids) > max_record_count:
118 # download in batches
119 object_ids.sort()
120 n_d = int(np.ceil((len(object_ids) / max_record_count)))
121 features = []
122 for i_d in tqdm(range(n_d)):
123 i_min = i_d * max_record_count
124 i_max = min(i_min + max_record_count - 1, len(object_ids) - 1)
125 where = "{}>={} and {}<={}".format(
126 object_id_field_name,
127 object_ids[i_min],
128 object_id_field_name,
129 object_ids[i_max],
130 )
131 params["where"] = where
132 data = _get_data(url_query, params, timeout=timeout, **kwargs)
133 features.extend(data["features"])
134 else:
135 # download all data in one go
136 data = _get_data(url_query, params, timeout=timeout, **kwargs)
137 features = data["features"]
138 if f == "json" or f == "pjson":
139 # Interpret the geometry field
140 data = []
141 for feature in features:
142 if "rings" in feature["geometry"]:
143 if len(feature["geometry"]) > 1:
144 raise (NotImplementedError("Multiple rings not supported yet"))
145 if len(feature["geometry"]["rings"]) == 1:
146 geometry = Polygon(feature["geometry"]["rings"][0])
147 else:
148 pols = [Polygon(xy) for xy in feature["geometry"]["rings"]]
149 keep = [0]
150 for i in range(1, len(pols)):
151 if pols[i].within(pols[keep[-1]]):
152 pols[keep[-1]] = pols[keep[-1]].difference(pols[i])
153 else:
154 keep.append(i)
155 if len(keep) == 1:
156 geometry = pols[keep[0]]
157 else:
158 geometry = MultiPolygon([pols[i] for i in keep])
159 elif (
160 len(feature["geometry"]) == 2
161 and "x" in feature["geometry"]
162 and "y" in feature["geometry"]
163 ):
164 geometry = Point(feature["geometry"]["x"], feature["geometry"]["y"])
165 else:
166 raise (Exception("Not supported yet"))
167 feature["attributes"]["geometry"] = geometry
168 data.append(feature["attributes"])
169 if len(data) == 0:
170 # Assigning CRS to a GeoDataFrame without a geometry column is not supported
171 gdf = gpd.GeoDataFrame()
172 else:
173 gdf = gpd.GeoDataFrame(data, crs=sr)
174 else:
175 # for geojson-data we can transform to GeoDataFrame right away
176 if len(features) == 0:
177 # Assigning CRS to a GeoDataFrame without a geometry column is not supported
178 gdf = gpd.GeoDataFrame()
179 else:
180 gdf = gpd.GeoDataFrame.from_features(features, crs=sr)
182 return gdf
185def _as_arcgis_polygon(polygon, sr=28992):
186 """Convert polygon input to ArcGIS REST polygon geometry."""
187 if isinstance(polygon, dict):
188 geometry = polygon.copy()
189 if "rings" not in geometry:
190 raise ValueError("polygon dict must contain 'rings'")
191 elif isinstance(polygon, Polygon):
192 rings = [list(polygon.exterior.coords)]
193 for interior in polygon.interiors:
194 rings.append(list(interior.coords))
195 geometry = {"rings": rings}
196 elif isinstance(polygon, MultiPolygon):
197 rings = []
198 for pol in polygon.geoms:
199 rings.append(list(pol.exterior.coords))
200 for interior in pol.interiors:
201 rings.append(list(interior.coords))
202 geometry = {"rings": rings}
203 else:
204 raise TypeError(
205 "polygon must be a shapely Polygon, MultiPolygon, or dict with 'rings'"
206 )
208 if "spatialReference" not in geometry:
209 geometry["spatialReference"] = {"wkid": sr}
210 return geometry
213def _get_data(url, params, timeout=5, **kwargs):
214 """get data using a request
216 Parameters
217 ----------
218 url : str
219 url
220 params : dict
221 request parameters
222 timeout : int, optional
223 timeout time of request. Default is 120.
225 Returns
226 -------
227 data
229 """
230 r = requests.get(url, params=params, timeout=timeout, **kwargs)
231 if not r.ok:
232 raise (HTTPError(f"Request not successful: {r.url}"))
233 data = r.json()
234 if "error" in data:
235 code = data["error"]["code"]
236 message = data["error"]["message"]
237 raise (Exception(f"Error code {code}: {message}"))
238 return data
241services = {
242 "lks_abo_rd": "Archeologisch booronderzoek",
243 # "lks_aeo_rd": "",
244 # "lks_aep_rd": "",
245 "lks_bhr_g_rd": "Geologisch booronderzoek (BRO)",
246 "lks_bhr_gt_rd": "Geotechnisch booronderzoek (BRO)",
247 "lks_bhr_rd": "Bodemkundig booronderzoek (BRO)",
248 "lks_cpt_rd": "Geotechnisch sondeeronderzoek (BRO)",
249 # "lks_dgm_rd": "",
250 # "lks_dgmdiep_rd": "",
251 "lks_gbo_rd": "Geologisch booronderzoek",
252 # "lks_geok_line_rd": "",
253 # "lks_geok_rd": "",
254 # "lks_geok_urania_rd": "",
255 # "lks_gmd5_rd": "",
256 # "lks_gmm_rd": "",
257 "lks_gmw_rd": "Grondwatermonitoringput",
258 "lks_gso_rd": "Geotechnisch sondeeronderzoek (GDN)",
259 # "lks_gtp_rd": "",
260 "lks_guf_rd": "Grondwatergebruiksystemen",
261 "lks_gwo_rd": "Put met onderzoekgegevens",
262 # "lks_mdl_rd": "",
263 # "lks_nzs_rd": "",
264 "lks_owo_rd": "Oppervlaktewateronderzoek",
265 # "lks_rgs_rd": "",
266 "lks_sfr_rd": "Bodemkundig wandonderzoek (BRO)",
267 # "lks_sgm2_rd": "",
268 "lks_vso_rd": "Verticaal elektrisch sondeeronderzoek",
269 "lks_wbo_rd": "Geologisch waterbodemonderzoek",
270 # "lks_wdm_brh_rd": "",
271 # "lks_wdm_ghg_rd": "",
272 # "lks_wdm_glg_rd": "",
273 # "lks_wdm_gvg_rd": "",
274 #'lks_wdm_gxg_rd": "",
275 #'lks_wdm_ref_rd": ""
276}
279def get_configuration(mapserver_url=None):
280 config = {}
281 if mapserver_url is None:
282 mapserver_url = "https://www.broloket.nl/standalone/rest/services"
284 # BRO
285 bro_mapserver_url = f"{mapserver_url}/uitgifteloket_bro"
286 config["Geologisch booronderzoek (BRO)"] = {
287 "mapserver": f"{bro_mapserver_url}/lks_bhr_g_rd_v1/MapServer",
288 "abbr": "bhrg",
289 "rest_url": "https://publiek.broservices.nl/sr/bhrg/v3",
290 }
291 config["Geotechnisch booronderzoek (BRO)"] = {
292 "mapserver": f"{bro_mapserver_url}/lks_bhr_gt_rd_v1/MapServer",
293 "abbr": "bhrgt",
294 "object": "BHR-GT",
295 "rest_url": "https://publiek.broservices.nl/sr/bhrgt/v2",
296 }
297 config["Bodemkundig booronderzoek (BRO)"] = {
298 "mapserver": f"{bro_mapserver_url}/lks_bhr_rd_v1/MapServer",
299 "abbr": "bhrp",
300 "object": "BHR_O",
301 "rest_url": "https://publiek.broservices.nl/sr/bhrp/v2",
302 }
303 config["Geotechnisch sondeeronderzoek (BRO)"] = {
304 "mapserver": f"{bro_mapserver_url}/lks_cpt_rd_v1/MapServer",
305 "abbr": "cpt",
306 "rest_url": "https://publiek.broservices.nl/sr/cpt/v1",
307 }
308 config["Bodemkundig wandonderzoek (BRO)"] = {
309 "mapserver": f"{bro_mapserver_url}/lks_sfr_rd_v1/MapServer",
310 "name_en": "pedological SoilFaceResearch",
311 "abbr": "sfr",
312 }
313 config["Grondwatermonitoringput"] = {
314 "mapserver": f"{bro_mapserver_url}/lks_gmw_rd_v1/MapServer",
315 "rest_url": "https://publiek.broservices.nl/gm/gmw/v1",
316 }
317 config["Grondwaterstandonderzoek"] = {
318 "reeks": "https://publiek.broservices.nl/gm/gld/v1/seriesAsCsv",
319 "rest_url": "https://publiek.broservices.nl/gm/gld/v1",
320 "abbr": "gld",
321 }
322 config["Grondwatersamenstellingsonderzoek"] = {
323 "abbr": "gar",
324 "rest_url": "https://publiek.broservices.nl/gm/gar/v1",
325 "object": "GAR_O",
326 }
327 config["Grondwatergebruiksystemen"] = {
328 "mapserver": f"{bro_mapserver_url}/lks_guf_rd_v1/MapServer",
329 "name_nl": "Grondwatergebruiksystemen",
330 "name_en": "Groundwater utilisation facility",
331 "class": "GroundwaterUtilisationFacility",
332 "abbr": "guf",
333 "object": "GUF_PO",
334 "rest_url": "https://publiek.broservices.nl/gu/guf/v1",
335 }
337 # DINO
338 gdn_mapserver_url = f"{mapserver_url}/uitgifteloket_gdn"
339 config["Geologisch booronderzoek"] = {
340 "mapserver": f"{gdn_mapserver_url}/lks_gbo_rd_v1/MapServer",
341 "download": "https://www.dinoloket.nl/uitgifteloket/api/brh/sampledescription/csv",
342 }
343 config["Boormonsterprofiel"] = config["Geologisch booronderzoek"].copy()
344 config["Boormonsterprofiel"]["greater_than_0"] = "MP_CNT"
346 config["Boormonsterfoto"] = config["Geologisch booronderzoek"].copy()
347 config["Boormonsterfoto"]["greater_than_0"] = "MF_CNT"
349 config["Boorgatmeting"] = config["Geologisch booronderzoek"].copy()
350 config["Boorgatmeting"]["greater_than_0"] = "BM_CNT"
352 config["Chemische analyse"] = config["Geologisch booronderzoek"].copy()
353 config["Chemische analyse"]["greater_than_0"] = "CA_CNT"
355 config["Korrelgrootte analyse"] = config["Geologisch booronderzoek"].copy()
356 config["Korrelgrootte analyse"]["greater_than_0"] = "KA_CNT"
358 config["Geologisch waterbodemonderzoek"] = {
359 "mapserver": f"{gdn_mapserver_url}/lks_wbo_rd_v1/MapServer",
360 "download": config["Geologisch booronderzoek"]["download"],
361 }
362 config["Archeologisch booronderzoek"] = {
363 "mapserver": f"{gdn_mapserver_url}/lks_abo_rd_v1/MapServer",
364 "download": config["Geologisch booronderzoek"]["download"],
365 }
366 config["Geotechnisch sondeeronderzoek (GDN)"] = {
367 "mapserver": f"{gdn_mapserver_url}/lks_gso_rd_v1/MapServer",
368 }
369 config["Put met onderzoekgegevens"] = {
370 "mapserver": f"{gdn_mapserver_url}/lks_gwo_rd_tiled_v1/MapServer",
371 }
372 config["Grondwatersamenstelling"] = {
373 "mapserver": config["Put met onderzoekgegevens"]["mapserver"],
374 "table": 2, # Grondwatersamenstelling
375 "greater_than_0": "SA_CNT",
376 "download": "https://www.dinoloket.nl/uitgifteloket/api/wo/gwo/qua/report",
377 }
378 config["Grondwaterstand"] = {
379 "mapserver": config["Put met onderzoekgegevens"]["mapserver"],
380 "table": 3, # Grondwaterstand
381 "greater_than_0": "ST_CNT",
382 "download": "https://www.dinoloket.nl/uitgifteloket/api/wo/gwo/full",
383 "details": "https://www.dinoloket.nl/uitgifteloket/api/wo/gwo/details",
384 }
385 config["Oppervlaktewateronderzoek"] = {
386 "mapserver": f"{gdn_mapserver_url}/lks_owo_rd_v1/MapServer"
387 }
388 config["Verticaal elektrisch sondeeronderzoek"] = {
389 "mapserver": f"{gdn_mapserver_url}/lks_vso_rd_v1/MapServer",
390 "download": "https://www.dinoloket.nl/uitgifteloket/api/ves/csv",
391 }
393 return config