4/8 Feature GISA-new impervious surface

pyVPRM/VPRM.py

@@ -1,3 +1,5 @@
+"""Preprocess satellite, land-cover, and static ancillary inputs for VPRM."""
+
 import warnings
 
 warnings.filterwarnings("ignore")
@@ -15,6 +17,7 @@
     to_esmf_grid,
     replace_inf_runs_ignore_nans
 )
+from pyVPRM.lib.regridding import conservative_regrid
 from scipy.ndimage import uniform_filter
 from pyproj import Transformer
 import copy
@@ -39,8 +42,20 @@
 
 
 class vprm_preprocessor:
-    """
-    Main class for the  Vegetation Photosynthesis and Respiration Model
+    """Prepare satellite and static spatial inputs for VPRM calculations.
+
+    Parameters
+    ----------
+    vprm_config_path : str or pathlib.Path
+        VPRM vegetation-class configuration file.
+    land_cover_map : satellite_data_manager, optional
+        Pre-calculated VPRM land-cover fractions.
+    verbose : bool, default=False
+        Enable verbose processing output.
+    n_cpus : int, default=1
+        Number of CPUs available to regridding operations.
+    flux_tower_instances : list, optional
+        Flux-tower data objects used for fitting workflows.
     """
 
     def __init__(
@@ -79,6 +94,7 @@ def __init__(
         self.empty_xr_lat_lon_grid = None
 
         self.land_cover_type = land_cover_map
+        self.impervious_surface_area = None
         # land_cover_type: tmin, topt, tmax
 
         with open(vprm_config_path, "r") as stream:
@@ -724,6 +740,88 @@ def add_land_cover_map(
                 self.land_cover_type.save(save_path)
         return
 
+    def add_impervious_surface_area(
+        self,
+        impervious_surface_area,
+        var_name="impervious_surface_percentage",
+        regridder_save_path=None,
+        n_cpus=None,
+        mpi=True,
+        logs=False,
+        max_source_cells=None,
+    ):
+        """Conservatively regrid an impervious-surface field to the satellite grid.
+
+        Parameters
+        ----------
+        impervious_surface_area : satellite_data_manager
+            Loaded and spatially cropped continuous impervious-surface data.
+            ``var_name`` must hold percentages from 0 to 100.
+        var_name : str, default="impervious_surface_percentage"
+            Name of the impervious-surface percentage variable.
+        regridder_save_path : str or pathlib.Path
+            Grid-specific path for conservative ESMF regridding weights.
+        n_cpus : int, optional
+            Number of MPI processes used while generating weights. Defaults to
+            :attr:`n_cpus`.
+        mpi : bool, default=True
+            Use ``mpirun`` when generating new ESMF weights.
+        logs : bool, default=False
+            Preserve ESMF log files when generating weights.
+        max_source_cells : int, optional
+            Maximum source cells in one ESMF invocation. Larger regular source
+            grids are split into disjoint y-axis partitions and their
+            conservative destination contributions are summed.
+
+        Returns
+        -------
+        None
+            The regridded percentage field is stored in
+            :attr:`impervious_surface_area` on the VPRM satellite grid.
+
+        Raises
+        ------
+        TypeError
+            If the supplied data are not managed by
+            :class:`satellite_data_manager`.
+        ValueError
+            If the requested variable or weight-cache path is missing.
+        RuntimeError
+            If ESMF weight generation fails.
+        """
+        if not isinstance(impervious_surface_area, satellite_data_manager):
+            raise TypeError(
+                "impervious_surface_area must be a satellite_data_manager instance."
+            )
+        if var_name not in impervious_surface_area.sat_img:
+            raise ValueError(
+                "Impervious-surface data do not contain {!r}.".format(var_name)
+            )
+        if regridder_save_path is None:
+            raise ValueError("regridder_save_path is required for ISA regridding.")
+
+        if n_cpus is None:
+            n_cpus = self.n_cpus
+        regridded = conservative_regrid(
+            impervious_surface_area.sat_img[[var_name]],
+            self.sat_imgs.sat_img,
+            weight_path=regridder_save_path,
+            source_mask=impervious_surface_area.sat_img[var_name] > 0,
+            max_source_cells=max_source_cells,
+            n_cpus=n_cpus,
+            mpi=mpi,
+            logs=logs,
+        )
+        regridded[var_name].attrs.update(
+            {
+                "long_name": "area-weighted percent impervious surface",
+                "units": "percent",
+                "source_product": type(impervious_surface_area).__name__,
+            }
+        )
+        self.impervious_surface_area = satellite_data_manager(sat_img=regridded)
+        return
+
     def calc_min_max_evi_lswi(self):
         """
         Calculate the minimim and maximum EVI and LSWI
@@ -1178,15 +1276,21 @@ def is_disjoint(self, this_sat_img):
         return dj
 
     def add_vprm_insts(self, vprm_insts, allow_reproject=True):
-        """
-        Merge one or more other vprm_preprocessor instances' satellite/land-cover
-        tiles into this one (e.g. combining adjacent spatial tiles into a
-        single larger domain).
-
-        Parameters:
-                vprm_insts (list): other vprm_preprocessor instances to merge in
-                allow_reproject (bool): reproject incoming tiles onto this
-                                        instance's grid/CRS if they don't already match
+        """Merge compatible VPRM preprocessors, including static ISA fields.
+
+        Parameters
+        ----------
+        vprm_insts : list[vprm_preprocessor]
+            Preprocessors whose spatial products will be tiled onto this
+            instance.
+        allow_reproject : bool, default=True
+            Reproject satellite data before merging when necessary.
+
+        Returns
+        -------
+        None
+            Satellite, land-cover, and impervious-surface products are merged
+            in place when present.
         """
         if isinstance(self.sat_imgs, satellite_data_manager):
             self.sat_imgs.add_tile(
@@ -1201,3 +1305,12 @@ def add_vprm_insts(self, vprm_insts, allow_reproject=True):
             self.land_cover_type.add_tile(
                 [v.land_cover_type for v in vprm_insts], reproject=False
             )
+        if self.impervious_surface_area is not None:
+            self.impervious_surface_area.add_tile(
+                [
+                    v.impervious_surface_area
+                    for v in vprm_insts
+                    if v.impervious_surface_area is not None
+                ],
+                reproject=False,
+            )

pyVPRM/lib/functions.py

@@ -247,7 +247,31 @@ def make_xesmf_grid(satellite_image, transformer = None):
     return pixel_grid
 
 
-def to_esmf_grid(sat_img):
+def to_esmf_grid(sat_img, grid_imask=None):
+    """Convert a regular grid to the SCRIP dataset consumed by ESMF.
+
+    Parameters
+    ----------
+    sat_img : xarray.Dataset or dict
+        Regular source or destination grid. Xarray inputs require
+        one-dimensional ``x`` and ``y`` coordinates and a rio CRS. Dictionary
+        inputs require ``lons`` and ``lats`` arrays.
+    grid_imask : xarray.DataArray or numpy.ndarray, optional
+        Integer or boolean ``y, x`` mask. Values of zero exclude cells from
+        ESMF weight generation; non-zero values retain them. If omitted, all
+        cells participate.
+
+    Returns
+    -------
+    xarray.Dataset
+        SCRIP-format grid description with centre coordinates, cell corners,
+        and ``grid_imask``.
+
+    Raises
+    ------
+    ValueError
+        If a supplied mask does not match the grid dimensions.
+    """
 
     if isinstance(sat_img, dict):
         x = sat_img["lons"]
@@ -292,7 +316,18 @@ def to_esmf_grid(sat_img):
         )
         x_center, y_center = t.transform(x_center, y_center)
         x_corner, y_corner = t.transform(x_corner, y_corner)
-    grid_imask = np.ones((ny, nx), dtype=np.int32)
+    if grid_imask is None:
+        grid_imask = np.ones((ny, nx), dtype=np.int32)
+    else:
+        if hasattr(grid_imask, "values"):
+            grid_imask = grid_imask.values
+        grid_imask = np.asarray(grid_imask, dtype=np.int32)
+        if grid_imask.shape != (ny, nx):
+            raise ValueError(
+                "grid_imask shape {} does not match y, x dimensions ({}, {}).".format(
+                    grid_imask.shape, ny, nx
+                )
+            )
 
     # generate output dataset
     dso = xr.Dataset()
@@ -877,4 +912,3 @@ def format_lon(lon):
 
     return float(chm_box.mean())
 
-