Uber hndexes 및 PostgreSQL을 사용한 래스터 분석-파이썬 튜토리얼-php.cn

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

목적

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

먼저 데이터를 찾아보자

래스터 데이터 다운로드

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

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

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

Raster Analysis Using Uber hndexes and PostgreSQL

네팔의 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

로그인 후 복사

재투영된 톱니바퀴가 정적 상태인지 확인/

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

로그인 후 복사

분석

폴리곤에서 get_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;

로그인 후 복사

쿼리 계획 :

Raster Analysis Using Uber hndexes and PostgreSQL

건물의 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

로그인 후 복사

Raster Analysis Using Uber hndexes and PostgreSQL

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

로그인 후 복사

Raster Analysis Using Uber hndexes and PostgreSQL

이제 타일 값이나 원하는 방식(예: maplibre)을 기반으로 MVT 타일을 시각화할 수 있습니다! 처리된 결과를 시각화하거나 다른 데이터세트와 결합할 수도 있습니다. 이것은 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/

위 내용은 Uber hndexes 및 PostgreSQL을 사용한 래스터 분석의 상세 내용입니다. 자세한 내용은 PHP 중국어 웹사이트의 기타 관련 기사를 참조하세요!