diff --git a/emiproc/grids.py b/emiproc/grids.py index 4c02a34..7906749 100644 --- a/emiproc/grids.py +++ b/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). diff --git a/emiproc/profiles/vprm.py b/emiproc/profiles/vprm.py index 65176fa..cc03c5e 100644 --- a/emiproc/profiles/vprm.py +++ b/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() diff --git a/emiproc/utilities.py b/emiproc/utilities.py index 97bba55..2d0cc1e 100644 --- a/emiproc/utilities.py +++ b/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. diff --git a/examples/scale_inventory_with_footprints.ipynb b/examples/scale_inventory_with_footprints.ipynb new file mode 100644 index 0000000..0d7ef00 --- /dev/null +++ b/examples/scale_inventory_with_footprints.ipynb @@ -0,0 +1,217 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Scaling inventory on footprints \n", + "\n", + "This example shows how we can retrieve how much emissions are being seen at \n", + "a measurement location, based on geospatial footprint of that site." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from pathlib import Path\n", + "import numpy as np\n", + "import xarray as xr\n", + "import rioxarray as rxr\n", + "import geopandas as gpd\n", + "from shapely.creation import polygons\n", + "from emiproc.grids import RegularGrid, GeoPandasGrid\n", + "from emiproc.regrid import remap_inventory\n", + "from emiproc.inventories.zurich import MapLuftZurich\n", + "from emiproc.inventories.utils import group_categories\n", + "\n", + "\n", + "\n", + "footprints_dir = Path('footprints')\n", + "\n", + "footprint_file = footprints_dir / 'zurich_footprint_raster_fclim2d_2303071400.txt'\n", + "x_file = footprints_dir / 'zurich_footprint_raster_x2d2303071400.txt'\n", + "y_file = footprints_dir / 'zurich_footprint_raster_y2d2303071400.txt'\n", + "\n", + "measurement_coordinates = 2680911.322, 1248390.798" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "\n", + "footprint = np.loadtxt(footprint_file)\n", + "x = np.loadtxt(x_file)\n", + "y = np.loadtxt(y_file)\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Get only the cells of the grid where footprint is greater than 0\n", + "mask = footprint > 1e-20\n", + "x_coords = x[mask] + measurement_coordinates[0]\n", + "y_coords = y[mask] + measurement_coordinates[1]\n", + "d2 = 5.\n", + "coords = np.array(\n", + " [\n", + " [x, y]\n", + " for x, y in zip(\n", + " [x_coords -d2, x_coords - d2, x_coords + d2, x_coords + d2],\n", + " [y_coords - d2, y_coords + d2, y_coords + d2, y_coords - d2],\n", + " )\n", + " ]\n", + ")\n", + "coords = np.rollaxis(coords, -1, 0)\n", + "# Create a gdf with the mask\n", + "gdf = gpd.GeoDataFrame(\n", + " {\n", + " 'footprint': footprint[mask],\n", + " },\n", + " geometry=polygons(coords),\n", + " crs=\"LV95\"\n", + ")\n", + "grid = GeoPandasGrid(gdf)\n", + "gdf" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "file = Path(\"mapLuft_2020_v2021.gdb\")\n", + "\n", + "zh_inv = MapLuftZurich(file, ['CO2'])\n", + "\n", + "# Group the categories \n", + "ZH_2_GNFR_CO2 = {\n", + " # PublicPower\n", + " \"GNFR_A\": [\n", + " \"c2201_BHKW_Emissionen_Kanton\", # ???\n", + " \"c2301_KHKWKehricht_Emissionen_Kanton\",\n", + " \"c2302_KHKWErdgas_Emissionen_Kanton\",\n", + " \"c2303_KHKWHeizoel_Emissionen_Kanton\",\n", + " ],\n", + " # Industry\n", + " \"GNFR_B\": [\n", + " \"c3201_Notstromanlagen_Emissionen_Kanton\",\n", + " \"c3301_Prozessenergie_Emissionen_Kanton\",\n", + " \"c3410_Bierbrauereien_Emissionen_Kanton\",\n", + " \"c3411_Brotproduktion_Emissionen_Kanton\",\n", + " ],\n", + " # Other stationary combustion (services, residential, agriculture)\n", + " \"GNFR_C\": [\n", + " \"c2101_Oelheizungen_Emissionen_Kanton\",\n", + " \"c2102_Gasheizungen_Emissionen_Kanton\",\n", + " \"c2103_HolzheizungenLokalisiert_Emissionen_Kanton\",\n", + " \"c2104_HolzheizungenDispers_Emissionen_Kanton\",\n", + " \"c2105_Warmwassererzeuger_Emissionen_Kanton\",\n", + " ],\n", + " # Fugitives\n", + " \"GNFR_D\": [\n", + " \"c3416_Tankstellen_Emissionen_Kanton\",\n", + " ],\n", + " # Solvents and product use\n", + " #\"GNFR_E\": [\n", + " #],\n", + " # Road transport\n", + " \"GNFR_F\": [\n", + " \"c1301_Personenwagen_Emissionen_Kanton\",\n", + " \"c1302_Lastwagen_Emissionen_Kanton\",\n", + " \"c1303_Motorraeder_Emissionen_Kanton\",\n", + " \"c1304_Linienbusse_Emissionen_Kanton\",\n", + " \"c1306_StartStopTankatmung_Emissionen_Kanton\",\n", + " \"c1307_Lieferwagen_Emissionen_Kanton\",\n", + " \"c1308_Reisebusse_Emissionen_Kanton\",\n", + " ],\n", + " # Shipping\n", + " \"GNFR_G\": [\n", + " \"c1101_Linienschiffe_Emissionen_Kanton\",\n", + " \"c1102_PrivaterBootsverkehr_Emissionen_Kanton\",\n", + " ],\n", + " # Offroad mobility\n", + " \"GNFR_I\": [\n", + " # c31xx are construction stuff\n", + " \"c3101_MaschinenHochbau_Emissionen_Kanton\",\n", + " \"c3104_MaschinenTiefbau_Emissionen_Kanton\",\n", + " \"c3419_IndustrielleFZ_Emissionen_Kanton\",\n", + " \"c4101_ForstwirtschaftlicheFZ_Emissionen_Kanton\",\n", + " \"c4201_LandwirtschaftlicheFZ_Emissionen_Kanton\",\n", + " ],\n", + " # Waste\n", + " \"GNFR_J\": [\n", + " \"c2401_Klaerschlammverwertung_Emissionen_Kanton\",\n", + " \"c3418_Vergaerwerk_Emissionen_Kanton\",\n", + " \"c3414_Krematorium_Emissionen_Kanton\",\n", + " \"c5201_Gruenabfallverbrennung_Emissionen_Kanton\",\n", + " \"c5301_HolzoefenKleingarten_Emissionen_Kanton\",\n", + " \"c5401_AbfallverbrennungHaus_Emissionen_Kanton\",\n", + " ],\n", + " # AgriLivestock\n", + " # AgriOther\n", + " #\"GNFR_L\": [\n", + " #],\n", + " # Others\n", + " \"GNFR_R\": [\n", + " \"c5601_Feuerwerke_Emissionen_Kanton\",\n", + " \"c5701_Tabakwaren_Emissionen_Kanton\",\n", + " \"c5801_BrandFeuerschaeden_Emissionen_Kanton\",\n", + " ],\n", + "}\n", + "\n", + "zh_groupped = group_categories(zh_inv, ZH_2_GNFR_CO2)\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "regrided = remap_inventory(zh_groupped, grid)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "regrided.total_emissions.plot.bar()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.2" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/tests/test_remap_inv.py b/tests/test_remap_inv.py index 1ba9bc0..50b8fb2 100644 --- a/tests/test_remap_inv.py +++ b/tests/test_remap_inv.py @@ -5,6 +5,7 @@ from emiproc.tests_utils.test_grids import regular_grid, gpd_grid from emiproc.inventories.utils import get_total_emissions +from emiproc.utilities import total_emissions_almost_equal @@ -22,7 +23,7 @@ def test_remap(): # Check the total emissions dictionaries are the same (grid is larger than inventory) total_inv = get_total_emissions(inv_with_pnt_sources) total_remapped = get_total_emissions(remaped_inv) - assert total_inv == total_remapped + assert total_emissions_almost_equal(total_inv, total_remapped) def test_remap_keep_shapes(): @@ -48,7 +49,7 @@ def test_remap_keep_shapes(): # Check the total emissions dictionaries are the same (grid is larger than inventory) total_inv = get_total_emissions(inv_with_pnt_sources) total_remapped = get_total_emissions(remaped_inv) - assert total_inv == total_remapped + assert total_emissions_almost_equal(total_inv, total_remapped) diff --git a/tests/test_utilities.py b/tests/test_utilities.py new file mode 100644 index 0000000..a58f5ea --- /dev/null +++ b/tests/test_utilities.py @@ -0,0 +1,40 @@ +import pytest +from emiproc.utilities import total_emissions_almost_equal + +def test_total_emissions_almost_equal(): + # Create a reference dictionary + ref_dict = { + 'sub1': {'cat1': 10.0, 'cat2': 20.0}, + 'sub2': {'cat1': 30.0, 'cat2': 40.0} + } + + # Test case 1: Total emissions are exactly equal + assert total_emissions_almost_equal(ref_dict, ref_dict) == True + + # Test case 2: Total emissions are almost equal within the default tolerance + total_dict_2 = { + 'sub1': {'cat1': 10.00001, 'cat2': 19.99999}, + 'sub2': {'cat1': 29.99999, 'cat2': 40.00001} + } + assert total_emissions_almost_equal(ref_dict, total_dict_2) == True + + # Test case 3: Total emissions are not almost equal within the default tolerance + total_dict_3 = { + 'sub1': {'cat1': 10.0001, 'cat2': 19.99}, + 'sub2': {'cat1': 29.9999, 'cat2': 40.01} + } + assert (total_emissions_almost_equal(ref_dict, total_dict_3) == False) + + # Test case 4: Total emissions have different subcategories + total_dict_4 = { + 'sub1': {'cat1': 10.0, 'cat2': 20.0}, + 'sub3': {'cat1': 30.0, 'cat2': 40.0} + } + pytest.raises(ValueError, total_emissions_almost_equal, ref_dict, total_dict_4) + + # Test case 5: Total emissions have different categories + total_dict_5 = { + 'sub1': {'cat1': 10.0, 'cat3': 20.0}, + 'sub2': {'cat1': 30.0, 'cat2': 40.0} + } + pytest.raises(ValueError, total_emissions_almost_equal, ref_dict, total_dict_5)