Update livestock sector implementation to use headcounts dataset

CLAUDE.md

@@ -64,7 +64,16 @@ Copy `.env.example` to `.env`. Key variables:
 
 ## Testing
 
-Tests live in `tests/unit/` and `tests/integration/`. Integration tests exercise the full pipeline with real (cached) data. Use `--cache` pytest fixtures to control data caching behaviour.
+Tests live in `tests/unit/` and `tests/integration/`. Tests follow pytest
+conventions and use pytest libraries.
+
+Integration tests exercise the full pipeline with real (cached) data. Use
+`input_files` pytest fixture to use cached data. For integration tests which
+should download DataSource assets from the internet, `input_files` should be
+omitted.
+
+Unit tests should avoid calling methods which read from the filesystem, and
+should construct example data to pass to methods.
 
 ## Code Style
 

changelog/183.feature.md

@@ -0,0 +1 @@
+Update livestock sector implementation to use headcounts dataset

docs/data-sources.md

@@ -140,12 +140,9 @@ Superpower Institute, with data synced using
 Some sectoral inventories are spatialised by identifying which Australian Land
 Use and Management (ALUM) Classification codes fall within those sectors, and
 then distributing inventory emissions to the areas where those ALUM codes are
-featured in the Land Use of Australia geographical dataset provided by the
-Australian Department of Agriculture.
-
-This is done using a combination of the official
+featured in the
 [Land Use of Australia](https://www.agriculture.gov.au/abares/aclump/land-use/land-use-of-australia-2010-11-to-2020-21)
-GeoTIFF raster, and a manual [mapping between ALUM codes and Open Methane sectors][4].
+geographical dataset provided by the Australian Department of Agriculture.
 
 Due to issues fetching the file directly from the agriculture.gov.au service
 in cloud-hosted services, the land use dataset is mirrored from the Open
@@ -200,33 +197,25 @@ substances which can act as a proxy for methane emission.
 
 ## Sector: Livestock
 
-Australian national dataset of CH4 flux estimates from enteric fermentation in livestock.
+Australian national dataset of livestock headcounts for cattle and sheep, used
+to distribute ANGA inventory emissions for enteric fermentation and manure for
+those animals.
 
 ### Sources
 
-Enteric fermentation emissions generated by CSIRO Ag. and Food using livestock
-census data and UNFCCC emissions factors. Underlying livestock numbers
-taken from: "Navarro, J. Marcos Martinez, R. (2021) Estimating long-term profits, fertiliser and pesticide use
-baselines in Australian agricultural regions. User Guide. CSIRO,
-Australia." 
+Livestock headcounts generated by CSIRO Agriculture and Food using livestock
+census data. Attribution TBC 
 
-- Dataset: [EntericFermentation.nc][2]
+- Dataset: TBC
 - Resolution: 0.01 degree
-- Period: 2020, single annual average
+- Period: Unknown
 - Updates: never
 
 ### Considerations
 
-- Spatial distribution might change during the year since farmers are moving their cows around
-- Doesn't appear to be available online
-
-### Alternative candidates
-
-- https://publications.csiro.au/publications/publication/PIcsiro:EP2022-1389/
-  - Resolution: 0.1 degree
-  - Period: 1972 - 2019, annual averages
-  - Global
-- Electronic ear tags might be a useful database
+- Real spatial distribution likely to be highly variable during the year as
+  grazing herds move through available pastures.
+- Dataset not yet publicly available.
 
 
 ## Sector: Termites
@@ -238,7 +227,7 @@ Global dataset of CH4 flux estimates from termites.
 Termite emissions used in [Saunois et al. 2020](https://essd.copernicus.org/articles/12/1561/2020/)
 supplied by Simona Castaldi and Sergio Noce.
 
-- Dataset: [termite_emissions_2010-2016.nc][5]
+- Dataset: [termite_emissions_2010-2016.nc][4]
 - Resolution: 0.5 degree
 - Period: mean of 2010 – 2016
 - Updates: never
@@ -305,7 +294,7 @@ and described in Zhang et al. (2023 under review) (could this be https://essd.co
 
 ## Sector: Agriculture
 
-Agricultural emissions (excluding livestock) reported in the
+Agricultural emissions (excluding cattle and sheep) reported in the
 [Australian UNFCCC Inventory](#Australian-UNFCCC-Inventory)
 are spatialised according to the [Land Use of Australia](#Land-Use-of-Australia) dataset.
 
@@ -475,5 +464,4 @@ spatialised according to the [Land Use of Australia](#Land-Use-of-Australia) dat
 [1]: https://openmethane.s3.amazonaws.com/prior/inputs/DLEM_totflux_CRU_diagnostic.nc
 [2]: https://openmethane.s3.amazonaws.com/prior/inputs/EntericFermentation.nc
 [3]: https://openmethane.s3.amazonaws.com/prior/inputs/ch4-electricity.csv
-[4]: https://openmethane.s3.amazonaws.com/prior/inputs/landuse-sector-map.csv
-[5]: https://openmethane.s3.amazonaws.com/prior/inputs/termite_emissions_2010-2016.nc
+[4]: https://openmethane.s3.amazonaws.com/prior/inputs/termite_emissions_2010-2016.nc

src/openmethane_prior/data_sources/landuse/__init__.py

@@ -1,2 +1 @@
-from .data import alum_sector_mapping_data_source, landuse_map_data_source
-from .alum import alum_codes_for_sector
+from .data import landuse_map_data_source

src/openmethane_prior/data_sources/landuse/data.py

@@ -23,11 +23,6 @@
 from openmethane_prior.lib import ConfiguredDataSource, DataSource
 from openmethane_prior.lib.data_manager.parsers import parse_csv
 
-alum_sector_mapping_data_source = DataSource(
-    name="alum-sector-mapping",
-    url="https://openmethane.s3.amazonaws.com/prior/inputs/landuse-sector-map.csv",
-    parse=parse_csv,
-)
 
 # filename of the GeoTIFF file inside the official zip file source
 NLUM_GEOTIFF_FILENAME="NLUM_v7_250_ALUMV8_2020_21_alb.tif"

src/openmethane_prior/lib/__init__.py

@@ -23,6 +23,7 @@
     DataManager,
     DataSource,
 )
+from .grid.grid import Grid
 from .grid.domain import Domain
 from .grid.regrid import regrid_data
 from .regrid import regrid_data_array_conservative

src/openmethane_prior/sectors/agriculture/sector.py

@@ -26,11 +26,7 @@
     inventory_data_source,
     inventory_domain_data_source,
 )
-from openmethane_prior.data_sources.landuse import (
-    alum_codes_for_sector,
-    alum_sector_mapping_data_source,
-    landuse_map_data_source,
-)
+from openmethane_prior.data_sources.landuse import landuse_map_data_source
 from openmethane_prior.lib import (
     kg_to_period_cell_flux,
     logger,
@@ -42,6 +38,17 @@
 
 logger = logger.get_logger(__name__)
 
+# ALUM Classification Version 8
+# Source: https://www.agriculture.gov.au/abares/aclump/land-use/alum-classification
+alum_codes_agriculture = [
+    520, # Intensive animal production
+    523, # Poultry farms
+    524, # Piggeries
+    525, # Aquaculture
+    526, # Horse studs
+    527, # Saleyards/stockyards
+    542, # Rural residential with agriculture
+]
 
 def process_emissions(
         sector: PriorSector,
@@ -50,11 +57,6 @@ def process_emissions(
 ):
     config = sector_config.prior_config
 
-    # Import a map of land use type numbers to emissions sectors
-    # make a dictionary of all landuse types corresponding to sectors in map
-    sector_mapping_asset = sector_config.data_manager.get_asset(alum_sector_mapping_data_source)
-    sector_alum_codes = alum_codes_for_sector(sector.name, sector_mapping_asset.data)
-
     # load the national inventory data, ready to calculate sectoral totals
     emissions_inventory = sector_config.data_manager.get_asset(inventory_data_source).data
 
@@ -80,7 +82,7 @@ def process_emissions(
     )
 
     # create a mask of pixels which match the sector code
-    sector_mask = np.isin(dataBand, sector_alum_codes)
+    sector_mask = np.isin(dataBand, alum_codes_agriculture)
     sector_xr = xr.DataArray(sector_mask, coords={ 'y': lu_y, 'x': lu_x  })
 
     # now aggregate to coarser resolution of the domain grid
@@ -111,8 +113,11 @@ def process_emissions(
 sector = PriorSector(
     name="agriculture",
     emission_category="anthropogenic",
-    unfccc_categories=[ # All Agriculture, except Enteric Fermentation
-        "3.B", # Manure Management
+    unfccc_categories=[ # All Agriculture, except cattle and sheep
+        "3.A.3", # Enteric Fermentation - Swine
+        "3.A.4", # Enteric Fermentation - Other Livestock
+        "3.B.3", # Manure Management - Swine
+        "3.B.4", # Manure Management - Other Livestock
         "3.C", # Rice Cultivation
         "3.D", # Agricultural Soils
         "3.E", # Prescribed Burning of Savannas

src/openmethane_prior/sectors/livestock/sector.py

@@ -15,64 +15,124 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 #
-
 import numpy as np
+import pandas as pd
+import pyproj
 import xarray as xr
 
+from openmethane_prior.data_sources.inventory import (
+    inventory_data_source,
+    get_sector_emissions_by_code,
+)
+from openmethane_prior.lib.data_manager.parsers import parse_csv
 from openmethane_prior.lib import (
     DataSource,
-    convert_to_timescale,
+    Grid,
     logger,
     PriorSector,
     PriorSectorConfig,
 )
+from openmethane_prior.lib.units import seconds_in_period
 
 logger = logger.get_logger(__name__)
 
-livestock_data_source = DataSource(
-    name="enteric-fermentation",
-    url="https://openmethane.s3.amazonaws.com/prior/inputs/EntericFermentation.nc",
+livestock_headcount_data_source = DataSource(
+    name="livestock-headcount",
+    url="https://openmethane.s3.amazonaws.com/prior/inputs/lmap_selected_cols_Nasimeh_171025.csv",
+    parse=parse_csv,
 )
 
+def gridded_livestock_headcounts(
+    livestock_df: pd.DataFrame,
+    domain_grid: Grid,
+) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
+    """Given a list of coordinates for livestock areas with headcounts of beef
+    cattle, dairy cattle and sheep, calculate the annual methane emissions in
+    each area, and return as gridded emissions over the domain of interest."""
+    # find the domain grid cell (ix, iy) each livestock area centre point is in
+    livestock_df["ix"], livestock_df["iy"], valid_cells = domain_grid.lonlat_to_cell_index(
+        lon=livestock_df["lon"],
+        lat=livestock_df["lat"],
+    )
+
+    # discard data outside the domain
+    livestock_df = livestock_df[valid_cells]
+
+    logger.debug(f"Mapping livestock headcounts to domain grid")
+    beef_gridded = np.zeros(domain_grid.shape)
+    dairy_gridded = np.zeros(domain_grid.shape)
+    sheep_gridded = np.zeros(domain_grid.shape)
+
+    row_idx = (livestock_df["iy"], livestock_df["ix"])
+    np.add.at(beef_gridded, row_idx, livestock_df["heads_mapped_mix_beef"])
+    np.add.at(dairy_gridded, row_idx, livestock_df["heads_mapped_dairy"])
+    np.add.at(sheep_gridded, row_idx, livestock_df["heads_mapped_mix_sheep"])
+
+    return beef_gridded, dairy_gridded, sheep_gridded
+
+
 def process_emissions(
     sector: PriorSector,
     sector_config: PriorSectorConfig,
     prior_ds: xr.Dataset,
 ):
     config = sector_config.prior_config
+    domain_grid = config.domain().grid
 
-    livestock_asset = sector_config.data_manager.get_asset(livestock_data_source)
-    with xr.open_dataset(livestock_asset.path) as lss:
-        ls = lss.load()
+    livestock_headcount_da = sector_config.data_manager.get_asset(livestock_headcount_data_source)
+    livestock_df: pd.DataFrame = livestock_headcount_da.data
 
-    domain_grid = config.domain().grid
+    # dataset X/Y coords are ESPF:4283 (GDA94), convert to EPSG:4326
+    transformer = pyproj.Transformer.from_crs(4283,4326)
+    livestock_df["lon"], livestock_df["lat"] = transformer.transform(
+        yy=livestock_df["Y"], xx=livestock_df["X"],
+    )
+
+    # place gridded headcounts on the grid
+    beef_np, dairy_np, sheep_np = gridded_livestock_headcounts(livestock_df, domain_grid)
+
+    # get inventory totals using animal-specific UNFCCC sectors
+    emissions_inventory = sector_config.data_manager.get_asset(inventory_data_source).data
+    get_emissions_params = dict(
+        emissions_inventory=emissions_inventory,
+        start_date=config.start_date,
+        end_date=config.end_date,
+    )
+    beef_total_emissions = get_sector_emissions_by_code(
+        category_codes=["3.A.1.b", "3.A.1.c", "3.B.1.b", "3.B.1.c"], **get_emissions_params,
+    )
+    dairy_total_emissions = get_sector_emissions_by_code(
+        category_codes=["3.A.1.a", "3.B.1.a"], **get_emissions_params,
+    )
+    sheep_total_emissions = get_sector_emissions_by_code(
+        category_codes=["3.A.2", "3.B.2"], **get_emissions_params,
+    )
+    logger.debug(f"Total emissions from beef: {beef_total_emissions / 1e6:.2f} kt, dairy: {dairy_total_emissions / 1e6:.2f} kt, sheep: {sheep_total_emissions / 1e6:.2f} kt")
 
-    # Re-project into domain coordinates
-    # - create meshgrids of the lats and lons
-    lonmesh, latmesh = np.meshgrid(ls.lon, ls.lat)
-    cell_x, cell_y, cell_valid = domain_grid.lonlat_to_cell_index(lonmesh, latmesh)
+    beef_total_headcount = livestock_df["heads_mapped_mix_beef"].sum()
+    dairy_total_headcount = livestock_df["heads_mapped_dairy"].sum()
+    sheep_total_headcount = livestock_df["heads_mapped_mix_sheep"].sum()
 
-    enteric_as_array = lss.CH4_total.to_numpy()
+    # distribute the national emission totals by headcount proportion
+    beef_ch4_np = beef_total_emissions * beef_np / beef_total_headcount
+    dairy_ch4_np = dairy_total_emissions * dairy_np / dairy_total_headcount
+    sheep_ch4_np = sheep_total_emissions * sheep_np / sheep_total_headcount
 
-    livestockCH4 = np.zeros(domain_grid.shape)
-    logger.info("Distribute livestock CH4 (long process)")
-    # we're accumulating emissions from fine to coarse grid
-    # accumulate in mass units and divide by area at end
-    for j in range(ls.lat.size):
-        ix, iy = cell_x[j,:], cell_y[j,:]
-        # input domain is bigger so mask indices out of range
-        mask = cell_valid[j, :]
-        if mask.any():
-            # the following needs to use .at method since iy,ix indices may be repeated and we need to acumulate
-            np.add.at(livestockCH4, (iy[mask], ix[mask]), enteric_as_array[j, mask])
+    total_ch4_np = beef_ch4_np + dairy_ch4_np + sheep_ch4_np
 
-    return convert_to_timescale(livestockCH4, domain_grid.cell_area)
+    # convert CH4 kg to CH4 kg/s/m2
+    return total_ch4_np / domain_grid.cell_area / seconds_in_period(config.start_date, config.end_date)
 
 
 sector = PriorSector(
     name="livestock",
     emission_category="anthropogenic",
-    unfccc_categories=["3.A"], # Enteric Fermentation
+    unfccc_categories=[
+        "3.A.1", # Enteric Fermentation - Cattle
+        "3.A.2", # Enteric Fermentation - Sheep
+        "3.B.1", # Manure Management - Cattle
+        "3.B.2", # Manure Management - Sheep
+    ],
     cf_standard_name="domesticated_livestock",
     create_estimate=process_emissions,
 )