628 skip plipoints containing all-nans when converting to forcingmodel

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)

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','')

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)

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
+ assert forcingmodel_object.forcing[0].name == 'abc_bnd_0003'
+ assert forcingmodel_object.forcing[1].name == 'abc_bnd_0004'
+
+
 @pytest.mark.systemtest
 def test_open_dataset_extra_correctdepths():
  """