Calculate spatial flux#

This example shows how to calculate a spatial greenhouse gas (GHG) flux based on modelled input from for instance SOMERS. Lusos calculates a spatial flux by taking a weighed average of the median fluxes of polygons that fall within each raster cell of a 2D grid. Following the previously explained concept, lusos computes the area of each polygon within a raster cell. The computed areas are used as the weights to assign the computed flux to each cell. This example shows the simple steps how to calculate to spatial flux. We will use the same example area from the coverage example.

import geopandas as gpd
from shapely import geometry as gmt

import lusos

# Bounding box for the example area
xmin, ymin, xmax, ymax = 111_000, 455_000, 116_000, 460_000
study_area = gpd.GeoDataFrame(geometry=[gmt.box(xmin, ymin, xmax, ymax)], crs=28992)

emissions = lusos.data.sample_emissions()
Downloading file 'emissions.geoparquet' from 'https://github.com/Deltares-research/lulucf-somers/raw/main/data/emissions.geoparquet' to '/home/runner/.cache/lusos'.

Let’s first checkout the emissions data:

print(emissions.head())
m = emissions.explore() # Display the emissions data on a map
study_area.explore(m=m, color="black", fill=False)
  PARCEL_ID      AREA    emission_t  \
0  100358_1  2.934975  15735.972375   
1   47019_1  1.460252  10960.169343   
2  440920_1  3.821042  17659.050978   
3  357317_1  1.992677  10410.785529   
4  442000_1  0.880378  12365.247052   

                                            geometry  
0  POLYGON Z ((115148.692 459977.124 0, 115160.20...  
1  POLYGON Z ((115534.479 459528.986 0, 115532.26...  
2  POLYGON Z ((115518.22 460249.249 0, 115516.605...  
3  POLYGON Z ((116024.58 459509.059 0, 116075.902...  
4  POLYGON Z ((115153.688 459076.022 0, 115150.11...  
Make this Notebook Trusted to load map: File -> Trust Notebook

The data contains emission values for polygons in ton/ha. Lusos expects this unit for the input data and translates this into a flux per m2. Using a different therefore means that the output will be in the wrong unit as well and must be translated to ton/m2.

Lusos needs two columns to be present in the input data: “parcel_id” and “median”. This is because lusos bases the calculation on the median emission and we need an id for each polygon to be able to identify it and correctly assign fluxes to cells. Without the presence of these columns, an error will be thrown. We can rename the column “emission_t” because this will be used as the median and make the other columns lowercase strings because “PARCEL_ID” is present in the data.

emissions.columns = emissions.columns.str.lower()
emissions.rename(columns={"emission_t": "median"}, inplace=True)
emissions.head()
parcel_id area median geometry
0 100358_1 2.934975 15735.972375 POLYGON Z ((115148.692 459977.124 0, 115160.20...
1 47019_1 1.460252 10960.169343 POLYGON Z ((115534.479 459528.986 0, 115532.26...
2 440920_1 3.821042 17659.050978 POLYGON Z ((115518.22 460249.249 0, 115516.605...
3 357317_1 1.992677 10410.785529 POLYGON Z ((116024.58 459509.059 0, 116075.902...
4 442000_1 0.880378 12365.247052 POLYGON Z ((115153.688 459076.022 0, 115150.11...

Now we can define our grid again for the calculation and calculate the flux.

grid = lusos.LassoGrid(xmin, ymin, xmax, ymax, 25, 25)
flux_per_m2 = lusos.calculate_somers_emissions(emissions, grid)
flux_per_m2
/home/runner/work/lulucf-somers/lulucf-somers/.pixi/envs/default/lib/python3.13/site-packages/xugrid/regrid/structured.py:606: FutureWarning: In a future version of xarray the default value for compat will change from compat='no_conflicts' to compat='override'. This is likely to lead to different results when combining overlapping variables with the same name. To opt in to new defaults and get rid of these warnings now use `set_options(use_new_combine_kwarg_defaults=True) or set compat explicitly.
  ds = xr.merge([ds_x, ds_y])
<xarray.DataArray (y: 200, x: 200)> Size: 320kB
array([[0.85330532, 0.86892862, 0.87292664, ..., 0.67917049, 0.74290027,
        0.74290027],
       [0.85824311, 0.87292664, 0.89821724, ..., 0.37916474, 0.592546  ,
        0.74290027],
       [0.86999863, 0.87292664, 0.97313272, ..., 0.39592677, 0.37290132,
        0.49532248],
       ...,
       [       nan,        nan,        nan, ...,        nan,        nan,
               nan],
       [       nan,        nan,        nan, ...,        nan,        nan,
               nan],
       [       nan,        nan,        nan, ...,        nan,        nan,
               nan]], shape=(200, 200))
Coordinates:
  * y            (y) float64 2kB 4.6e+05 4.6e+05 4.599e+05 ... 4.55e+05 4.55e+05
  * x            (x) float64 2kB 1.11e+05 1.11e+05 ... 1.16e+05 1.16e+05
    spatial_ref  int64 8B 0