927 fix incorrect number of returned points by dfmtinterp uds to plip…

…oints with outofbounds plipoints (#948)

* updated minimal xugrid version

* updated interp_uds_to_plipoints by removing error

* passing non-topology vars

* updated/improved test

* updated example file

* updated whatsnew

dfm_tools/interpolate_grid2bnd.py

@@ -465,25 +465,30 @@ def interp_uds_to_plipoints(uds:xu.UgridDataset, gdf:geopandas.GeoDataFrame) ->
     # TODO: this function requires gdf with points, but the interp_regularnc_to_plipoints requires paths to plifiles (and others also)
 
     facedim = uds.grid.face_dimension
+    edgedim = uds.grid.edge_dimension
+    nodedim = uds.grid.node_dimension
     ncbnd_construct = get_ncbnd_construct()
     dimn_point = ncbnd_construct['dimn_point']
     varn_pointname = ncbnd_construct['varn_pointname']
 
-    ds = uds.ugrid.sel_points(x=gdf.geometry.x, y=gdf.geometry.y)
-    #TODO: drop mesh2d_face_x and mesh2d_face_y variables
-
-    if len(gdf)!=ds.sizes[facedim]: #TODO: check this until https://github.com/Deltares/xugrid/issues/100 is solved, after that, make a testcase that checks only this if-statement
-        ds_points = geopandas.points_from_xy(ds.x,ds.y)
-        gdfpoint_inds_bool = pd.Series(index=range(len(gdf)))
-        gdfpoint_inds_bool[:] = True
-        for iR, gdf_row in gdf.iterrows():
-            gdf_point = gdf_row.geometry
-            if gdf_point in ds_points:
-                gdfpoint_inds_bool.iloc[iR] = False
-        gdf_stats = gdf.copy()
-        gdf_stats['missing'] = gdfpoint_inds_bool
-        raise ValueError(f'requested {len(gdf)} points but resulted in ds with {ds.sizes[facedim]} points, missing points are probably outside of the uds model domain:\n{gdf_stats}')
-
+    # drop node/edge variables since they are not interpolated but 
+    # get an additional face dimension if some points are out of bounds
+    # TODO: revert after fixing https://github.com/Deltares/xugrid/issues/274
+    vars_without_facedim = []
+    for varn in uds.variables:
+        if facedim not in uds[varn].dims:
+            vars_without_facedim.append(varn)
+    uds_face = uds.drop(vars_without_facedim)
+
+    # interpolate to provided points
+    ds = uds_face.ugrid.sel_points(x=gdf.geometry.x, y=gdf.geometry.y)
+
+    # re-add removed variables again, sometimes important for e.g. depth
+    # TODO: remove after fixing https://github.com/Deltares/xugrid/issues/274
+    for varn in vars_without_facedim:
+        if edgedim not in uds[varn].dims and nodedim not in uds[varn].dims:
+            ds[varn] = uds[varn]
+
     # rename station dimname and varname (is index, are both mesh2d_nFaces to start with)
     ds = ds.rename({facedim:dimn_point}) # rename mesh2d_nFaces to plipoints
     ds = ds.rename_vars({dimn_point:varn_pointname})

docs/whats-new.md

@@ -8,6 +8,7 @@
 - update private functions under `dfmt.download_ERA5()` to CDS-Beta (requires ECMWF apikey instead) in [#925](https://github.com/Deltares/dfm_tools/pull/925)
 - simplified prevention of int dtypes in `dfmt.preprocess_ERA5()` in [#943](https://github.com/Deltares/dfm_tools/pull/943)
 - simplified `dfmt.open_dataset_extra()` by dropping multi-quantity support in [#946](https://github.com/Deltares/dfm_tools/pull/946)
+- improved `dfmt.interp_uds_to_plipoints()` by supporting outofbound points and new xugrid version in [#948](https://github.com/Deltares/dfm_tools/pull/948)
 
 ### Fix
 - also apply convert_360to180 to longitude variable in `dfmt.open_dataset_curvilinear()` in [#913](https://github.com/Deltares/dfm_tools/pull/913)

pyproject.toml

@@ -36,8 +36,8 @@ dependencies = [
 	"netcdf4>=1.5.4",
 	#bottleneck<1.3.3 pip install fails in py39
 	"bottleneck>=1.3.3",
-	#xugrid<0.11.0 does not use all meshkernel attrs (https://github.com/Deltares/xugrid/pull/250) and is not rechunking partitions (https://github.com/Deltares/xugrid/pull/253)
-	"xugrid>=0.11.0",
+	#xugrid<0.11.1 drops outofbounds plipoints in uds.sel_points() https://github.com/Deltares/xugrid/issues/100
+	"xugrid>=0.11.1",
 	#cdsapi<0.7.0 supports credentials from environment variables
 	"cdsapi>=0.7.0",
 	#pydap<3.4.0 is from May 2017 and does not support newer python versions

tests/examples/postprocess_interpolate_uds_toplipoints.py

@@ -42,7 +42,7 @@
 ds_atdepths = ds_atdepths.rename({'depth_from_z0':'depth'})
 
 #interpolate to plipoints
-ds_plipoints = dfmt.interp_uds_to_plipoints(uds=ds_atdepths, gdf=gdf.drop(4)) #gdf.drop(4) #workaround for plipoints out of the model domain
+ds_plipoints = dfmt.interp_uds_to_plipoints(uds=ds_atdepths, gdf=gdf)
 fig,ax = plt.subplots()
 uds.mesh2d_flowelem_bl.ugrid.plot(ax=ax,center=False)
 gdf.plot(ax=ax)

tests/test_interpolate_grid2bnd.py

@@ -529,8 +529,9 @@ def test_interp_uds_to_plipoints():
     file_nc = dfmt.data.fm_grevelingen_map(return_filepath=True)
     uds = dfmt.open_partitioned_dataset(file_nc)
 
-    gdf = gpd.GeoDataFrame(geometry=gpd.points_from_xy([51500,55000],[418800,421500]))
-    gdf[varn_pointname] = ['pt_0001','pt_0002']
+    # adding xy=[1,1] deliberately to test if there are nans included
+    gdf = gpd.GeoDataFrame(geometry=gpd.points_from_xy([51500,55000,1],[418800,421500,1]))
+    gdf[varn_pointname] = ['pt_0001','pt_0002','pt_0003']
 
     # fig, ax = plt.subplots()
     # uds.mesh2d_flowelem_bl.ugrid.plot(ax=ax)
@@ -543,12 +544,14 @@ def test_interp_uds_to_plipoints():
     #interpolate to plipoints
     ds_plipoints = dfmt.interp_uds_to_plipoints(uds=ds_atdepths, gdf=gdf) #workaround for plipoints out of the model domain
     assert varn_pointname in ds_plipoints.coords
+    assert 'depth' in ds_plipoints.coords # TODO: maybe should be z like ncbnd_construct['varn_depth'], although get_Dataset_atdepths() makes depth variable
 
     retrieved = ds_plipoints.mesh2d_sa1.isel(time=-1).to_numpy()
     expected = np.array([[29.00111981, 28.99379263],
-                         [28.95170139, 28.63716083]])
+                         [28.95170139, 28.63716083],
+                         [np.nan, np.nan]])
 
-    assert (np.abs(retrieved - expected) < 1e-8).all()
+    assert np.allclose(retrieved, expected, equal_nan=True)
 
 
 def test_ext_add_boundary_object_per_polyline_onepolyline(tmp_path):