Olá, Neste blog falaremos sobre como podemos fazer análises Raster com facilidade usando índices h3.

Objetivo

Para aprendizagem, faremos cálculos para descobrir quantos edifícios existem na área de assentamento determinada pela ESRI Land Cover. Vamos buscar dados de nível nacional para vetores e raster.

Vamos primeiro encontrar os dados

Baixar dados raster

Eu baixei a área de assentamento da Esri Land Cover.

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

Vamos baixar o ano de 2023, com tamanho aproximado de 362 MB

Baixe Edifícios OSM do Nepal

Fonte: http://download.geofabrik.de/asia/nepal.html

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

Pré-processar os dados

Vamos aplicar algum pré-processamento aos dados antes dos cálculos reais da célula h3
Estaremos usando o programa de linha de comando gdal para esta etapa. Instale o gdal na sua máquina

Conversão para Geotiff otimizado para nuvem

Se você não conhece cog, confira aqui: https://www.cogeo.org/

• Verifique se gdal_translate está disponível

gdal_translate --version

Deve imprimir a versão gdal que você está usando

• Reprojetar raster para 4326

Seu raster pode ter uma fonte srs diferente, altere-a de acordo

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

• Agora vamos converter tif em geotif otimizado para nuvem

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

Demorou aproximadamente um minuto para converter tiff reprojetado em geotiff

Insira dados osm na tabela postgresql

Estamos usando osm2pgsql para inserir dados osm em nossa tabela

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

osm2pgsql levou 274s (4m 34s) no geral.

Você também pode usar arquivos geojson se tiver algum usando ogr2ogr

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

ogro2gr tem ampla gama de suporte para drivers, então você é bastante flexível sobre qual é sua opinião. A saída é a tabela Postgresql

Calcular índices h3

Postgresql

Instalar

pip install pgxnclient cmake
pgxn install h3

Crie uma extensão em seu banco de dados

create extension h3;
create extension h3_postgis CASCADE;

Agora vamos criar a tabela de edifícios

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

Inserir dados na tabela

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

Registro e tempo:

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

Agora vamos calcular os índices h3 para esses edifícios usando centroid . Aqui 8 é a resolução h3 na qual estou trabalhando. Saiba mais sobre resoluções aqui

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

Operações raster

Instalar

pip install h3 h3ronpy rasterio asyncio asyncpg aiohttp

Certifique-se de que a engrenagem reprojetada esteja em estática/

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

Execute o script fornecido no repositório para criar células h3 a partir de raster. Estou reamostrando pelo método de modo: isso depende do tipo de dados que você possui. Para dados categóricos, o modo se ajusta melhor. Saiba mais sobre métodos de reamostragem aqui

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

Registro :

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

Análise

Vamos criar uma função para get_h3_indexes em um polígono

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;

Vamos pegar todos os prédios identificados como área construída em nossa área de interesse

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;

Plano de consulta:

Isso pode ser melhorado ainda mais se for adicionado o índice na coluna h3_ix de edifícios

create index on buildings(h3_ix);

Quando a contagem de tiros: havia 24.416 edifícios na minha área com classe construída classificada como ESRI

Verificação

Vamos verificar se a saída é verdadeira: vamos obter os edifícios como 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;

Vamos obter células h3 também

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

A precisão pode ser aumentada após aumentar a resolução h3 e também dependerá da entrada e da técnica de reamostragem

Limpar

Abandone as tabelas que não precisamos

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

Opcional - blocos de vetor

Para visualizar os blocos, vamos construir rapidamente blocos vetoriais usando pg_tileserv

• Baixar pg_tileserv Baixe do link fornecido acima ou use o docker
• Exportar credenciais

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

• Conceder permissão de seleção de tabela

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

• Vamos criar geometria em índices h3 para visualização (você pode fazer isso diretamente da consulta se estiver construindo blocos de st_asmvt manualmente)

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

• Crie um índice naquele geom h3 para visualizá-lo de forma eficaz

CREATE INDEX idx_esri_landcover_geometry 
ON esri_landcover 
USING GIST (geometry);

• Correr

  ./pg_tileserv

• Agora você pode visualizar esses blocos MVT com base no valor dos blocos ou da maneira que desejar, por exemplo: maplibre! Você também pode visualizar o resultado processado ou combiná-lo com outros conjuntos de dados. Esta é a visualização que fiz nos índices h3 com base em seu cell_value no QGIS

Repositório de origem: https://github.com/kshitijrajsharma/raster-análise-using-h3

Referências:

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