Merge branch 'coli' of https://github.com/C2SM-RCM/emiproc into coli

emiproc/grids.py

@@ -4,18 +4,20 @@
 grids used in different emissions inventories.
 """
 from __future__ import annotations
+
+import math
+import warnings
 from functools import cache, cached_property
 from typing import Iterable
-import warnings
-import numpy as np
-import xarray as xr
+
 import geopandas as gpd
+import numpy as np
 import pyproj
-import math
-
+import xarray as xr
 from netCDF4 import Dataset
-from shapely.geometry import Polygon, Point, box, LineString, MultiPolygon
+from shapely.geometry import LineString, MultiPolygon, Point, Polygon, box
 from shapely.ops import split
+from shapely.creation import polygons
 
 WGS84 = 4326
 WGS84_PROJECTED = 3857
@@ -217,13 +219,37 @@ def __init__(
         self.nx, self.ny = nx, ny
         self.dx, self.dy = (xmax - xmin) / nx, (ymax - ymin) / ny
 
-        self.lon_range = np.linspace(xmin, xmax, nx) + self.dx / 2
-        self.lat_range = np.linspace(ymin, ymax, ny) + self.dy / 2
+        self.lon_range = np.arange(xmin, xmax, self.dx) + self.dx / 2
+        self.lat_range = np.arange(ymin, ymax, self.dy) + self.dy / 2
 
         if name is None:
             name = f"RegularGrid_x({xmin},{xmax})_y({ymin},{ymax})_nx({nx})_ny({ny})"
 
         super().__init__(name, crs)
+
+    @cached_property
+    def cells_as_polylist(self) -> list[Polygon]:
+
+        x_coords, y_coords = np.meshgrid(
+            self.lon_range - self.dx / 2.,
+            self.lat_range - self.dy / 2.
+        )
+        # Reshape to 1D (order set for backward compatibility)
+        x_coords = x_coords.flatten(order="F")
+        y_coords = y_coords.flatten(order="F")
+        dx = float(self.dx)
+        dy = float(self.dy)
+        coords = np.array(
+            [
+                [x, y]
+                for x, y in zip(
+                    [x_coords, x_coords, x_coords + dx, x_coords + dx],
+                    [y_coords, y_coords + dy, y_coords + dy, y_coords],
+                )
+            ]
+        )
+        coords = np.rollaxis(coords, -1, 0)
+        return polygons(coords)
 
     def cell_corners(self, i, j):
         """Return the corners of the cell with indices (i,j).

emiproc/profiles/vprm.py

@@ -74,7 +74,7 @@ def calculate_vegetation_indices(
 
 
 def calculate_vprm_emissions(df: pd.DataFrame, df_vprm: pd.DataFrame) -> pd.DataFrame:
-    """Caluculate the emissions using the VPRM model.
+    """Calculate the emissions using the VPRM model.
 
     This function uses timeseries of vegetation indices, temperature and radiation 
     to calculate the
@@ -88,26 +88,30 @@ def calculate_vprm_emissions(df: pd.DataFrame, df_vprm: pd.DataFrame) -> pd.Data
 
     The equations used are the following:
 
-    #. PAR (Photosynthetically Active Radiation) is calculated from the shortwave radiation:
+    PAR (Photosynthetically Active Radiation) is calculated from the shortwave radiation:
+    
     .. math::
         \\mathrm{PAR} = \\frac{\\mathrm{RAD}}{0.505}
     
-    #. Respiration is calculated from the temperature:
+    Respiration is calculated from the temperature:
+    
     .. math::
         \\mathrm{Resp} = \\alpha * T + \\beta
     
-    #. The Gross Ecosystem Exchange (GEE) is calculated from the temperature, PAR and the vegetation indices:
+    The Gross Ecosystem Exchange (GEE) is calculated from the temperature, PAR and the vegetation indices:
+    
     .. math::
-        \\mathrm{GEE} = \\lambda * T_{scale} * P_{scale} * W_{scale} * \\mathrm{EVI} * \\mathrm{PAR} * \\frac{1}{1 + \\frac{\\mathrm{PAR}}{PAR0}}
+        \\mathrm{GEE} = \\lambda * T_{scale} * P_{scale} * W_{scale} * \\mathrm{EVI} * \\frac{ \\mathrm{PAR} }{1 + \\frac{\\mathrm{PAR}}{PAR0}}
 
     where the different scales are:
 
     - :math:`T_{scale}`: Temperature scale
-    if :math:`T < T_{min}` then :math:`T_{scale} = 0`
-    else:
+
     .. math::
-        T_{scale} = \\frac{(T - T_{min}) * (T - T_{max})}{(T - T_{min}) * (T - T_{max}) + (T - T_{opt})^2}
-    
+        T_{\\text{scale}} = \\frac{(T - T_{\\text{min}}) \\cdot (T - T_{\\text{max}})}{(T - T_{\\text{min}}) \\cdot (T - T_{\\text{max}}) + (T - T_{\\text{opt}})^2} \\text{if } T \\geq T_{\\text{min}} \\text{ else } 0
+        
+
+
     - :math:`P_{scale}`: Photosynthesis scale
     .. math::
         P_{scale} = \\frac{1 + \\mathrm{LSWI}}{2}
@@ -116,52 +120,65 @@ def calculate_vprm_emissions(df: pd.DataFrame, df_vprm: pd.DataFrame) -> pd.Data
     .. math::
         W_{scale} = \\frac{1 + \\mathrm{LSWI}}{1 + \\mathrm{LSWI}_{max}}
 
-    #. The Net Ecosystem Exchange (NEE) is calculated from the respiration and GEE:
+    The Net Ecosystem Exchange (NEE) is calculated from the respiration and GEE.
+    
     .. math::
         \\mathrm{NEE} = \\mathrm{Resp} + \\mathrm{GEE}
 
+        
+    Units for all fluxes (NEE, GEE, Resp, ...) are
+
+    .. math:: 
+        \\frac{\\mu mol_{\\mathrm{CO2}}}{m^2 * s}
 
         
     Urban modifications
 
-    The VPRM model can be extended to urban areas.
+    The VPRM model can be extended to urban areas according to  [Urban_VPRM]_ .
     
-    #. The urban temperature is used instead of the global temperature.
-    #. The formula for :math:`P_{scale}` is modified to:
+    - A "urban temperature" is used instead of the global temperature to represent
+        the urban heat island phenomenon.
+    - The formula for :math:`P_{scale}` is modified to
     .. math::
         P_{scale} = \\frac{\\mathrm{EVI} - \\mathrm{EVI}_{min}}{\\mathrm{EVI}_{max} - \\mathrm{EVI}_{min}}
-    #. The respiration is calculated differently:
+    - The respiration is calculated differently
     .. math::
-        \\mathrm{Resp} = \\frac{\\mathrm{resp_e_init}}{2} * (1 - \\mathrm{ISA}) + \\frac{\\mathrm{EVI} + \\mathrm{EVI}_{min} * \\mathrm{ISA}}{\\mathrm{EVI}_{ref}} * \\frac{\\mathrm{resp_e_init}}{2}
+        \\mathrm{Resp} = \\frac{\\mathrm{Resp_{e-init}}}{2} * (1 - \\mathrm{ISA}) + \\frac{\\mathrm{EVI} + \\mathrm{EVI}_{min} * \\mathrm{ISA}}{\\mathrm{EVI}_{ref}} * \\frac{\\mathrm{Resp_{e-init}}}{2}
+    
+    where :math:`\\mathrm{Resp_{e-init}}` is the basic vprm respiration and :math:`\\mathrm{ISA}` is the impervious surface area at the vegetation location.
+    
     .. warning::
         The urban VPRM model is currently not fully implemented.
     
     
     :param df: Dataframe with the observations. It must be a multiindex dataframe with the following columns:
-        - RAD: Shortwave radiation in W/m2
+    
+        - `RAD`: Shortwave radiation in W/m2
         - ('T', 'global'): Temperature in degC
         - (vegetation_type, 'lswi'): Land Surface Water Index 
         - (vegetation_type, 'evi'): Enhanced Vegetation Index
         - ('T', 'urban'): Optional for urban VPRM. Temperature in degC (urban area)
+    
     :param df_vprm: Dataframe with the VPRM parameters.
         Each row must correspond to a vegetation type and have the following columns:
-        - alpha: Respiration parameter
-        - beta: Respiration parameter
-        - lambda: Photosynthesis parameter
-        - Tmin: Minimum temperature for the vegetation type
-        - Topt: Optimal temperature for the vegetation type
-        - Tmax: Maximum temperature for the vegetation type
-        - Tlow: Low temperature for the vegetation type
-        - PAR0: Photosynthetically Active Radiation parameter
-        - is_urban: Boolean indicating if the vegetation type is urban (optional, default is False)
-    :return: Dataframe with the emissions. Columns added are:
-        Units are:math:`\frac{\mu mol_{\mathrm{CO2}}}{m^2 * s}`
+        
+        - `alpha`: Respiration parameter
+        - `beta`: Respiration parameter
+        - `lambda`: Photosynthesis parameter
+        - `Tmin`: Minimum temperature for photosynthesis
+        - `Topt`: Optimal temperature for photosynthesis
+        - `Tmax`: Maximum temperature for photosynthesis
+        - `Tlow`: Low temperature for photosynthesis
+        - `PAR0`: Photosynthetically Active Radiation parameter
+        - `is_urban`: Boolean indicating if the vegetation type is urban (optional, default is False)
+   
+    :return: Dataframe with the emissions. Some columns are added
+
         - (vegetation_type, 'resp_min'): Respiration at the minimum temperature
         - (vegetation_type, 'resp_max'): Respiration at the maximum temperature
         - (vegetation_type, 'resp'): Respiration
         - (vegetation_type, 'gee'): Gross Ecosystem Exchange
         - (vegetation_type, 'nee'): Net Ecosystem Exchange (nee = gee - resp)
-
     """
     logger = logging.getLogger(__name__)
     df = df.copy()

emiproc/utilities.py

@@ -520,6 +520,44 @@ def confirm_prompt(question: str) -> bool:
     return reply == "y"
 
 
+def total_emissions_almost_equal(
+    total_dict_1: dict[str, dict[str, float]],
+    total_dict_2: dict[str, dict[str, float]],
+    relative_tolerance: float = 1e-5,
+) -> bool:
+    """Test that the total emissions of two inventories are almost equal.
+    
+    :arg total_dict_1: The first total emissions dictionary
+    :arg total_dict_2: The second total emissions dictionary
+    :arg relative_tolerance: The relative tolerance to use for the comparison.
+        The comparison is done as follows:
+        abs(total_dict_1 - total_dict_2) < relative_tolerance * (total_dict_1 + total_dict_2) / 2
+        
+    :returns: True if the total emissions are almost equal, False otherwise
+    :raises ValueError: If the two dictionaries have different keys.
+    """
+    for sub in total_dict_1.keys() | total_dict_2.keys():
+        if sub not in total_dict_1 or sub not in total_dict_2:
+            raise ValueError(
+                f"Subcategory {sub} is missing in one of the dictionaries"
+            )
+        for cat in total_dict_1[sub].keys() | total_dict_2[sub].keys():
+            if cat not in total_dict_1[sub] or cat not in total_dict_2[sub]:
+                raise ValueError(
+                    f"Category {cat} is missing in one of the dictionaries for substance {sub}"
+                )
+            # Get a very small proportion of the total emissions
+            err_allowed = (
+                relative_tolerance
+                * (total_dict_1[sub][cat] + total_dict_2[sub][cat])
+                / 2
+            )
+
+            if abs(total_dict_1[sub][cat] - total_dict_2[sub][cat]) > err_allowed:
+                return False
+    return True
+
+
 class ProgressIndicator:
     """Used to show progress for long operations.