Cropland Data Layer (CDL) — ежегодный общедоступный категориальный растровый продукт Национальной службы сельскохозяйственной статистики Министерства сельского хозяйства США (USDA NASS). Он описывает сельскохозяйственные культуры и другие категории земного покрова США на уровне пикселей. В этой коллекции содержится по одному одноканальному GeoTIFF на каждый год 2008-2025; все файлы имеют пространственное разрешение 30 м и проекцию NAD83 / Conus Albers (EPSG:5070). Значения являются кодами классов, а не измеряемой непрерывной величиной; код 0 обозначает фон или отсутствие данных. Методика CDL и входные данные менялись по годам, поэтому точность и классификатор следует брать из метаданных соответствующего годового выпуска, а не считать едиными для всей временной серии.

Работа с USDA Cropland Data Layer (CDL) в Zarr на S3

1. Что это за данные

USDA National Agricultural Statistics Service Cropland Data Layer — ежегодная тематическая карта земного покрова и сельскохозяйственных культур США. В опубликованном Zarr переменная cdl хранит коды классов для измерений year, y, x на сетке EPSG:5070 (NAD83 / Conus Albers). Ноль означает фон или отсутствие данных, а не отдельный класс.

В хранилище также есть ovr_2x/cdl: двухкратный обзор для быстрого просмотра. Для статистики классов и экспорта используйте исходную переменную cdl, поскольку категории нельзя корректно усреднять.

2. Zarr, Xarray и ленивое чтение

Zarr хранит массивы чанками, поэтому можно выбрать год и пространственное окно до чтения значений. xarray сохраняет имена координат, CRS и легенду из метаданных.

Открытие куба не загружает его целиком: чтение начинается при .compute(), .load(), обращении к .values или визуализации.

3. Подготовка окружения

uv (рекомендуется)

uv venv
source .venv/bin/activate
uv pip install "xarray[io]" dask s3fs rioxarray matplotlib

pip

python -m venv .venv
source .venv/bin/activate
python -m pip install "xarray[io]" dask s3fs rioxarray matplotlib

4. Подключение к CDL на S3

import s3fs
import xarray as xr

s3_endpoint_url = "https://s3.esoil.ru"
store_url = "s3://datasets/cdl.zarr"

fs = s3fs.S3FileSystem(
    anon=True,
    client_kwargs={"endpoint_url": s3_endpoint_url},
)

ds = xr.open_zarr(
    fs.get_mapper(store_url),
    chunks="auto",
    consolidated=None,
)

consolidated=None не предполагает наличие старого .zmetadata и позволяет Xarray прочитать фактическую конфигурацию Zarr v3. ds пока содержит только метаданные и отложенный план чтения.

5. Быстрый осмотр и легенда

print(ds)
print(ds.year.values)

classes = ds["cdl"]
print(classes.attrs["long_name"])

labels = dict(zip(classes.attrs["flag_values"], classes.attrs["flag_meanings"].split()))
print(labels[1])

flag_values и flag_meanings задают легенду. В CDL имена классов разделены пробелами и используют подчёркивания, например corn или developed_low_intensity. Не добавляйте ноль в эту легенду: он обозначает фон или отсутствующие данные.

6. Базовые операции

Выбор года и области

Координаты имеют единицы метров в EPSG:5070. Сначала преобразуйте границы исследуемой области в эту CRS, затем выберите окно. Направление координаты y может быть убывающим, поэтому границы передаются через min и max.

import geopandas as gpd

study_area = gpd.read_file("study_area.geojson").to_crs("EPSG:5070")
minx, miny, maxx, maxy = study_area.union_all().bounds

cdl_2024 = classes.sel(
    year=2024,
    x=slice(minx, maxx),
    y=slice(maxy, miny),
)

Если окно оказалось пустым, проверьте порядок y в ds.y.values[:2] и поменяйте границы в slice местами.

Чтение результата

cdl_2024 = cdl_2024.compute()
print(cdl_2024.shape)

В память попадает только выбранный год и окно, а не весь многолетний куб.

Подсчёт классов

valid = cdl_2024.where(cdl_2024 != 0)
class_counts = valid.to_series().value_counts().sort_index()
class_counts.index = class_counts.index.astype(int)
print(class_counts)

7. Двухкратный обзор

ovr_2x/cdl предназначен для интерактивного обзора и быстрого чернового графика. Он содержит те же измерения и не заменяет основной массив для площадей, подсчёта пикселей или анализа мелких объектов.

overview_ds = xr.open_zarr(
    fs.get_mapper(store_url),
    group="ovr_2x",
    chunks="auto",
    consolidated=None,
)
overview = overview_ds["cdl"].sel(year=2024)
overview.isel(x=slice(0, 2000), y=slice(0, 2000)).compute().plot()

8. Ленивые вычисления и Dask

Например, долю пикселей с данными по годам можно получить одним отложенным вычислением:

coverage_by_year = (classes != 0).mean(dim=("y", "x"))
print(coverage_by_year.compute())

Для больших задач сначала выбирайте годы и область, затем стройте выражение и вызывайте .compute() один раз. Не вызывайте .values для всего куба.

9. Визуализация и экспорт

import rioxarray  # регистрирует accessor .rio

small = classes.sel(year=2024).isel(x=slice(0, 1000), y=slice(0, 1000))
small.compute().plot()

cdl_2024.rio.to_raster("cdl_2024.tif", compress="ZSTD")

Экспортируйте один выбранный год за раз: у DataArray не должно остаться измерения year. Для карт классов применяйте дискретную палитру, а не непрерывную шкалу.

10. Практические рекомендации

• Отбрасывайте код 0 из статистики классов.
• Выполняйте пространственный выбор в EPSG:5070.
• Используйте ovr_2x/cdl только для просмотра; для анализа берите cdl.
• Ограничивайте запрос годом и пространственным окном до .compute().