Данные созданы командой дистанционного зондирования Земли (Earth Observation Team) Научно-технического филиала (STB) Министерства сельского хозяйства и агропродовольствия Канады (AAFC). Работа ведётся с 2009 года в рамках операционной системы картографирования типов посевов. Исходные данные: оптические данные (Landsat-5/7/8/9, Sentinel‑2, AWiFS, DMC, SPOT, Gaofen‑1) и радарные данные (RADARSAT‑2, с 2021 г. – RADARSAT Constellation Mission), обучение и валидация проводится на данных провинциальных страховых компаний (самый точный источник), полевые наблюдения, собранные сотрудниками AAFC, региональными партнёрами и подрядчиками (десятки тысяч точек ежегодно), а также референсные слои земного покрова для несельскохозяйственных классов. Классификация проводится алгоритмом Деревья Решений с точностью более 85%. Данные обновляются ежегодно и разделены по провинциям.

Работа с Annual Crop Inventory (ACI) в Zarr на S3

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

Это ежегодная канадская карта сельскохозяйственных культур Annual Crop Inventory (ACI). Исходные данные — региональные GeoTIFF с категориальными кодами классов. После конвертации они образуют единый Zarr v3-куб с осями (year, y, x).

В наборе две переменные типа uint8:

Переменная	Значение
aci	Выбранный код класса ACI.
aci_alternative	Самый сильный другой класс в зоне перекрытия источников; 0 означает, что альтернативы нет.

Значение 0 в aci означает отсутствие исходных данных (fill value). Смысл всех ненулевых кодов задаёт легенда в атрибутах переменных: flag_values, flag_meanings и class_labels.

Куб использует одну общую сетку для всех лет: CRS первого доступного 30-метрового источника, разрешение 30 м и единое выравнивание пикселей. Метаданные координат и grid mapping записаны в CF-совместимом виде.

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

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

Открытие датасета не скачивает весь ACI в память. Данные читаются, когда нужен результат: например, после .compute(), .load(), получения .values или построения графика.

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

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

uv venv
source .venv/bin/activate        # Linux/macOS
# .venv\Scripts\activate         # Windows PowerShell
uv pip install xarray rioxarray zarr s3fs dask

pip

python -m venv .venv
source .venv/bin/activate        # Linux/macOS
pip install xarray rioxarray zarr s3fs dask

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

import s3fs
import xarray as xr

s3_endpoint_url = "https://s3.esoil.ru"
store_url = "s3://datasets/aci.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 позволяет Xarray использовать фактическую конфигурацию опубликованного Zarr, не предполагая наличие старого .zmetadata. После этой команды ds — описание данных и отложенный план чтения, а не полный массив в памяти.

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

print(ds)
print(ds.data_vars)       # aci, aci_alternative
print(ds.coords)          # year, y, x
print(ds.attrs)

aci = ds["aci"]
print(aci.attrs["flag_values"])
print(aci.attrs["class_labels"])

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

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

Сначала ограничьте год и небольшое окно в координатах CRS куба. Направление оси y определяется affine-преобразованием; обычно для северной точки первой задают большее значение y.

window = ds.sel(
    year=2024,
    x=slice(500_000, 530_000),
    y=slice(5_600_000, 5_570_000),
)

classes = window["aci"]
alternatives = window["aci_alternative"]

Числа выше — только пример координат. Перед реальной выборкой посмотрите ds.x, ds.y и CRS/grid mapping набора, а не подставляйте широту и долготу вслепую.

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

classes_loaded = classes.compute()
# или, когда массив действительно нужен в NumPy:
values = classes.values

Подсчёт покрытия и классов

valid = classes.where(classes != 0)
class_counts = valid.to_series().value_counts().sort_index().compute()
print(class_counts)

Значение 0 следует исключать из статистики классов: это no-data, а не класс культуры.

7. Перекрытия источников

aci_alternative полезна для аудита участков, где региональные источники дали разные ненулевые классы.

conflicts = window.where(
    (window["aci"] != 0) & (window["aci_alternative"] != 0)
)

conflict_pixels = (conflicts["aci_alternative"] != 0).sum().compute()
print(f"Пикселей с альтернативным классом: {int(conflict_pixels)}")
print(ds.attrs.get("conflict_resolution"))

Для выбора aci и aci_alternative был использован следующий алгоритм:

1. Для каждого выходного пикселя брали все ACI-источники, которые его покрывают.
2. Рассматривали только ненулевые коды классов из всех источников.
3. Для каждого такого кода в каждом источнике считали вес: сколько раз этот код встречается в окрестности 5×5 пикселей вокруг данного пикселя.
4. Вес каждого кода суммировали между всеми источниками.
5. В aci записывали класс с наибольшим суммарным весом.
6. В aci_alternative записывали второй по весу класс или 0, если пересечений не было.
7. При одинаковом весе побеждал первый встреченный класс.

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

Все обычные операции Xarray можно строить лениво. Например, доля валидных пикселей по годам:

coverage_by_year = (ds["aci"] != 0).mean(dim=("y", "x"))
print(coverage_by_year.compute())

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

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

Сначала загрузите требуемое окно, затем примените дискретную палитру, согласованную с class_labels.

small = ds["aci"].sel(year=2024).isel(x=slice(0, 1000), y=slice(0, 1000))
small.compute().plot()

Для экспорта в GeoTIFF вызовите метод .rio.to_raster. rioxarray использует CRS и transform, записанные в Zarr. В примере ниже экспортируются значения всего 2024 года:

ds["aci"].sel(year=2024).rio.to_raster("aci_2024.tif")

Для уменьшения размера сохраняемого файла используйте сжатие:

ds["aci"].sel(year=2024).rio.to_raster("aci_2024.tif", compress="ZSTD")

Перед экспортом проверьте, что у выбранного DataArray не осталось измерения year; если выбрано несколько лет, экспортируйте их по одному.

Если вы получаете ошибку

AttributeError: 'DataArray' object has no attribute 'rio'

значит вы не импортировали модуль rioxarray:

import rioxarray

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

1. Сначала проверяйте координаты, CRS и доступные годы через ds.coords.
2. Для статистики классов всегда явно исключайте код 0.
3. Ограничивайте запросы по year, x и y до вызова .compute().
4. Используйте aci_alternative только для анализа конфликтов мозаики.