Привет! В этом блоге мы поговорим о том, как можно легко выполнить растровый анализ с использованием индексов h3.

Цель

Для обучения мы выполним расчет, чтобы выяснить, сколько зданий находится на территории поселения, определенной ESRI Land Cover. Давайте нацелимся на данные национального уровня как для векторных, так и для растровых данных.

Давайте сначала найдем данные

Загрузить растровые данные

Я загрузил территорию поселения с сайта Esri Land Cover.

• https://livingatlas.arcgis.com/landcover/

Давайте скачаем 2023 год, размером около 362 МБ

Скачать OSM Здания Непала

Источник: http://download.geofabrik.de/asia/nepal.html

wget http://download.geofabrik.de/asia/nepal-latest.osm.pbf

Предварительная обработка данных

Давайте применим некоторую предварительную обработку к данным перед фактическими вычислениями ячеек h3.
Для этого шага мы будем использовать программу командной строки gdal. Установите gdal на свой компьютер

Преобразование в оптимизированный для облака Geotiff

Если вы не знаете о cog, оформите заказ здесь: https://www.cogeo.org/

• Проверьте, доступен ли gdal_translate

gdal_translate --version

Он должен напечатать версию gdal, которую вы используете

• Перепроектировать растр до 4326

Ваш растр может иметь другой источник srs, измените его соответствующим образом

gdalwarp esri-settlement-area-kathmandu-grid.tif esri-landcover-4326.tif -s_srs EPSG:32645 -t_srs EPSG:4326

• Теперь давайте конвертируем tif в геотиф, оптимизированный для облака.

gdal_translate -of COG esri-landcover-4326.tif esri-landcover-cog.tif

Преобразование перепроецированного tiff в геотифф заняло около минуты

Вставьте данные osm в таблицу postgresql

Мы используем osm2pgsql для вставки данных osm в нашу таблицу

osm2pgsql --create nepal-latest.osm.pbf -U postgres

osm2pgsql в целом заняло 274 секунды (4 минуты 34 секунды).

Вы также можете использовать файлы geojson, если они у вас есть, используя ogr2ogr

ogr2ogr -f PostgreSQL  PG:"dbname=postgres user=postgres password=postgres" buildings_polygons_geojson.geojson -nln buildings

ogro2gr имеет широкий спектр поддержки драйверов, поэтому вы можете достаточно гибко выбирать вводные данные. Вывод — таблица Postgresql

Рассчитать индексы h3

Постгреск

Установить

pip install pgxnclient cmake
pgxn install h3

Создайте расширение в своей базе данных

create extension h3;
create extension h3_postgis CASCADE;

Теперь создадим таблицу зданий

CREATE TABLE buildings (
  id SERIAL PRIMARY KEY,
  osm_id BIGINT,
  building VARCHAR,
  geometry GEOMETRY(Polygon, 4326)
);

Вставить данные в таблицу

INSERT INTO buildings (osm_id, building, geometry)
SELECT osm_id, building, way
FROM planet_osm_polygon pop
WHERE building IS NOT NULL;

Журнал и время:

Updated Rows    8048542
Query   INSERT INTO buildings (osm_id, building, geometry)
    SELECT osm_id, building, way
    FROM planet_osm_polygon pop
    WHERE building IS NOT NULL
Start time  Mon Aug 12 08:23:30 NPT 2024
Finish time Mon Aug 12 08:24:25 NPT 2024

Теперь давайте рассчитаем индексы h3 для этих зданий, используя centroid . Здесь 8 — разрешение h3, над которым я работаю. Подробнее о разрешениях можно узнать здесь

ALTER TABLE buildings ADD COLUMN h3_index h3index GENERATED ALWAYS AS (h3_lat_lng_to_cell(ST_Centroid(geometry), 8)) STORED;

Растровые операции

Установить

pip install h3 h3ronpy rasterio asyncio asyncpg aiohttp

Убедитесь, что перепроецируемый винтик находится в статическом состоянии/

mv esri-landcover-cog.tif ./static/

Запустите скрипт из репозитория, чтобы создать ячейки h3 из растра . Я выполняю повторную выборку по методу режима: это зависит от типа имеющихся у вас данных. Для категориальных данных режим подходит лучше. Подробнее о методах повторной выборки можно узнать здесь

python cog2h3.py --cog esri-landcover-cog.tif --table esri_landcover --res 8 --sample_by mode

Бревно :

2024-08-12 08:55:27,163 - INFO - Starting processing
2024-08-12 08:55:27,164 - INFO - COG file already exists: static/esri-landcover-cog.tif
2024-08-12 08:55:27,164 - INFO - Processing raster file: static/esri-landcover-cog.tif
2024-08-12 08:55:41,664 - INFO - Determined Min fitting H3 resolution: 13
2024-08-12 08:55:41,664 - INFO - Resampling original raster to : 1406.475763m
2024-08-12 08:55:41,829 - INFO - Resampling Done
2024-08-12 08:55:41,831 - INFO - New Native H3 resolution: 8
2024-08-12 08:55:41,967 - INFO - Converting H3 indices to hex strings
2024-08-12 08:55:42,252 - INFO - Raster calculation done in 15 seconds
2024-08-12 08:55:42,252 - INFO - Creating or replacing table esri_landcover in database
2024-08-12 08:55:46,104 - INFO - Table esri_landcover created or updated successfully in 3.85 seconds.
2024-08-12 08:55:46,155 - INFO - Processing completed

Анализ

Давайте создадим функцию для получения_h3_indexes в многоугольнике

CREATE OR REPLACE FUNCTION get_h3_indexes(shape geometry, res integer)
  RETURNS h3index[] AS $$
DECLARE
  h3_indexes h3index[];
BEGIN
  SELECT ARRAY(
    SELECT h3_polygon_to_cells(shape, res)
  ) INTO h3_indexes;

  RETURN h3_indexes;
END;
$$ LANGUAGE plpgsql IMMUTABLE;

Давайте получим все те здания, которые определены как застроенная территория в интересующей нас области

WITH t1 AS (
  SELECT *
  FROM esri_landcover el
  WHERE h3_ix = ANY (
    get_h3_indexes(
      ST_GeomFromGeoJSON('{
        "coordinates": [
          [
            [83.72922006065477, 28.395029869336483],
            [83.72922006065477, 28.037312312532066],
            [84.2367635433626, 28.037312312532066],
            [84.2367635433626, 28.395029869336483],
            [83.72922006065477, 28.395029869336483]
          ]
        ],
        "type": "Polygon"
      }'), 8
    )
  ) AND cell_value = 7
)
SELECT *
FROM buildings bl
JOIN t1 ON bl.h3_ix = t1.h3_ix;

План запроса:

Это можно улучшить, если добавить индекс в столбец зданий h3_ix

create index on buildings(h3_ix);

При подсчете съемки: в моем районе было 24416 зданий с классом постройки, присвоенным ESRI

Проверка

Давайте проверим, верен ли вывод: получим здания в формате geojson

WITH t1 AS (
  SELECT *
  FROM esri_landcover el
  WHERE h3_ix = ANY (
    get_h3_indexes(
      ST_GeomFromGeoJSON('{
        "coordinates": [
          [
            [83.72922006065477, 28.395029869336483],
            [83.72922006065477, 28.037312312532066],
            [84.2367635433626, 28.037312312532066],
            [84.2367635433626, 28.395029869336483],
            [83.72922006065477, 28.395029869336483]
          ]
        ],
        "type": "Polygon"
      }'), 8
    )
  ) AND cell_value = 7
)
SELECT jsonb_build_object(
  'type', 'FeatureCollection',
  'features', jsonb_agg(ST_AsGeoJSON(bl.*)::jsonb)
)
FROM buildings bl
JOIN t1 ON bl.h3_ix = t1.h3_ix;

Давайте тоже получим ячейки h3

with t1 as (
  SELECT *, h3_cell_to_boundary_geometry(h3_ix)
  FROM esri_landcover el
  WHERE h3_ix = ANY (
    get_h3_indexes(
      ST_GeomFromGeoJSON('{
        "coordinates": [
          [
            [83.72922006065477, 28.395029869336483],
            [83.72922006065477, 28.037312312532066],
            [84.2367635433626, 28.037312312532066],
            [84.2367635433626, 28.395029869336483],
            [83.72922006065477, 28.395029869336483]
          ]
        ],
        "type": "Polygon"
      }'), 8
    )
  ) AND cell_value = 7
)
SELECT jsonb_build_object(
  'type', 'FeatureCollection',
  'features', jsonb_agg(ST_AsGeoJSON(t1.*)::jsonb)
)
FROM t1

Точность может быть увеличена после увеличения разрешения h3, а также будет зависеть от ввода и метода повторной выборки.

Очистка

Удалить ненужные нам таблицы

drop table planet_osm_line;
drop table planet_osm_point;
drop table planet_osm_polygon;
drop table planet_osm_roads;
drop table osm2pgsql_properties;

Необязательно — векторные плитки

Для визуализации тайлов позволяет быстро создавать векторные тайлы с помощью pg_tileserv

• Загрузить pg_tileserv Загрузите ссылку, указанную выше, или используйте docker
• Экспорт учетных данных

export DATABASE_URL=postgresql://postgres:postgres@localhost:5432/postgres

• Предоставить нашей таблице разрешение на выбор

GRANT SELECT ON buildings to postgres;
GRANT SELECT ON esri_landcover to postgres;

• Давайте создадим геометрию по индексам h3 для визуализации (вы можете сделать это непосредственно из запроса, если вы создаете тайлы из st_asmvt вручную)

ALTER TABLE esri_landcover 
ADD COLUMN geometry geometry(Polygon, 4326) 
GENERATED ALWAYS AS (h3_cell_to_boundary_geometry(h3_ix)) STORED;

• Создайте индекс для этой геометрии h3, чтобы эффективно ее визуализировать.

CREATE INDEX idx_esri_landcover_geometry 
ON esri_landcover 
USING GIST (geometry);

• Бегать

  ./pg_tileserv

• Теперь вы можете визуализировать эти плитки MVT на основе значения плитки или любым удобным для вас способом, например: Maplibre! Вы также можете визуализировать обработанный результат или объединить его с другими наборами данных. Это визуализация, которую я сделал для индексов h3 на основе их значения cell_value в QGIS.

Репозиторий исходного кода: https://github.com/kshitijrajsharma/raster-anaанализ-using-h3

Ссылки:

• https://blog.rustprooflabs.com/2022/04/postgis-h3-intro
• https://jsonsingh.com/blog/uber-h3/
• https://h3ronpy.readthedocs.io/en/latest/