CROME (Crop Map of England) – ежегодно обновляемый тематический слой землепользования Англии, разрабатываемый Rural Payments Agency (RPA). Продукт построен на основе классификации Random Forest по данным спутников Sentinel‑1 (радарная поляризация VV/VH) и Sentinel‑2 (мультиспектральные каналы RGBNIR). Пространственная модель представлена гексагональными ячейками площадью 0,41 га, разделенная на 46 английских графств. Каждая ячейка содержит уникальный идентификатор, код землепользования (LUCODE), вероятность класса и дату классификации. Данные доступны ежегодно с 2016 по 2024 год после уборочного сезона. **поле с кодом культуры - „lucode“

Работа с CROME GeoParquet на S3

Что понадобится

CROME — ежегодные границы сельскохозяйственных участков Англии. Файлы DuckDB — встраиваемая аналитическая база: отдельный сервер и загрузка всех данных в память не нужны. DuckDB позволяет работать со всеми годами через Python, не поднимая отдельную базу данных.

Установите Python-пакет одним из способов:

uv add duckdb geopandas
# или: python -m pip install duckdb geopandas

Подробности: DuckDB Python API, работа с S3 и пространственные функции.

Подключение ко всему датасету

Сохраните код как crome.py и запустите uv run python crome.py.

import duckdb


def open_crome() -> duckdb.DuckDBPyConnection:
    con = duckdb.connect()
    con.execute("""
        INSTALL httpfs;
        LOAD httpfs;
        INSTALL spatial;
        LOAD spatial;

        CREATE OR REPLACE SECRET esoil_s3 (
            TYPE s3,
            ENDPOINT 's3.esoil.ru',
            URL_STYLE 'path',
            USE_SSL true,
            SCOPE 's3://datasets/'
        );

        CREATE OR REPLACE VIEW crome AS
        SELECT *
        FROM read_parquet(
            's3://datasets/crome/year=*/crome.parquet',
            hive_partitioning = true,
            hive_types = {'year': INTEGER}
        );
    """)
    return con


con = open_crome()
print(con.sql("SELECT min(year), max(year), count(*) AS parcels FROM crome").df())

crome — это view, то есть сохранённое имя для удалённых Parquet-файлов. Создание view не скачивает весь датасет. Результаты .df() — обычные Pandas DataFrame; вызывайте их только после достаточно узкого запроса.

Анализ культур по всем годам

lucode — код культуры. shape_area измеряется в квадратных метрах.

result = con.sql("""
    SELECT year,
           lucode,
           count(*) AS parcels,
           sum(shape_area) / 1_000_000 AS area_km2
    FROM crome
    WHERE lucode IS NOT NULL
    GROUP BY year, lucode
    ORDER BY year, area_km2 DESC
""")

print(result.limit(30).df())

Чтобы получить только 2023 год, меняется только фильтр, а не путь к файлу:

parcels_2023 = con.sql("""
    SELECT cromeid, lucode, county, shape_area
    FROM crome
    WHERE year = 2023 AND county = 'Norfolk'
""")
print(parcels_2023.limit(20).df())

Выборка по ROI-полигону

Предположим, имеется файл study_area.geojson с полигоном территории исследования. Геометрия файла может быть в любой CRS с корректно указанным метаданным; GeoPandas перепроецирует её в EPSG:27700 для CROME. Сначала применяется дешёвый фильтр bbox, затем точный ST_Intersects по WKB.

import geopandas as gpd

roi_file = gpd.read_file("study_area.geojson").to_crs("EPSG:27700")
roi_geometry = roi_file.geometry.union_all()
minx, miny, maxx, maxy = roi_geometry.bounds
roi_wkt = roi_geometry.wkt

roi = con.execute(
    """
    SELECT year, cromeid, lucode, county, shape_area,
           ST_AsText(ST_GeomFromWKB(geometry)) AS geometry_wkt
    FROM crome
    WHERE year BETWEEN ? AND ?
      AND geometry IS NOT NULL
      AND bbox.xmin <= ? AND bbox.xmax >= ?
      AND bbox.ymin <= ? AND bbox.ymax >= ?
      AND ST_Intersects(ST_GeomFromWKB(geometry), ST_GeomFromText(?))
    """,
    [2022, 2024, maxx, minx, maxy, miny, roi_wkt],
).fetchdf()

print(roi)

Для чисто табличных задач не выбирайте geometry: WKB — самый тяжёлый столбец.