diff --git a/nlmod/read/ahn.py b/nlmod/read/ahn.py index 25f0370f..ec01c295 100644 --- a/nlmod/read/ahn.py +++ b/nlmod/read/ahn.py @@ -1,12 +1,16 @@ import datetime as dt import logging +import requests +from requests.exceptions import HTTPError import geopandas as gpd import matplotlib.pyplot as plt import numpy as np import pandas as pd +import shapely import rasterio import rioxarray +from rioxarray.merge import merge_arrays import xarray as xr from rasterio import merge from rasterio.io import MemoryFile @@ -15,14 +19,14 @@ from .. import cache from ..dims.grid import get_extent from ..dims.resample import structured_da_to_ds -from ..util import get_ds_empty +from ..util import get_ds_empty, extent_to_polygon from .webservices import arcrest, wcs logger = logging.getLogger(__name__) @cache.cache_netcdf(coords_2d=True) -def get_ahn(ds=None, identifier="AHN4_DTM_5m", method="average", extent=None): +def get_ahn(ds=None, identifier="AHN4_5M_M", method="average", extent=None, **kwargs): """Get a model dataset with ahn variable. Parameters @@ -31,11 +35,14 @@ def get_ahn(ds=None, identifier="AHN4_DTM_5m", method="average", extent=None): dataset with the model information. identifier : str, optional Possible values for the different AHN-versions are (casing is important): - AHN1: 'ahn1_5m' - AHN2: 'ahn2_05m_i', 'ahn2_05m_n', 'ahn2_05m_r' or 'ahn2_5m' - AHN3: 'AHN3_05m_DSM', 'AHN3_05m_DTM', 'AHN3_5m_DSM' or 'AHN3_5m_DTM' - AHN4: 'AHN4_DTM_05m', 'AHN4_DTM_5m', 'AHN4_DSM_05m' or 'AHN4_DSM_5m' - The default is 'AHN4_DTM_5m'. + AHN1: 'AHN1_5M' + AHN2: 'AHN2_05M_I', 'AHN2_05M_N', 'AHN2_05M_R' or 'AHN2_5M_M' + AHN3: 'AHN3_05M_M', 'AHN3_05M_R', 'AHN3_5M_M' or 'AHN3_5M_R' + AHN4: 'AHN4_05M_M', 'AHN4_05M_R', 'AHN4_5M_M' or 'AHN4_5M_R' + AHN5: 'AHN5_5M_M', 'AHN5_5M_R', 'AHN5_05M_M' or 'AHN5_05M_R' + The identifier determines the resolution (05M for 0.5 m and 5M for 5 m) and the + type of height data (M = DTM = surface level, R = DSM = also other features). + The default is 'AHN4_5M_M'. method : str, optional Method used to resample ahn to grid of ds. See documentation of nlmod.resample.structured_da_to_ds for possible values. The default is @@ -50,18 +57,7 @@ def get_ahn(ds=None, identifier="AHN4_DTM_5m", method="average", extent=None): """ if extent is None and ds is not None: extent = get_extent(ds) - version = int(identifier[3]) - if version == 1: - ahn_ds_raw = get_ahn1(extent, identifier=identifier) - elif version == 2: - ahn_ds_raw = get_ahn2(extent, identifier=identifier) - elif version == 3: - ahn_ds_raw = get_ahn3(extent, identifier=identifier) - elif version == 4: - ahn_ds_raw = get_ahn4(extent, identifier=identifier) - else: - raise (ValueError(f"Unknown ahn-version: {version}")) - + ahn_ds_raw = _get_ahn_ellipsis(extent, identifier=identifier, **kwargs) ahn_ds_raw = ahn_ds_raw.drop_vars("band") if ds is None: @@ -303,15 +299,80 @@ def get_ahn4_tiles(extent=None): return gdf +def _get_tiles_ellipsis( + extent=None, + crs=28992, + timeout=120.0, + base_url="https://api.ellipsis-drive.com/v3", + path_id="a9d410ad-a2f6-404c-948a-fdf6b43e77a6", + timestamp_id="87a21a71-c39f-4e92-a43b-207bc7dfe714", +): + url = f"{base_url}/path/{path_id}/vector/timestamp/{timestamp_id}/listFeatures" + + r = requests.get(url, timeout=timeout) + if not r.ok: + raise (HTTPError(f"Request not successful: {r.url}")) + gdf = gpd.GeoDataFrame.from_features(r.json()["result"]["features"], crs=4326) + gdf = gdf.to_crs(crs) + # remove small digits becuase of crs-transformation + gdf.geometry = gdf.geometry.apply(_round_coordinates, ndigits=0) + + gdf = gdf.set_index("AHN") + + if extent is not None: + gdf = gdf.loc[gdf.intersection(extent_to_polygon(extent)).area > 0] + return gdf + + +def _round_coordinates(geom, ndigits=2): + + def _round_coords(x, y, z=None): + x = round(x, ndigits) + y = round(y, ndigits) + + if z is not None: + z = round(x, ndigits) + return (x, y, z) + else: + return (x, y) + + return shapely.ops.transform(_round_coords, geom) + + @cache.cache_netcdf() -def get_ahn1(extent, identifier="ahn1_5m", as_data_array=True): +def get_ahn1(extent, identifier="AHN1_5M", as_data_array=None, **kwargs): """Download AHN1. Parameters ---------- extent : list, tuple or np.array extent - identifier : TYPE, optional + identifier : str, optional + Only allowed value is 'AHN1_5M' (or the equivalent ahn1_5m). The default is + "AHN1_5M". + + Returns + ------- + xr.DataArray + DataArray of the AHN + """ + _assert_as_data_array_is_none(as_data_array) + identifier = _rename_identifier(identifier) + da = _get_ahn_ellipsis(extent, identifier, **kwargs) + # original data is in cm. Convert the data to m, which is the unit of other ahns + da = da / 100 + return da + + +@cache.cache_netcdf() +def get_ahn1_legacy(extent, identifier="ahn1_5m", as_data_array=True): + """Download AHN1. + + Parameters + ---------- + extent : list, tuple or np.array + extent + identifier : str, optional Only allowed value is 'ahn1_5m'. The default is "ahn1_5m". as_data_array : bool, optional return the data as as xarray DataArray if true. The default is True. @@ -331,14 +392,38 @@ def get_ahn1(extent, identifier="ahn1_5m", as_data_array=True): @cache.cache_netcdf() -def get_ahn2(extent, identifier="ahn2_5m", as_data_array=True): +def get_ahn2(extent, identifier="AHN2_5M_M", as_data_array=None, **kwargs): """Download AHN2. Parameters ---------- extent : list, tuple or np.array extent - identifier : TYPE, optional + identifier : str, optional + The identifier determines the resolution (05M for 0.5 m and 5M for 5 m) and the + type of height data (M = DTM = surface level, R = DSM = also other features). + Possible values are 'AHN2_05M_I', 'AHN2_05M_N', 'AHN2_05M_R' and 'AHN2_5M_M'. + The default is "AHN2_5M_M". + + Returns + ------- + xr.DataArray + DataArray of the AHN + """ + _assert_as_data_array_is_none(as_data_array) + identifier = _rename_identifier(identifier) + return _get_ahn_ellipsis(extent, identifier, **kwargs) + + +@cache.cache_netcdf() +def get_ahn2_legacy(extent, identifier="ahn2_5m", as_data_array=True): + """Download AHN2. + + Parameters + ---------- + extent : list, tuple or np.array + extent + identifier : str, optional Possible values are 'ahn2_05m_i', 'ahn2_05m_n', 'ahn2_05m_r' and 'ahn2_5m'. The default is "ahn2_5m". as_data_array : bool, optional @@ -355,14 +440,38 @@ def get_ahn2(extent, identifier="ahn2_5m", as_data_array=True): @cache.cache_netcdf() -def get_ahn3(extent, identifier="AHN3_5m_DTM", as_data_array=True): +def get_ahn3(extent, identifier="AHN3_5M_M", as_data_array=None, **kwargs): """Download AHN3. Parameters ---------- extent : list, tuple or np.array extent - identifier : TYPE, optional + identifier : str, optional + The identifier determines the resolution (05M for 0.5 m and 5M for 5 m) and the + type of height data (M = DTM = surface level, R = DSM = also other features). + Possible values are 'AHN3_05M_M', 'AHN3_05M_R', 'AHN3_5M_M' and 'AHN3_5M_R'. + The default is "AHN3_5M_M". + + Returns + ------- + xr.DataArray + DataArray of the AHN + """ + _assert_as_data_array_is_none(as_data_array) + identifier = _rename_identifier(identifier) + return _get_ahn_ellipsis(extent, identifier, **kwargs) + + +@cache.cache_netcdf() +def get_ahn3_legacy(extent, identifier="AHN3_5m_DTM", as_data_array=True): + """Download AHN3. + + Parameters + ---------- + extent : list, tuple or np.array + extent + identifier : str, optional Possible values are 'AHN3_05m_DSM', 'AHN3_05m_DTM', 'AHN3_5m_DSM' and 'AHN3_5m_DTM'. The default is "AHN3_5m_DTM". as_data_array : bool, optional @@ -378,14 +487,38 @@ def get_ahn3(extent, identifier="AHN3_5m_DTM", as_data_array=True): @cache.cache_netcdf() -def get_ahn4(extent, identifier="AHN4_DTM_5m", as_data_array=True): +def get_ahn4(extent, identifier="AHN4_5M_M", as_data_array=None, **kwargs): """Download AHN4. Parameters ---------- extent : list, tuple or np.array extent - identifier : TYPE, optional + identifier : str, optional + The identifier determines the resolution (05M for 0.5 m and 5M for 5 m) and the + type of height data (M = DTM = surface level, R = DSM = also other features). + Possible values are 'AHN4_05M_M', 'AHN4_05M_R', 'AHN4_5M_M' and 'AHN4_5M_R'. + The default is "AHN4_5M_M". + + Returns + ------- + xr.DataArray + DataArray of the AHN + """ + _assert_as_data_array_is_none(as_data_array) + identifier = _rename_identifier(identifier) + return _get_ahn_ellipsis(extent, identifier, **kwargs) + + +@cache.cache_netcdf() +def get_ahn4_legacy(extent, identifier="AHN4_DTM_5m", as_data_array=True): + """Download AHN4. + + Parameters + ---------- + extent : list, tuple or np.array + extent + identifier : str, optional Possible values are 'AHN4_DTM_05m', 'AHN4_DTM_5m', 'AHN4_DSM_05m' and 'AHN4_DSM_5m'. The default is "AHN4_DTM_5m". as_data_array : bool, optional @@ -400,6 +533,67 @@ def get_ahn4(extent, identifier="AHN4_DTM_5m", as_data_array=True): return _download_and_combine_tiles(tiles, identifier, extent, as_data_array) +@cache.cache_netcdf() +def get_ahn5(extent, identifier="AHN5_5M_M", **kwargs): + """Download AHN5. + + Parameters + ---------- + extent : list, tuple or np.array + extent + identifier : str, optional + The identifier determines the resolution (05M for 0.5 m and 5M for 5 m) and the + type of height data (M = DTM = surface level, R = DSM = also other features). + Possible values are 'AHN5_5M_M', 'AHN5_5M_R', 'AHN5_05M_M' and'AHN5_05M_R'. + The default is "AHN5_5M_M". + + Returns + ------- + xr.DataArray + DataArray of the AHN + """ + + return _get_ahn_ellipsis(extent, identifier, **kwargs) + + +@cache.cache_netcdf() +def _get_ahn_ellipsis(extent, identifier="AHN5_5M_M", **kwargs): + """Download AHN5. + + Parameters + ---------- + extent : list, tuple or np.array + extent + identifier : str, optional + Possible values are 'AHN5_5M_M' (dtm), 'AHN5_5M_R' (dsm), 'AHN5_05M_M' (dtm) and + 'AHN5_05M_R' (dsm). The default is "AHN5_5M_M". + + Returns + ------- + xr.DataArray + DataArray of the AHN + """ + tiles = _get_tiles_ellipsis(extent=extent, **kwargs) + if identifier not in tiles.columns: + raise (ValueError(f"Unknown ahn-identifier: {identifier}")) + tiles = tiles[~tiles[identifier].isna()] + das = [] + for tile in tqdm(tiles.index, desc=f"Downloading tiles of {identifier}"): + url = tiles.at[tile, identifier] + if url.endswith(".zip"): + path = url.split("/")[-1].replace(".zip", ".TIF") + if path.lower().endswith(".tif.tif"): + path = path[:-4] + da = rioxarray.open_rasterio( + rasterio.open(f"zip+{url}!/{path}"), mask_and_scale=True + )[0] + else: + da = rioxarray.open_rasterio(url, mask_and_scale=True)[0] + da = da.sel(x=slice(extent[0], extent[1]), y=slice(extent[3], extent[2])) + das.append(da) + return merge_arrays(das) + + def _download_and_combine_tiles(tiles, identifier, extent, as_data_array): """Internal method to download and combine ahn-data.""" if tiles.empty: @@ -423,3 +617,38 @@ def _download_and_combine_tiles(tiles, identifier, extent, as_data_array): da = da.sel(x=slice(extent[0], extent[1]), y=slice(extent[3], extent[2])) return da return memfile + + +def _rename_identifier(identifier): + rename = { + "ahn1_5m": "AHN1_5M", + "ahn2_05m_i": "AHN2_05M_I", + "ahn2_05m_n": "AHN2_05M_N", + "ahn2_05m_r": "AHN2_05M_R", + "ahn2_5m": "AHN2_5M_M", + "AHN3_05m_DTM": "AHN3_05M_M", + "AHN3_05m_DSM": "AHN3_05M_R", + "AHN3_5m_DTM": "AHN3_5M_M", + "AHN3_5m_DSM": "AHN3_5M_R", + "AHN4_DTM_05m": "AHN4_05M_M", + "AHN4_DSM_05m": "AHN4_05M_R", + "AHN4_DTM_5m": "AHN4_5M_M", + "AHN4_DSM_5m": "AHN4_5M_R", + } + if identifier in rename: + id_new = rename[identifier] + logger.warning(f"The identifier {identifier} is deprecated. Rename to {id_new}") + identifier = id_new + return identifier + + +def _assert_as_data_array_is_none(as_data_array): + if as_data_array is not None: + raise ( + DeprecationWarning( + "The as_data_array-argument has been removed from the ahn-" + "methods, and these methods now allways return a DataArray. " + "Remove the as_data_array-argument or use the legcay ahn-" + "methods." + ) + )