-
Notifications
You must be signed in to change notification settings - Fork 3
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Buisdrainage #276
Buisdrainage #276
Changes from 9 commits
97bbc0a
8d9f0d1
1b0053f
b98c2da
338b628
bff68d8
c9a2c41
5dc7313
8d9ab6f
3c7609f
c60c2d6
8fc0101
1f7d70a
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
Original file line number | Diff line number | Diff line change |
---|---|---|
|
@@ -13,6 +13,7 @@ | |
rws, | ||
waterboard, | ||
webservices, | ||
nhi, | ||
) | ||
from .geotop import get_geotop | ||
from .regis import get_regis |
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,88 @@ | ||
import os | ||
import requests | ||
import logging | ||
import numpy as np | ||
import rioxarray | ||
from ..dims.resample import structured_da_to_ds | ||
|
||
logger = logging.getLogger(__name__) | ||
|
||
|
||
def download_file(url, pathname, filename=None, overwrite=False, timeout=120): | ||
if filename is None: | ||
filename = url.split("/")[-1] | ||
fname = os.path.join(pathname, filename) | ||
if overwrite or not os.path.isfile(fname): | ||
logger.info(f"Downloading {filename}") | ||
r = requests.get(url, allow_redirects=True, timeout=timeout) | ||
with open(fname, "wb") as file: | ||
file.write(r.content) | ||
return fname | ||
|
||
|
||
def download_buisdrainage(pathname, overwrite=False): | ||
url_bas = "https://thredds.data.nhi.nu/thredds/fileServer/opendap/models/nhi3_2/25m" | ||
|
||
# download resistance | ||
url = f"{url_bas}/buisdrain_c_ras25/buisdrain_c_ras25.nc" | ||
fname_c = download_file(url, pathname, overwrite=overwrite) | ||
|
||
# download drain depth | ||
url = f"{url_bas}/buisdrain_d_ras25/buisdrain_d_ras25.nc" | ||
fname_d = download_file(url, pathname, overwrite=overwrite) | ||
|
||
return fname_c, fname_d | ||
|
||
|
||
def add_buisdrainage( | ||
ds, | ||
pathname=None, | ||
cond_var="buisdrain_cond", | ||
depth_var="buisdrain_depth", | ||
cond_method="average", | ||
depth_method="mode", | ||
): | ||
if pathname is None: | ||
pathname = ds.cachedir | ||
# download files if needed | ||
fname_c, fname_d = download_buisdrainage(pathname) | ||
|
||
# make sure crs is set of ds | ||
if ds.rio.crs is None: | ||
ds = ds.rio.write_crs(28992) | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. should we use the EPSG definition in nlmod for this? There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Somehow this does not give the same problems as with maps (I believe, but did not thoroughly check). It might break stuff, if we change this. If other data sources are also in 28992, they do not need to be transformed, but if we change the crs of the model Dataset to our own definition, an unwanted transformation is applied. The buisdrainage data is in WGS84, so you may be right. Let's check this in the future. |
||
|
||
# use method = "average" on conductance, to account for locations without pipe | ||
# drainage, where the conductance is 0 | ||
buisdrain_c = rioxarray.open_rasterio(fname_c, mask_and_scale=True)[0] | ||
# calculate a conductance (per m2) from a resistance | ||
cond = 1 / buisdrain_c | ||
# set conductance to 0 where resistance is infinite or 0 | ||
cond = cond.where(~(np.isinf(cond) | np.isnan(cond)), 0.0) | ||
cond = cond.rio.write_crs(buisdrain_c.rio.crs) | ||
# resample to model grid | ||
ds[cond_var] = structured_da_to_ds(cond, ds, method=cond_method) | ||
# multiply by area to get a conductance | ||
ds[cond_var] = ds[cond_var] * ds["area"] | ||
|
||
# use mthod = "mode" on depth, to retreive the depth that occurs most in each cell | ||
mask_and_scale = False | ||
buisdrain_d = rioxarray.open_rasterio(fname_d, mask_and_scale=mask_and_scale)[0] | ||
if mask_and_scale: | ||
nodata = np.nan | ||
else: | ||
nodata = buisdrain_d.attrs["_FillValue"] | ||
# set buisdrain_d to nodata where it is 0 | ||
mask = buisdrain_d != 0 | ||
buisdrain_d = buisdrain_d.where(mask, nodata).rio.write_crs(buisdrain_d.rio.crs) | ||
# resample to model grid | ||
ds[depth_var] = structured_da_to_ds( | ||
buisdrain_d, ds, method=depth_method, nodata=nodata | ||
) | ||
if not mask_and_scale: | ||
# set nodata values to NaN | ||
ds[depth_var] = ds[depth_var].where(ds[depth_var] != nodata) | ||
|
||
# from cm to m | ||
ds[depth_var] = ds[depth_var] / 100.0 | ||
|
||
return ds |
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,22 @@ | ||
import os | ||
import numpy as np | ||
import tempfile | ||
import nlmod | ||
import pytest | ||
|
||
tmpdir = tempfile.gettempdir() | ||
|
||
|
||
@pytest.mark.slow | ||
def test_buidrainage(): | ||
model_ws = os.path.join(tmpdir, "buidrain") | ||
ds = nlmod.get_ds([110_000, 130_000, 435_000, 445_000], model_ws=model_ws) | ||
ds = nlmod.read.nhi.add_buisdrainage(ds) | ||
|
||
# assert that all locations with a specified depth also have a positive conductance | ||
mask = ~ds["buisdrain_depth"].isnull() | ||
assert np.all(ds["buisdrain_cond"].data[mask] > 0) | ||
|
||
# assert that all locations with a positive conductance also have a specified depth | ||
mask = ds["buisdrain_cond"] > 0 | ||
assert np.all(~np.isnan(ds["buisdrain_depth"].data[mask])) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Why set seasons to None? The keyword summer_months sort of suggests winter/summer being the only possible two seasons?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I changed this as a tuple gave a bug,
seasons
(time_series_namerecord
) somehow needed to be a list, and it is not allowed to to add a list as an optional arggument. Therefore whenseasons
is None, I set it toseasons = ["winter", "summer"]
. Maybe better to add two arguments, calledwinter_name
andsummer_name
?