forcing: user can set fill_mode and fill_value in ATMInterpolator

thetis/forcing.py

@@ -123,7 +123,9 @@ def __init__(self, function_space, wind_stress_field,
                  ncfile_pattern, init_date, target_coordsys=None,
                  vect_rotator=None,
                  east_wind_var_name='uwind', north_wind_var_name='vwind',
-                 pressure_var_name='prmsl', verbose=False):
+                 pressure_var_name='prmsl', fill_mode=None,
+                 fill_value=numpy.nan,
+                 verbose=False):
         """
         :arg function_space: Target (scalar) :class:`FunctionSpace` object onto
             which data will be interpolated.
@@ -144,14 +146,20 @@ def __init__(self, function_space, wind_stress_field,
             to target function space (optional).
         :kwarg east_wind_var_name, north_wind_var_name, pressure_var_name:
             wind component and pressure field names in netCDF file.
+        :kwarg fill_mode: Determines how points outside the source grid will be
+            treated. If 'nearest', value of the nearest source point will be
+            used. Otherwise a constant fill value will be used (default).
+        :kwarg float fill_value: Set the fill value (default: NaN)
         :kwarg bool verbose: Se True to print debug information.
         """
         self.function_space = function_space
         self.wind_stress_field = wind_stress_field
         self.atm_pressure_field = atm_pressure_field
 
         # construct interpolators
-        self.grid_interpolator = interpolation.NetCDFLatLonInterpolator2d(self.function_space, to_latlon)
+        self.grid_interpolator = interpolation.NetCDFLatLonInterpolator2d(
+            self.function_space, to_latlon, fill_mode=fill_mode,
+            fill_value=fill_value)
         var_list = [east_wind_var_name, north_wind_var_name, pressure_var_name]
         self.reader = interpolation.NetCDFSpatialInterpolator(
             self.grid_interpolator, var_list)

thetis/interpolation.py

@@ -363,11 +363,16 @@ class SpatialInterpolator2d(SpatialInterpolator):
     __metaclass__ = ABCMeta
 
     @PETSc.Log.EventDecorator("thetis.SpatialInterpolator2d.__init__")
-    def __init__(self, function_space, to_latlon):
+    def __init__(self, function_space, to_latlon, fill_mode=None,
+                 fill_value=numpy.nan):
         """
         :arg function_space: target Firedrake FunctionSpace
         :arg to_latlon: Python function that converts local mesh coordinates to
             latitude and longitude: 'lat, lon = to_latlon(x, y)'
+        :kwarg fill_mode: Determines how points outside the source grid will be
+            treated. If 'nearest', value of the nearest source point will be
+            used. Otherwise a constant fill value will be used (default).
+        :kwarg float fill_value: Set the fill value (default: NaN)
         """
         assert function_space.ufl_element().value_shape() == ()
 
@@ -389,6 +394,8 @@ def __init__(self, function_space, to_latlon):
         lat, lon = to_latlon(*coords)
         self.mesh_lonlat = numpy.array([lon, lat]).T
 
+        self.fill_mode = fill_mode
+        self.fill_value = fill_value
         self._initialized = False
 
     @PETSc.Log.EventDecorator("thetis.SpatialInterpolator2d._create_interpolator")
@@ -408,7 +415,9 @@ def _create_interpolator(self, lat_array, lon_array):
         subset_lat = lat_array[self.ind_lon, self.ind_lat].ravel()
         subset_lon = lon_array[self.ind_lon, self.ind_lat].ravel()
         subset_lonlat = numpy.array((subset_lon, subset_lat)).T
-        self.grid_interpolator = GridInterpolator(subset_lonlat, self.mesh_lonlat)
+        self.grid_interpolator = GridInterpolator(
+            subset_lonlat, self.mesh_lonlat, fill_mode=self.fill_mode,
+            fill_value=self.fill_value)
         self._initialized = True
 
         # debug: plot subsets