Hola, en este blog hablaremos sobre cómo podemos realizar análisis ráster con facilidad utilizando índices h3.

Objetivo

Para aprender, haremos cálculos para determinar cuántos edificios hay en el área de asentamiento determinada por ESRI Land Cover. Apuntemos a datos a nivel nacional tanto para vectores como para ráster.

Primero busquemos los datos.

Descargar datos ráster

He descargado el área de asentamiento de Esri Land Cover.

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

Descarguemos el año 2023, de un tamaño aproximado de 362 MB

Descargar OSM Edificios de Nepal

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

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

Preprocesar los datos

Apliquemos algo de preprocesamiento a los datos antes de realizar los cálculos reales de la celda h3
Usaremos el programa de línea de comandos gdal para este paso. Instale gdal en su máquina

Conversión a Geotiff optimizado en la nube

Si no conoce cog, consulte aquí: https://www.cogeo.org/

• Comprueba si gdal_translate está disponible

gdal_translate --version

Debería imprimir la versión gdal que estás usando

• Reproyectar ráster a 4326

Su ráster puede tener srs de origen diferentes, cámbielo en consecuencia

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

• Ahora vamos a convertir tif a geotif optimizado para la nube

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

Se tardó aproximadamente un minuto en convertir tiff reproyectado a geotiff

Insertar datos de osm en la tabla postgresql

Estamos usando osm2pgsql para insertar datos de osm en nuestra tabla

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

osm2pgsql tardó 274s (4m 34s) en total.

También puedes usar archivos geojson si tienes alguno que use ogr2ogr

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

ogro2gr tiene una amplia gama de soporte para controladores, por lo que eres bastante flexible en cuanto a tus aportaciones. La salida es la tabla Postgresql

Calcular índices h3

postgresql

Instalar

pip install pgxnclient cmake
pgxn install h3

Crear extensión en tu base de datos

create extension h3;
create extension h3_postgis CASCADE;

Ahora creemos la tabla de edificios

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

Insertar datos a la tabla

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

Registro y sincronización:

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

Ahora calculemos los índices h3 para aquellos edificios que usan centroide. Aquí 8 es la resolución h3 en la que estoy trabajando. Obtenga más información sobre las resoluciones aquí

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

Operaciones ráster

Instalar

pip install h3 h3ronpy rasterio asyncio asyncpg aiohttp

Asegúrate de que el engranaje reproyectado esté en estático/

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

Ejecute el script proporcionado en el repositorio para crear celdas h3 a partir de un ráster. Estoy remuestreando por método de modo: esto depende del tipo de datos que tenga. Para el modo de datos categóricos se adapta mejor. Obtenga más información sobre los métodos de remuestreo aquí

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álisis

Creemos una función para get_h3_indexes en un 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;

Consigamos todos aquellos edificios que estén identificados como área construida en nuestra área de interés

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;

Plan de consulta:

Esto se puede mejorar aún más si se agrega un índice en la columna h3_ix de edificios

create index on buildings(h3_ix);

Al realizar el recuento de disparos: había 24416 edificios en mi área con clase de construcción clasificada según ESRI

Verificación

Verifiquemos si el resultado es verdadero: obtengamos los edificios 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;

Obtengamos celdas h3 también

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

La precisión se puede aumentar después de aumentar la resolución h3 y también dependerá de la entrada y la técnica de remuestreo

Limpieza

Eliminar las tablas que no necesitamos

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: mosaicos vectoriales

Para visualizar los mosaicos, construyamos rápidamente mosaicos vectoriales usando pg_tileserv

• Descargar pg_tileserv Descargue desde el enlace proporcionado arriba o use la ventana acoplable
• Exportar credenciales

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

• Otorgar permiso de selección de mesa

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

• Creemos geometría en índices h3 para visualización (puede hacerlo directamente desde la consulta si está creando mosaicos desde st_asmvt manualmente)

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

• Cree un índice en esa geom h3 para visualizarla de manera efectiva

CREATE INDEX idx_esri_landcover_geometry 
ON esri_landcover 
USING GIST (geometry);

• Correr

  ./pg_tileserv

• Ahora puedes visualizar esos mosaicos MVT según el valor de los mosaicos o de cualquier forma que desees, por ejemplo: ¡maplibre! También puede visualizar el resultado procesado o combinarlo con otros conjuntos de datos. Esta es la visualización que hice en los índices h3 según su valor de celda en QGIS

Repositorio fuente: https://github.com/kshitijrajsharma/raster-analysis-using-h3

Referencias:

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