feat: ingest gpw livestock dataset by jeffseif · Pull Request #3 · cornerstone-data/luc

jeffseif · 2026-09-23T23:17:36Z

Background

We now ingest GPW_LIVESTOCK which is a ~100m raster for the head-per-area density of buffalo, cattle, goat, horse, and sheep globally.

This is step one towards being able to attribute LUC emissions to livestock commodities, like beef. Step two will be extending the FAOSTAT production dataset to ingest livestock production (today it only covers crops). And step three will be bringing these into the statistical attribution model.

Unrelated commits

I noticed that the literal we have been using for CC BY 4.0 has been inconsistent -- so I adopted the version on their website.
I noticed that our specification of EPSG's was inconsistent, so I adopted int throughout.

Notes

We reproject to EPSG:4326 during ingestion so that harmonize can remain a simple resample (not a reprojection)
The source is provided in ESRI:54052 at ~1km resolution, so in order for this data to land in a non-lossy way in EPSG:4326 we had to add a new grid whose resolution is no coarser (at ~900m). We also took the opportunity to use a grid which tiles into MAPSPAM 10:1.
Because averages applied during resampling ignore nan's by default, we need to map them to 0, filling the entire destination pixels with data. Otherwise the operation would be biased high within pixels which straddle data/nodata like the coast.
We convert from head/km² to head/ha to align with the other intensive datasets we ingest.
We only ingest 2000, 2005, 2010, and 2020 because they match those of MAPSPAM -- and thus are the only epochs useable in our statistical allocation. The dataset covers through 2022, so if we wanted those (e.g., for JD or D attribution at other assessment years), we could extend that, bump the dataset version, re-ingest (~minutes), and then re-run downstreams.
The source dataset had three pixels with large negative headcount values in CATTLE/2020. Presumably this was an error aggregating over the source's nodata values of -32_000:int16. As a result, we clip the data with min=0.

Three non-physical pixels

This is fully ingested over the 280 GNW tiles.

gcloud storage ls gs://cornerstone-ingest-us-central1/raster/gpw/livestock/v0/ten-degree-tile/ | wc -l
     280

gdalinfo

gdalinfo /vsigs/cornerstone-ingest-us-central1/raster/gpw/livestock/v0/ten-degree-tile/10S_050W.tif
Driver: GTiff/GeoTIFF
Files: /vsigs/cornerstone-ingest-us-central1/raster/gpw/livestock/v0/ten-degree-tile/10S_050W.tif
Size is 1200, 1200
Coordinate System is:
GEOGCRS["WGS 84",
    ENSEMBLE["World Geodetic System 1984 ensemble",
        MEMBER["World Geodetic System 1984 (Transit)"],
        MEMBER["World Geodetic System 1984 (G730)"],
        MEMBER["World Geodetic System 1984 (G873)"],
        MEMBER["World Geodetic System 1984 (G1150)"],
        MEMBER["World Geodetic System 1984 (G1674)"],
        MEMBER["World Geodetic System 1984 (G1762)"],
        MEMBER["World Geodetic System 1984 (G2139)"],
        MEMBER["World Geodetic System 1984 (G2296)"],
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]],
        ENSEMBLEACCURACY[2.0]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["Horizontal component of 3D system."],
        AREA["World."],
        BBOX[-90,-180,90,180]],
    ID["EPSG",4326]]
Data axis to CRS axis mapping: 2,1
Origin = (-50.000000000000000,-10.000000000000000)
Pixel Size = (0.008333333333333,-0.008333333333333)
Metadata:
  OVERVIEW_RESAMPLING=NEAREST
  watershed-data-version=v0
  watershed-processing-time=2026-09-23T22:28:35.566482+00:00
  watershed-processing-version=jeffseif/ingest-gpw-livestock
  watershed-product-name=livestock
  watershed-remote-url=git@github.com:cornerstone-data/luc.git
  watershed-source-name=gpw
  AREA_OR_POINT=Area
Image Structure Metadata:
  LAYOUT=COG
  COMPRESSION=DEFLATE
  INTERLEAVE=PIXEL
  OVERVIEW_RESAMPLING=NEAREST
Corner Coordinates:
Upper Left  ( -50.0000000, -10.0000000) ( 50d 0' 0.00"W, 10d 0' 0.00"S)
Lower Left  ( -50.0000000, -20.0000000) ( 50d 0' 0.00"W, 20d 0' 0.00"S)
Upper Right ( -40.0000000, -10.0000000) ( 40d 0' 0.00"W, 10d 0' 0.00"S)
Lower Right ( -40.0000000, -20.0000000) ( 40d 0' 0.00"W, 20d 0' 0.00"S)
Center      ( -45.0000000, -15.0000000) ( 45d 0' 0.00"W, 15d 0' 0.00"S)
Band 1 Block=512x512 Type=Float32, ColorInterp=Gray
  Description = buffalo:heads-per-ha:2000
  Overviews: 600x600, 300x300, 150x150, 75x75
Band 2 Block=512x512 Type=Float32, ColorInterp=Undefined
  Description = buffalo:heads-per-ha:2005
  Overviews: 600x600, 300x300, 150x150, 75x75
Band 3 Block=512x512 Type=Float32, ColorInterp=Undefined
  Description = buffalo:heads-per-ha:2010
  Overviews: 600x600, 300x300, 150x150, 75x75
Band 4 Block=512x512 Type=Float32, ColorInterp=Undefined
  Description = buffalo:heads-per-ha:2020
  Overviews: 600x600, 300x300, 150x150, 75x75
Band 5 Block=512x512 Type=Float32, ColorInterp=Undefined
  Description = cattle:heads-per-ha:2000
  Overviews: 600x600, 300x300, 150x150, 75x75
Band 6 Block=512x512 Type=Float32, ColorInterp=Undefined
  Description = cattle:heads-per-ha:2005
  Overviews: 600x600, 300x300, 150x150, 75x75
Band 7 Block=512x512 Type=Float32, ColorInterp=Undefined
  Description = cattle:heads-per-ha:2010
  Overviews: 600x600, 300x300, 150x150, 75x75
Band 8 Block=512x512 Type=Float32, ColorInterp=Undefined
  Description = cattle:heads-per-ha:2020
  Overviews: 600x600, 300x300, 150x150, 75x75
Band 9 Block=512x512 Type=Float32, ColorInterp=Undefined
  Description = goat:heads-per-ha:2000
  Overviews: 600x600, 300x300, 150x150, 75x75
Band 10 Block=512x512 Type=Float32, ColorInterp=Undefined
  Description = goat:heads-per-ha:2005
  Overviews: 600x600, 300x300, 150x150, 75x75
Band 11 Block=512x512 Type=Float32, ColorInterp=Undefined
  Description = goat:heads-per-ha:2010
  Overviews: 600x600, 300x300, 150x150, 75x75
Band 12 Block=512x512 Type=Float32, ColorInterp=Undefined
  Description = goat:heads-per-ha:2020
  Overviews: 600x600, 300x300, 150x150, 75x75
Band 13 Block=512x512 Type=Float32, ColorInterp=Undefined
  Description = horse:heads-per-ha:2000
  Overviews: 600x600, 300x300, 150x150, 75x75
Band 14 Block=512x512 Type=Float32, ColorInterp=Undefined
  Description = horse:heads-per-ha:2005
  Overviews: 600x600, 300x300, 150x150, 75x75
Band 15 Block=512x512 Type=Float32, ColorInterp=Undefined
  Description = horse:heads-per-ha:2010
  Overviews: 600x600, 300x300, 150x150, 75x75
Band 16 Block=512x512 Type=Float32, ColorInterp=Undefined
  Description = horse:heads-per-ha:2020
  Overviews: 600x600, 300x300, 150x150, 75x75
Band 17 Block=512x512 Type=Float32, ColorInterp=Undefined
  Description = sheep:heads-per-ha:2000
  Overviews: 600x600, 300x300, 150x150, 75x75
Band 18 Block=512x512 Type=Float32, ColorInterp=Undefined
  Description = sheep:heads-per-ha:2005
  Overviews: 600x600, 300x300, 150x150, 75x75
Band 19 Block=512x512 Type=Float32, ColorInterp=Undefined
  Description = sheep:heads-per-ha:2010
  Overviews: 600x600, 300x300, 150x150, 75x75
Band 20 Block=512x512 Type=Float32, ColorInterp=Undefined
  Description = sheep:heads-per-ha:2020
  Overviews: 600x600, 300x300, 150x150, 75x75

Manual inspection

Here is the 2020 GPW grassland (= 1; cultivated rangeland) mask and the GPW cattle density in Mato Grosso, Brazil (10S_060W):

And here is a zoomed in view showing the different overlaps and resolutions:

(vibecoded) comparison with FAOSTAT national totals

Claude code wrote this script for me:

Spot check GPW vs FAOSTAT national totals

"""Spot-check ingested GPW livestock tiles against FAOSTAT stocks, country by country.

For each country: sums heads/ha x hectares-per-pixel over the country's pixels in every ingested
tile it touches, then sets that beside FAOSTAT's stocks (element 5111) for the same species and
year. GPW is scaled to FAOSTAT country by country, so cattle should land near 1.0; sparse sheep
and goats read low, because GPW's integer heads per km2 round them to zero.

Run from the repo root, since Config finds .env by walking up from the working directory:

  uv run python ~/claude-scratch/livestock-luc/spot_check_gpw_livestock.py HND NIC BLZ BRA ARG PRY

Tiles not yet ingested are reported and left out of the sums, which then read low. Pass
--faostat-zip to reuse a local copy of the FAOSTAT bulk archive instead of fetching 32 MiB.
"""

import argparse
import collections
import os
import tempfile

import numpy
import pandas
import rasterio
import rasterio.features
import xarray

from jdluc import config, emit, storage, utils
from jdluc.datasets import faostat_production, gpw_livestock, worldbank_jurisdictions

STOCKS_ELEMENT_CODE = "5111"
STOCKS_ITEM_CODE_TO_SPECIES_NAME = {
    "946": gpw_livestock.Species.BUFFALO.name,
    "866": gpw_livestock.Species.CATTLE.name,
    "1016": gpw_livestock.Species.GOAT.name,
    "1096": gpw_livestock.Species.HORSE.name,
    "976": gpw_livestock.Species.SHEEP.name,
}


def get_hectares(dataset: rasterio.DatasetReader) -> numpy.ndarray:
    lat = dataset.bounds.top - (numpy.arange(dataset.height) + 0.5) * dataset.res[1]
    lon = dataset.bounds.left + (numpy.arange(dataset.width) + 0.5) * dataset.res[0]
    grid = xarray.DataArray(
        numpy.zeros(dataset.shape), coords={"y": lat, "x": lon}, dims=("y", "x")
    )
    return emit.get_hectares_per_pixel(darray=grid).values


def get_gpw_heads(iso_3166s: list[str]) -> dict[tuple[str, str, int], float]:
    national = worldbank_jurisdictions.AdminLevel.NATIONAL
    gdf = worldbank_jurisdictions.get_jurisdiction_for_admin_level(admin_level=national)
    root = config.Config.from_dot_env().ingest_root
    heads: dict[tuple[str, str, int], float] = collections.defaultdict(float)
    for iso_3166 in iso_3166s:
        geometry = gdf.loc[iso_3166].geometry
        tile_ids = worldbank_jurisdictions.get_ten_degree_tile_ids_for_admin_id(
            admin_id=iso_3166, admin_level=national
        )
        for tile_id in sorted(tile_ids):
            uri = storage.join_uri(
                prefix=gpw_livestock.DATASET.get_prefix(tile_id=tile_id), root=root
            )
            if storage.path_exists(uri=uri):
                with rasterio.open(storage.to_gdal_path(uri=uri)) as dataset:
                    mask = rasterio.features.rasterize(
                        [(geometry, 1)],
                        dtype="uint8",
                        fill=0,
                        out_shape=dataset.shape,
                        transform=dataset.transform,
                    ).astype(bool)
                    hectares = get_hectares(dataset=dataset)[mask]
                    for band_idx, (species, year) in enumerate(
                        gpw_livestock.SPECIES_YEARS, start=1
                    ):
                        heads_per_ha = dataset.read(band_idx)[mask]
                        heads[(iso_3166, species.name, year)] += float(
                            numpy.nansum(heads_per_ha * hectares)
                        )
            else:
                print(
                    f"{iso_3166:s}: {tile_id:s} is not ingested; its heads are left out"
                )
    return heads


def get_faostat_stocks(
    iso_3166s: list[str], path_to_zip: str
) -> dict[tuple[str, str, int], float]:
    years = {str(year) for year in gpw_livestock.YEARS}
    stocks: dict[tuple[str, str, int], float] = {}
    for row in faostat_production.iter_rows(path_to_zip=path_to_zip):
        if (
            row["Element Code"] == STOCKS_ELEMENT_CODE
            and row["Item Code"] in STOCKS_ITEM_CODE_TO_SPECIES_NAME
            and row["Year"] in years
            and row["Value"]
            and row["Flag"] != faostat_production.MISSING_FLAG
        ):
            iso_3166 = faostat_production.get_iso_3166(m49_code=row["Area Code (M49)"])
            if iso_3166 in iso_3166s:
                species_name = STOCKS_ITEM_CODE_TO_SPECIES_NAME[row["Item Code"]]
                stocks[(iso_3166, species_name, int(row["Year"]))] = float(row["Value"])
    return stocks


def main() -> int:
    parser = argparse.ArgumentParser(
        description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter
    )
    parser.add_argument("iso_3166s", nargs=argparse.ONE_OR_MORE)
    parser.add_argument(
        "--faostat-zip", help="a local copy of the FAOSTAT bulk archive"
    )
    args = parser.parse_args()
    iso_3166s = [worldbank_jurisdictions.iso_3166_str(s=s) for s in args.iso_3166s]

    heads = get_gpw_heads(iso_3166s=iso_3166s)
    if args.faostat_zip:
        stocks = get_faostat_stocks(iso_3166s=iso_3166s, path_to_zip=args.faostat_zip)
    else:
        with tempfile.TemporaryDirectory() as local_dir:
            path_to_zip = os.path.join(local_dir, "data.zip")
            utils.save_remote_url_to_local_path(
                local_path=path_to_zip,
                params={},
                remote_url=faostat_production.BULK_URL,
            )
            stocks = get_faostat_stocks(iso_3166s=iso_3166s, path_to_zip=path_to_zip)

    names = ["iso_3166", "species", "year"]
    # NB: a left join, so a species FAOSTAT doesn't report shows NaN beside its GPW heads
    df = (
        pandas.Series(heads, name="gpw_heads")
        .rename_axis(index=names)
        .to_frame()
        .join(pandas.Series(stocks, name="faostat_stocks").rename_axis(index=names))
        .sort_index()
        .reset_index()
    )
    df["ratio"] = df["gpw_heads"] / df["faostat_stocks"]
    with pandas.option_context("display.width", 160, "display.max_rows", None):
        print(
            df.to_string(
                formatters={
                    "gpw_heads": "{:,.0f}".format,
                    "faostat_stocks": "{:,.0f}".format,
                    "ratio": "{:.3f}".format,
                },
                index=False,
            )
        )
    return 0


if __name__ == "__main__":
    raise SystemExit(main())

Output for ARG and BRA

iso_3166 species  year   gpw_heads faostat_stocks ratio
     ARG BUFFALO  2000          17            NaN   NaN
     ARG BUFFALO  2005          13            NaN   NaN
     ARG BUFFALO  2010          13            NaN   NaN
     ARG BUFFALO  2020          26            NaN   NaN
     ARG  CATTLE  2000  47,743,509     48,674,400 0.981
     ARG  CATTLE  2005  55,791,987     57,033,528 0.978
     ARG  CATTLE  2010  47,799,558     48,949,744 0.977
     ARG  CATTLE  2020  53,171,717     54,460,799 0.976
     ARG    GOAT  2000   2,716,811      3,490,200 0.778
     ARG    GOAT  2005   3,184,263      4,200,000 0.758
     ARG    GOAT  2010   3,163,577      4,037,036 0.784
     ARG    GOAT  2020   3,555,958      4,695,830 0.757
     ARG   HORSE  2000   3,326,297      3,600,000 0.924
     ARG   HORSE  2005   3,341,833      3,655,000 0.914
     ARG   HORSE  2010   3,400,512      3,600,000 0.945
     ARG   HORSE  2020   2,287,264      2,586,855 0.884
     ARG   SHEEP  2000  12,452,020     13,561,600 0.918
     ARG   SHEEP  2005  14,577,032     15,497,000 0.941
     ARG   SHEEP  2010  14,318,489     15,024,744 0.953
     ARG   SHEEP  2020  13,572,032     14,571,766 0.931
     BRA BUFFALO  2000     818,125      1,102,551 0.742
     BRA BUFFALO  2005     845,861      1,173,629 0.721
     BRA BUFFALO  2010     848,857      1,184,511 0.717
     BRA BUFFALO  2020   1,461,884      1,502,254 0.973
     BRA  CATTLE  2000 167,492,283    169,875,524 0.986
     BRA  CATTLE  2005 204,656,777    207,156,696 0.988
     BRA  CATTLE  2010 206,979,967    209,541,109 0.988
     BRA  CATTLE  2020 215,289,014    217,836,282 0.988
     BRA    GOAT  2000   7,832,884      9,346,813 0.838
     BRA    GOAT  2005   9,170,949     10,306,722 0.890
     BRA    GOAT  2010   8,587,504      9,312,784 0.922
     BRA    GOAT  2020  10,862,679     12,101,686 0.898
     BRA   HORSE  2000   3,594,788      5,831,817 0.616
     BRA   HORSE  2005   3,229,805      5,787,249 0.558
     BRA   HORSE  2010   3,420,421      5,514,253 0.620
     BRA   HORSE  2020   3,648,054      5,959,922 0.612
     BRA   SHEEP  2000  13,233,246     14,784,958 0.895
     BRA   SHEEP  2005  13,969,105     15,588,041 0.896
     BRA   SHEEP  2010  15,717,606     17,380,581 0.904
     BRA   SHEEP  2020  18,733,284     20,623,064 0.908

CATTLE has a >97% ratio; SHEEP >90%; GOAT >70%; HORSE ~90% for ARG but 50-60% for BRA; and BUFFALO is absent from ARG entirely.

juliannedeangelo9

lgtm!

jeffseif added 3 commits September 23, 2026 12:39

use consistent formatting for 'CC BY 4.0' license throughout

13104d8

use EPSG: int's throughout

29a52b3

Add GPW_LIVESTOCK dataset

6d92d1c

jeffseif requested a review from juliannedeangelo9 September 23, 2026 23:17

update docs and citations

78c13ac

jeffseif mentioned this pull request Sep 24, 2026

feat: ingest FAOSTAT livestock production and head data #4

Merged

juliannedeangelo9 approved these changes Sep 24, 2026

View reviewed changes

jeffseif merged commit 35a4cc0 into main Sep 24, 2026
2 checks passed

jeffseif deleted the jeffseif/ingest-gpw-livestock branch September 24, 2026 22:10

jeffseif mentioned this pull request Sep 25, 2026

feat: attribute PASTURE emissions to grazing livestock -- starting with BEEF_CATTLE #5

Open

Navigation Menu

feat: ingest gpw livestock dataset - #3

feat: ingest gpw livestock dataset#3
jeffseif merged 4 commits into
mainfrom
jeffseif/ingest-gpw-livestock

jeffseif commented Sep 23, 2026 •

edited

Loading

Uh oh!

juliannedeangelo9 left a comment

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

Conversation

jeffseif commented Sep 23, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Background

Unrelated commits

Notes

Manual inspection

(vibecoded) comparison with FAOSTAT national totals

Uh oh!

juliannedeangelo9 left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

jeffseif commented Sep 23, 2026 •

edited

Loading