Set points outside grid to undefined (#56)

* Set points outside grid to undefined

* Linting

* Correct unit test

* Linting again

* Linting again

* Linting again

* Linting again

* Lint

Author: trygveasp

pysurfex/interpolation.py

@@ -326,7 +326,14 @@ def gridpos2points(
     )
 
 
-def grid2points(grid, points, grid_values, operator="bilinear", elev_gradient=None):
+def grid2points(
+    grid,
+    points,
+    grid_values,
+    operator="bilinear",
+    elev_gradient=None,
+    max_distance=25000.0,
+):
     """Convert a grid to points.
 
     Args:
@@ -335,6 +342,8 @@ def grid2points(grid, points, grid_values, operator="bilinear", elev_gradient=No
         grid_values (np.ndarray): Grid values
         operator (str, optional): Interpolation operator. Defaults to "bilinear".
         elev_gradient (float, optional): Elevation gradient for downscaler
+        max_distance (float, optional): Maximum distance from grid points.
+                                        Defaults to 25000.0.
 
     Raises:
         NotImplementedError: Operator not implemented
@@ -355,20 +364,24 @@ def grid2points(grid, points, grid_values, operator="bilinear", elev_gradient=No
             values = gridpp.simple_gradient(
                 grid.grid, points.points, grid_values, elev_gradient, gridpp.Bilinear
             )
+        in_grid = np.array(points.inside_grid(grid, distance=max_distance))
+        values[in_grid == False] = np.nan  # noqa: E712
     elif operator == "nearest":
         if elev_gradient is None:
             values = gridpp.nearest(grid.grid, points.points, grid_values)
         else:
             values = gridpp.simple_gradient(
                 grid.grid, points.points, grid_values, elev_gradient, gridpp.Nearest
             )
+        in_grid = np.array(points.inside_grid(grid, distance=max_distance))
+        values[in_grid == False] = np.nan  # noqa: E712
     else:
         raise NotImplementedError(f"Operator {operator} not implemented!")
     return values
 
 
 def inside_grid(grid_lons, grid_lats, p_lons, p_lats, distance=2500.0):
-    """Get number of neighbours.
+    """Check if inside grid.
 
     Args:
         grid_lons (np.ndarray): Grid longitudes

tests/unit/test_verification.py

@@ -430,28 +430,28 @@ def fixture_data_surfex_nc_file(tmp_path_factory):
 
  BETA = 0 ;
 
- LATORI = 59.1044427299632 ;
+ LATORI = 59.0 ;
 
- LONORI = 9.09806611390404 ;
+ LONORI = 9.1 ;
 
  IMAX = 2 ;
 
  JMAX = 3 ;
 
   XX =
-  10000, 20000, 30000, 40000, 50000, 60000;
+  40000, 80000, 40000, 80000, 40000, 80000;
 
   YY =
-  10000, 20000, 30000, 40000, 50000, 60000;
+  0, 40000, 80000, 0, 40000, 80000;
 
   DX =
-  10000, 10000, 10000, 10000, 10000, 10000;
+  40000, 40000, 40000, 40000, 40000, 40000;
 
   DY =
-  10000, 10000, 10000, 10000, 10000, 10000;
+  40000, 40000, 40000, 40000, 40000, 40000;
 
   T2M =
-  285, 285, 285, 285, 285, 285;
+  281, 282, 283, 284, 285, 286;
 
   HU2M =
   0.1, 0.2, 0.3, 0.4, 0.5, 0.6;
@@ -488,7 +488,7 @@ def test_verif_nc(tmp_path_factory, stationlist_file, data_surfex_nc_file):
     with working_directory(tmp_path_factory.getbasetemp()):
         converter2ds(argv=argv)
         ds = xr.open_dataset(vfilename, engine="netcdf4")
-        assert ds.fcst.data[0][0][1] == 285
+        assert ds.fcst.data[0][0][1] == 283.0
 
 
 @pytest.mark.usefixtures("_mockers")