Bonjour, Dans ce blog, nous expliquerons comment effectuer facilement une analyse raster à l'aide des index h3.

Objectif

Pour l'apprentissage, nous effectuerons des calculs pour déterminer le nombre de bâtiments présents dans la zone de peuplement déterminée par ESRI Land Cover. Visons les données au niveau national pour les vecteurs et les raster.

Trouvons d'abord les données

Télécharger des données raster

J'ai téléchargé la zone de peuplement depuis Esri Land Cover.

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

Téléchargeons l'année 2023, d'une taille d'environ 362 Mo

Télécharger Bâtiments OSM du Népal

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

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

Prétraiter les données

Appliquons un prétraitement aux données avant les calculs réels des cellules h3
Nous utiliserons le programme de ligne de commande gdal pour cette étape. Installez gdal sur votre machine

Conversion vers Geotiff optimisé pour le cloud

Si vous ne connaissez pas Cog, commandez ici : https://www.cogeo.org/

• Vérifiez si gdal_translate est disponible

gdal_translate --version

Il devrait imprimer la version gdal que vous utilisez

• Reprojeter le raster vers 4326

Votre raster peut avoir des sources sources différentes, modifiez-le en conséquence

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

• Convertissons maintenant le tif en géotif optimisé pour le cloud

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

Il a fallu environ une minute pour convertir le tiff reprojeté en géotiff

Insérer les données osm dans la table postgresql

Nous utilisons osm2pgsql pour insérer des données osm dans notre table

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

osm2pgsql a pris 274 secondes (4 min 34 s) au total.

Vous pouvez également utiliser des fichiers geojson si vous en avez en utilisant ogr2ogr

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

ogro2gr propose une large gamme de supports pour les pilotes, vous êtes donc assez flexible quant à la nature de votre entrée. La sortie est une table Postgresql

Calculer les indices h3

Postgresql

Installer

pip install pgxnclient cmake
pgxn install h3

Créez une extension dans votre base de données

create extension h3;
create extension h3_postgis CASCADE;

Créons maintenant la table des bâtiments

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

Insérer des données dans le tableau

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

Journal et timing :

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

Calculons maintenant les indices h3 pour ces bâtiments en utilisant centroid . Ici, 8 est la résolution h3 sur laquelle je travaille. En savoir plus sur les résolutions ici

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

Opérations raster

Installer

pip install h3 h3ronpy rasterio asyncio asyncpg aiohttp

Assurez-vous que le rouage reprojeté est en statique/

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

Exécutez le script fourni dans le référentiel pour créer des cellules h3 à partir du raster. Je rééchantillonne par méthode de mode : cela dépend du type de données dont vous disposez. Pour les données catégorielles, le mode convient mieux. En savoir plus sur les méthodes de rééchantillonnage ici

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

Enregistrer :

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

Analyse

Créons une fonction pour get_h3_indexes dans un polygone

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;

Obtenons tous les bâtiments identifiés comme zone construite dans notre zone d'intérêt

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 requête :

Cela peut encore être amélioré si un index est ajouté sur la colonne h3_ix des bâtiments

create index on buildings(h3_ix);

Lors du décompte des tirs : il y avait 24416 bâtiments dans ma région avec une classe de construction classée selon ESRI

Vérification

Vérifions si le résultat est vrai : obtenons les bâtiments au format 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;

Obtenons également les cellules 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

La précision peut être augmentée après avoir augmenté la résolution h3 et dépendra également de la technique d'entrée et de rééchantillonnage

Nettoyage

Supprimez les tables dont nous n'avons pas besoin

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

Facultatif - Tuiles vectorielles

Pour visualiser les tuiles, construisons rapidement des tuiles vectorielles à l'aide de pg_tileserv

• Télécharger pg_tileserv Téléchargez à partir du lien fourni ci-dessus ou utilisez Docker
• Exporter les identifiants

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

• Accordez à notre autorisation de sélection de table

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

• Créons une géométrie sur les index h3 pour la visualisation (vous pouvez le faire directement à partir d'une requête si vous créez manuellement des tuiles à partir de st_asmvt)

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

• Créez un index sur cette géom h3 pour la visualiser efficacement

CREATE INDEX idx_esri_landcover_geometry 
ON esri_landcover 
USING GIST (geometry);

• Courir

  ./pg_tileserv

• Vous pouvez maintenant visualiser ces tuiles MVT en fonction de la valeur des tuiles ou comme vous le souhaitez, par exemple : maplibre ! Vous pouvez également visualiser le résultat traité ou le combiner avec d'autres ensembles de données. C'est la visualisation que j'ai faite sur les index h3 en fonction de leur cell_value dans QGIS

Repo source : https://github.com/kshitijrajsharma/raster-analysis-using-h3

Références :

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