Карта создана лабораторией Global Land Analysis&Discovery (GLAD) под руководством Мэрилендского университета (США) на основе спутниковых снимков Sentinel‑2 (10 м) и PlanetScope (3 м), а также полевых данных с 450 проверочных пикселей и 18 тыс. точек наблюдений. Алгоритм включает обучение двух моделей Random Forest, блочное картографирование и калибровку вероятностных порогов, чтобы площади по карте совпали с несмещёнными оценками по выборке. Итог — первая открытая 10-метровая карта Китая для кукурузы и сои с общей точностью 91,8 %; оценки площадей согласуются с официальной статистикой на уровне провинций. Карта доступна в Google Earth Engine. Расшифровка пикселей: class_names={ 1: «corn», 2: «soybean» }

Работа с China Crop Map в Zarr на S3

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

China Crop Map — категориальная карта культур Китая за 2019 год. Zarr хранит переменную china_crop_map на исходной сетке EPSG:4326 с размером пикселя около 30 м. Значение 0 означает отсутствие покрытия источником.

Это один снимок за 2019 год, а не временной куб: массив имеет только измерения y и x, без year и date.

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

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

Или создайте обычное виртуальное окружение и установите те же пакеты через python -m pip.

3. Подключение к данным на S3

import s3fs
import xarray as xr

fs = s3fs.S3FileSystem(
    anon=True,
    client_kwargs={"endpoint_url": "https://s3.esoil.ru"},
)
ds = xr.open_zarr(
    fs.get_mapper("s3://datasets/china_crop_map"),
    chunks="auto",
    consolidated=None,
)

consolidated=None позволяет Xarray прочитать фактическую конфигурацию Zarr v3, не ожидая старого .zmetadata. Открытие набора не загружает всю карту в память: данные будут прочитаны при .compute(), .load(), .values или построении графика.

4. Быстрый осмотр

classes = ds["china_crop_map"]

print(ds)
print(classes.dims)                   # ('y', 'x')
print(ds.attrs["temporal_coverage"])  # 2019

Не используйте .sel(year=2019): такого измерения в наборе нет. Выбирайте пространственное окно непосредственно из classes.

5. Выбор области и чтение результата

Координаты x и y — долгота и широта в EPSG:4326. Если область исследования записана в другой CRS, преобразуйте её в EPSG:4326.

import geopandas as gpd

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

# y уменьшается от севера к югу.
subset = classes.sel(x=slice(minx, maxx), y=slice(maxy, miny))
subset = subset.compute()
print(subset.shape)

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

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

Исходный GeoTIFF не содержит встроенной легенды. Анализируйте коды как численные значения, пока внешняя легенда не будет добавлена отдельно.

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

Ноль исключается из статистики: это no-data, а не класс культуры.

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

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

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

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

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

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

subset.rio.to_raster("china_crop_map_2019.tif", compress="ZSTD")

Для категориальных кодов используйте дискретную палитру. Без внешней легенды нельзя надёжно сопоставить числовые коды названиям культур.

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

• Карта описывает только 2019 год, но не имеет временного измерения.
• Выполняйте пространственный выбор в EPSG:4326.
• Исключайте 0 из статистики классов.
• Ограничивайте запрос пространственным окном до .compute().