From b54557623d6d1a62003b7ed7ca4601adf59b5146 Mon Sep 17 00:00:00 2001 From: veenstrajelmer Date: Fri, 3 Nov 2023 20:09:33 +0100 Subject: [PATCH 1/5] skip all-nan plipoint --- dfm_tools/interpolate_grid2bnd.py | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/dfm_tools/interpolate_grid2bnd.py b/dfm_tools/interpolate_grid2bnd.py index 865930531..1fdd37434 100644 --- a/dfm_tools/interpolate_grid2bnd.py +++ b/dfm_tools/interpolate_grid2bnd.py @@ -570,11 +570,16 @@ def plipointsDataset_to_ForcingModel(plipointsDataset): #select data for this point, ffill nans, concatenating time column, constructing T3D/TimeSeries and append to hcdfm.ForcingModel() datablock_xr_onepoint = plipointsDataset.isel({dimn_point:iP}) plipoint_name = str(datablock_xr_onepoint[varn_pointname].to_numpy()) + plipoint_onlynan = False for quan in quantity_list: datablock_xr_onepoint[quan].attrs['locationname'] = plipoint_name #TODO: is there a nicer way of passing this data? if datablock_xr_onepoint[quan].isnull().all(): # check if all values of plipoint are nan (on land) - warnings.warn(UserWarning(f'Plipoint "{plipoint_name}" might be on land since it only contain nan values. Consider altering your plifile or using plipointsDataset.ffill(dim="plipoints").bfill(dim="plipoints"). Nan values replaced with 0 to avoid bc-writing errors')) #TODO: maybe fill along plipoints dimension, but beware on nans in deep water that are filled with neighbours - datablock_xr_onepoint[quan] = datablock_xr_onepoint[quan].fillna(0) + plipoint_onlynan = True + warnings.warn(UserWarning(f'Plipoint "{plipoint_name}" might be on land since it only contain nan values. This point is skipped to avoid bc-writing errors. Consider altering your PolyFile or extrapolate the data.')) + + # skip this point if quantity has only-nan values + if plipoint_onlynan: + continue if dimn_depth in plipointsDataset.dims: ts_one = Dataset_to_T3D(datablock_xr_onepoint) From 844c629552f7d538830eec593c404f2142514e6a Mon Sep 17 00:00:00 2001 From: veenstrajelmer Date: Fri, 3 Nov 2023 21:07:01 +0100 Subject: [PATCH 2/5] added unittest --- dfm_tools/modelbuilder.py | 2 +- tests/test_interpolate_grid2bnd.py | 47 +++++++++++++++++++++++++++++- 2 files changed, 47 insertions(+), 2 deletions(-) diff --git a/dfm_tools/modelbuilder.py b/dfm_tools/modelbuilder.py index cf8590fbb..41843d242 100644 --- a/dfm_tools/modelbuilder.py +++ b/dfm_tools/modelbuilder.py @@ -28,7 +28,7 @@ ] -def cmems_nc_to_bc(ext_bnd, list_quantities, tstart, tstop, file_pli, dir_pattern, dir_output, refdate_str): +def cmems_nc_to_bc(ext_bnd, list_quantities, tstart, tstop, file_pli, dir_pattern, dir_output, refdate_str=None): #input examples in https://github.com/Deltares/dfm_tools/blob/main/tests/examples/preprocess_interpolate_nc_to_bc.py file_bc_basename = os.path.basename(file_pli).replace('.pli','') diff --git a/tests/test_interpolate_grid2bnd.py b/tests/test_interpolate_grid2bnd.py index 764cb4ac5..cc6c85d89 100644 --- a/tests/test_interpolate_grid2bnd.py +++ b/tests/test_interpolate_grid2bnd.py @@ -12,12 +12,15 @@ import datetime as dt import xarray as xr import shapely +import pandas as pd import geopandas as gpd from dfm_tools.interpolate_grid2bnd import (read_polyfile_as_gdf_points, tidemodel_componentlist, components_translate_upper, - get_ncbnd_construct + get_ncbnd_construct, + interp_regularnc_to_plipointsDataset, ) +from dfm_tools.hydrolib_helpers import PolyFile_to_geodataframe_points import hydrolib.core.dflowfm as hcdfm @@ -167,6 +170,48 @@ def test_interpolate_nc_to_bc(): assert forcing0.quantityunitpair[1].quantity == 'salinitybnd' +@pytest.mark.systemtest +def test_plipointsDataset_to_ForcingModel_drop_allnan_points(): + #construct polyfile gdf + point_x = [-71.5, -71.5, -71.5, -71.5,] + point_y = [12.5, 12.6, 12.7, 12.8] + polyobject_pd = pd.DataFrame(dict(x=point_x, y=point_y)) + poly_object = dfmt.DataFrame_to_PolyObject(polyobject_pd, name="abc_bnd") + polyfile_object = hcdfm.PolyFile() + polyfile_object.objects.append(poly_object) + gdf_points = PolyFile_to_geodataframe_points(polyfile_object) + + # actual cmems data + no3_values = [[[ np.nan, np.nan, np.nan, np.nan, + 5.2622618e-05], + [ np.nan, np.nan, np.nan, np.nan, + 5.1448391e-05], + [9.0356334e-05, 1.8063841e-04, np.nan, 7.0045113e-05, + 4.9546972e-05], + [4.0657567e-05, 4.6050434e-05, 4.7762936e-05, 4.2734890e-05, + 4.0186114e-05], + [2.6492798e-05, 2.8967228e-05, 3.0298394e-05, 3.1432221e-05, + 3.2138258e-05], + [2.4448847e-05, 2.6705173e-05, 2.6929094e-05, 2.5548252e-05, + 2.5827681e-05]]] + ds = xr.Dataset() + ds['longitude'] = xr.DataArray([-72. , -71.75, -71.5 , -71.25, -71. ], dims='longitude') + ds['latitude'] = xr.DataArray([12. , 12.25, 12.5 , 12.75, 13. , 13.25], dims='latitude') + ds['time'] = xr.DataArray([639204.],dims='time').assign_attrs({'units': 'hours since 1950-01-01'}) + ds['tracerbndNO3'] = xr.DataArray(no3_values, dims=('time','latitude','longitude')).assign_attrs({'units': 'mmol m-3'}) + ds = xr.decode_cf(ds) + + # interpolate to points and convert to forcingmodel + data_interp = interp_regularnc_to_plipointsDataset(data_xr_reg=ds, gdf_points=gdf_points) + ForcingModel_object = dfmt.plipointsDataset_to_ForcingModel(plipointsDataset=data_interp) + + # check if the resulting forcingmodel does not have 4 but 2 points + # the first two points were skipped because they were all nan + assert len(ForcingModel_object.forcing) == 2 + ForcingModel_object.forcing[0].name == 'abc_bnd_0003' + ForcingModel_object.forcing[1].name == 'abc_bnd_0004' + + @pytest.mark.systemtest def test_open_dataset_extra_correctdepths(): """ From b9911ca7a114a3c889feccd45701f17cb970e962 Mon Sep 17 00:00:00 2001 From: veenstrajelmer Date: Fri, 3 Nov 2023 21:11:32 +0100 Subject: [PATCH 3/5] updated whatsnew --- docs/whats-new.md | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/whats-new.md b/docs/whats-new.md index 8745731cb..8f764ba21 100644 --- a/docs/whats-new.md +++ b/docs/whats-new.md @@ -12,6 +12,7 @@ - increased buffer in `dfmt.download_ERA5()` by [@veenstrajelmer](https://github.com/veenstrajelmer) in [#612](https://github.com/Deltares/dfm_tools/pull/612) - support for Polygon geometries in `dfmt.geodataframe_to_PolyFile()` by [@veenstrajelmer](https://github.com/veenstrajelmer) in [#610](https://github.com/Deltares/dfm_tools/pull/610) - fill nan-values in initial salinity/temperature netcdf dataset in `dfmt.preprocess_ini_cmems_to_nc()` by [@veenstrajelmer](https://github.com/veenstrajelmer) in [#617](https://github.com/Deltares/dfm_tools/pull/617) +- skip all-nan boundary support points instead of converting to zeros in `dfmt.plipointsDataset_to_ForcingModel()` by [@veenstrajelmer](https://github.com/veenstrajelmer) in [#637](https://github.com/Deltares/dfm_tools/pull/637) ## 0.15.0 (2023-10-19) From b35978587af63ea6c601112722f12280f275239f Mon Sep 17 00:00:00 2001 From: veenstrajelmer Date: Fri, 3 Nov 2023 21:12:32 +0100 Subject: [PATCH 4/5] fixed code bugs --- tests/test_interpolate_grid2bnd.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/test_interpolate_grid2bnd.py b/tests/test_interpolate_grid2bnd.py index cc6c85d89..c6f751883 100644 --- a/tests/test_interpolate_grid2bnd.py +++ b/tests/test_interpolate_grid2bnd.py @@ -208,8 +208,8 @@ def test_plipointsDataset_to_ForcingModel_drop_allnan_points(): # check if the resulting forcingmodel does not have 4 but 2 points # the first two points were skipped because they were all nan assert len(ForcingModel_object.forcing) == 2 - ForcingModel_object.forcing[0].name == 'abc_bnd_0003' - ForcingModel_object.forcing[1].name == 'abc_bnd_0004' + assert ForcingModel_object.forcing[0].name == 'abc_bnd_0003' + assert ForcingModel_object.forcing[1].name == 'abc_bnd_0004' @pytest.mark.systemtest From 6e381805caf0890b25574a18791eab8d20a43339 Mon Sep 17 00:00:00 2001 From: veenstrajelmer Date: Fri, 3 Nov 2023 21:13:06 +0100 Subject: [PATCH 5/5] fixed code smells --- tests/test_interpolate_grid2bnd.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/tests/test_interpolate_grid2bnd.py b/tests/test_interpolate_grid2bnd.py index c6f751883..062f66216 100644 --- a/tests/test_interpolate_grid2bnd.py +++ b/tests/test_interpolate_grid2bnd.py @@ -203,13 +203,13 @@ def test_plipointsDataset_to_ForcingModel_drop_allnan_points(): # interpolate to points and convert to forcingmodel data_interp = interp_regularnc_to_plipointsDataset(data_xr_reg=ds, gdf_points=gdf_points) - ForcingModel_object = dfmt.plipointsDataset_to_ForcingModel(plipointsDataset=data_interp) + forcingmodel_object = dfmt.plipointsDataset_to_ForcingModel(plipointsDataset=data_interp) # check if the resulting forcingmodel does not have 4 but 2 points # the first two points were skipped because they were all nan - assert len(ForcingModel_object.forcing) == 2 - assert ForcingModel_object.forcing[0].name == 'abc_bnd_0003' - assert ForcingModel_object.forcing[1].name == 'abc_bnd_0004' + assert len(forcingmodel_object.forcing) == 2 + assert forcingmodel_object.forcing[0].name == 'abc_bnd_0003' + assert forcingmodel_object.forcing[1].name == 'abc_bnd_0004' @pytest.mark.systemtest