Merge pull request #73 from C2SM-RCM/coli

unit conversion for wrf and some other day capabilities

emiproc/exports/wrf.py

@@ -14,11 +14,13 @@
 from shapely.creation import polygons
 
 from emiproc.exports.utils import get_temporally_scaled_array
-from emiproc.grids import WGS84, Grid
+from emiproc.grids import WGS84, Grid, RegularGrid
 from emiproc.inventories import Inventory
+from emiproc.utilities import HOUR_PER_DAY, get_day_per_year
+from emiproc.utils.constants import get_molar_mass
 
 
-class WRF_Grid(Grid):
+class WRF_Grid(RegularGrid):
     """Grid of the wrf model.
 
     The grid is a pseudo regular grid, in the sense that the grid is regular
@@ -41,7 +43,7 @@ def __init__(self, grid_filepath: PathLike):
         """
 
         grid_filepath = Path(grid_filepath)
-        super().__init__(name=grid_filepath.stem, crs=WGS84)
+        Grid.__init__(self, name=grid_filepath.stem, crs=WGS84)
 
         ds = xr.open_dataset(grid_filepath, engine="netcdf4")
         self.attributes = ds.attrs
@@ -127,6 +129,8 @@ def export_wrf_hourly_emissions(
 ) -> Path:
     """Export the inventory to WRF chemi files.
 
+    Output units are in mole/km2/hour.
+
     .. note::
         When running this function on Windows, the files will be saved with the
         `:` replaced by `-` in the file name, because Windows does not allow `:` in
@@ -157,6 +161,18 @@ def export_wrf_hourly_emissions(
 
     da = get_temporally_scaled_array(inv, time_range, sum_over_cells=False)
 
+    # Molar mass conversion mol / kg
+    mm_factor = 1 / (
+        xr.DataArray([get_molar_mass(sub) for sub in inv.substances], dims="substance")
+        * 1e-3
+    )
+    # year / hour
+    temporal_conversion = 1 / float(get_day_per_year(inv.year) * HOUR_PER_DAY)
+    # km2 / cell (km2/m2 * m2/cell)
+    spatial_conversion = xr.DataArray(1e-6 / grid.cell_areas, dims="cell")
+
+    da = da * mm_factor * temporal_conversion * spatial_conversion
+
     # Unstack the datarray to get on the regular 2D grid
     shape = grid.shape
     x_index = np.arange(shape[0])
@@ -199,6 +215,7 @@ def export_wrf_hourly_emissions(
                 {
                     "emiproc": f"This file was created by emiproc on {datetime.now()}",
                     "emiproc_history": f"Created from the inventory {inv.name} with {inv.history=}",
+                    "unit": "moles/km2/h",
                 }
             )
         )

emiproc/inventories/__init__.py

@@ -256,6 +256,7 @@ def copy(self, no_gdfs: bool = False, profiles: bool = True) -> Inventory:
         inv = Inventory()
         inv.__class__ = self.__class__
         inv.history = deepcopy(self.history)
+        inv.year = self.year
         if hasattr(self, "grid"):
             inv.grid = self.grid
 

emiproc/utilities.py

@@ -34,6 +34,17 @@
 HOUR_PER_YR = DAY_PER_YR * HOUR_PER_DAY
 
 
+def get_day_per_year(year: int | None) -> int | float:
+    """Get the number of days in a year accounting for leap years."""
+    if year is None:
+        logger = logging.getLogger("emiproc.get_day_per_year")
+        logger.warning("Year is None, using 365.25")
+        return DAY_PER_YR
+    if year % 4 == 0 and (year % 100 != 0 or year % 400 == 0):
+        return 366
+    return 365
+
+
 class Units(Enum):
     """Units for emissions."""
 

emiproc/utils/constants.py

@@ -0,0 +1,18 @@
+# Molar mass in g / mol
+MOLAR_MASSES_ = {
+    "CH4": 16.04,
+    "CO2": 44.01,
+    # This is a test value
+    "test": 1.0,
+    "test2": 2.0,
+}
+
+
+def get_molar_mass(substance: str) -> float:
+    """Get the molar mass of a substance in g/mol."""
+    if substance not in MOLAR_MASSES_:
+        raise ValueError(
+            f"Unknown molar mass for substance `{substance}`."
+            f"Please add it to the MOLAR_MASSES_ dictionary in {__file__}."
+        )
+    return MOLAR_MASSES_[substance]

tests/utils/test_day_per_year.py

@@ -0,0 +1,26 @@
+from emiproc.utilities import get_day_per_year
+
+
+def test_get_day_per_year():
+    days = get_day_per_year(2000)
+    assert days == 366
+
+
+def test_get_day_per_year_2():
+    days = get_day_per_year(2001)
+    assert days == 365
+
+
+def test_get_day_per_year_none():
+    days = get_day_per_year(None)
+    assert days == 365.25
+
+
+def test_get_day_per_year_century():
+    days_2100 = get_day_per_year(2100)
+    days_2000 = get_day_per_year(2000)
+    assert days_2100 == 365
+    assert days_2000 == 366
+
+    days_2400 = get_day_per_year(2400)
+    assert days_2400 == 366

tests/utils/test_molar_mass.py

@@ -0,0 +1,15 @@
+import pytest
+from emiproc.utils.constants import get_molar_mass
+
+
+def test_molar_mass():
+
+    mm_ch4 = get_molar_mass("CH4")
+
+    assert mm_ch4 == 16.04
+
+
+def test_no_molar_mass_raises():
+
+    with pytest.raises(ValueError):
+        get_molar_mass("SOMETHING UNKNOWN")