안녕하세요. 이번 블로그에서는 h3 인덱스를 사용하여 래스터 분석을 쉽게 수행하는 방법에 대해 설명하겠습니다.

목적

학습을 위해 ESRI Land Cover가 결정한 거주 지역에 건물이 몇 개 있는지 계산해 보겠습니다. 벡터와 래스터 모두에 대한 국가 수준의 데이터를 목표로 삼겠습니다.

먼저 데이터를 찾아보자

래스터 데이터 다운로드

Esri Land Cover에서 정착지 지역을 다운로드했습니다.

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

약 362MB 크기의 2023년을 다운로드합니다.

네팔의 OSM 건물 다운로드

출처 : http://download.geofabrik.de/asia/nepal.html

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

데이터 전처리

실제 h3 셀 계산 전에 데이터에 일부 전처리를 적용해 보겠습니다.
이 단계에서는 gdal 명령줄 프로그램을 사용합니다. 컴퓨터에 gdal을 설치하세요

클라우드에 최적화된 Geotiff로 변환

코그에 대해 모르신다면 여기에서 확인하세요: 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를 클라우드에 최적화된 geotif로 변환할 수 있습니다.

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

재투영된 TIFF를 Geotiff로 변환하는 데 약 1분이 걸렸습니다.

postgresql 테이블에 osm 데이터 삽입

우리는 osm2pgsql을 사용하여 테이블에 osm 데이터를 삽입하고 있습니다.

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

osm2pgsql은 전체적으로 274초(4분 34초)가 걸렸습니다.

ogr2ogr을 사용하는 경우 geojson 파일도 사용할 수 있습니다.

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

ogro2gr은 드라이버에 대한 광범위한 지원을 제공하므로 입력 내용에 대해 매우 유연합니다. 출력은 Postgresql 테이블입니다.

h3 인덱스 계산

포스트그레SQL

설치하다

pip install pgxnclient cmake
pgxn install h3

DB에 확장 프로그램 만들기

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

이제 centroid 를 사용하여 해당 건물의 h3 지수를 계산해 보겠습니다. 여기 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

재투영된 톱니바퀴가 static/
상태인지 확인하세요.

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);

촬영 횟수: 내 지역에는 ESRI로 분류된 건축 등급의 건물이 24416개 있었습니다.

확인

출력이 true인지 확인합니다. 건물을 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! 처리된 결과를 시각화하거나 다른 데이터세트와 결합할 수도 있습니다. 이것은 QGIS의 cell_value를 기반으로 h3 인덱스에 대해 수행한 시각화입니다.

소스 저장소 : https://github.com/kshitijrajsharma/raster-analytic-using-h3

참고자료 :

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