Hallo, in diesem Blog werden wir darüber sprechen, wie wir mit h3-Indizes ganz einfach Rasteranalysen durchführen können.

Objektiv

Zum Lernen werden wir eine Berechnung durchführen, um herauszufinden, wie viele Gebäude sich in der von ESRI Land Cover ermittelten Siedlungsfläche befinden. Ziel ist es, Daten auf nationaler Ebene sowohl für Vektor- als auch für Rasterdaten zu erhalten.

Lassen Sie uns zunächst die Daten finden

Laden Sie Rasterdaten herunter

Das Siedlungsgebiet habe ich von Esri Land Cover heruntergeladen.

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

Laden wir das Jahr 2023 mit einer Größe von ca. 362 MB herunter

Laden Sie OSM Buildings of Nepal herunter

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

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

Verarbeiten Sie die Daten vor

Lassen Sie uns vor den eigentlichen h3-Zellenberechnungen eine Vorverarbeitung auf die Daten anwenden
Für diesen Schritt verwenden wir das Befehlszeilenprogramm gdal. Installieren Sie gdal auf Ihrem Computer

Konvertierung in Cloud-optimiertes Geotiff

Wenn Sie cog nicht kennen, schauen Sie hier vorbei: https://www.cogeo.org/

• Überprüfen Sie, ob gdal_translate verfügbar ist

gdal_translate --version

Es sollte die GDAL-Version drucken, die Sie verwenden

• Raster auf 4326 umprojektieren

Ihr Raster hat möglicherweise andere Quell-SRS, ändern Sie es entsprechend

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

• Jetzt können wir TIF in Cloud-optimiertes Geotif konvertieren

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

Die Konvertierung der neu projizierten TIFF-Datei in Geotiff-Datei dauerte etwa eine Minute.

Fügen Sie OSM-Daten in die Postgresql-Tabelle ein

Wir verwenden osm2pgsql, um OSM-Daten in unsere Tabelle einzufügen

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

osm2pgsql benötigte insgesamt 274 Sekunden (4 Minuten und 34 Sekunden).

Sie können GeoJSON-Dateien auch verwenden, wenn Sie welche haben, indem Sie ogr2ogr verwenden

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

ogro2gr bietet eine breite Palette an Treiberunterstützung, sodass Sie bei Ihren Eingaben ziemlich flexibel sind. Ausgabe ist Postgresql-Tabelle

Berechnen Sie h3-Indizes

Postgresql

Installieren

pip install pgxnclient cmake
pgxn install h3

Erweiterung in Ihrer Datenbank erstellen

create extension h3;
create extension h3_postgis CASCADE;

Jetzt erstellen wir die Gebäudetabelle

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

Daten in Tabelle einfügen

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

Protokoll und 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

Jetzt berechnen wir die h3-Indizes für diese Gebäude mithilfe von centroid . Hier ist 8 die H3-Auflösung, an der ich arbeite. Erfahren Sie hier mehr über Auflösungen

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

Rasteroperationen

Installieren

pip install h3 h3ronpy rasterio asyncio asyncpg aiohttp

Stellen Sie sicher, dass das neu projizierte Zahnrad statisch ist/

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

Führen Sie das im Repository bereitgestellte Skript aus, um h3-Zellen aus Raster zu erstellen. Ich resample nach Modusmethode: Dies hängt von der Art der Daten ab, die Sie haben. Für kategoriale Daten ist der Modus besser geeignet. Erfahren Sie hier mehr über Resampling-Methoden

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

Protokoll :

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

Erstellen wir eine Funktion zum get_h3_indexes in einem Polygon

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;

Lassen Sie uns alle Gebäude abrufen, die in unserem Interessengebiet als bebaute Fläche identifiziert wurden.

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;

Abfrageplan:

Dies kann weiter verbessert werden, wenn ein Index zur h3_ix-Gebäudespalte hinzugefügt wird

create index on buildings(h3_ix);

Zählung der Aufnahme: In meiner Gegend gab es 24416 Gebäude mit der Bauklasse, die von ESRI klassifiziert wurde

Überprüfung

Lasst uns überprüfen, ob die Ausgabe wahr ist: Lasst uns die Gebäude als Geojson abrufen

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;

Lass uns auch h3-Zellen bekommen

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

Die Genauigkeit kann nach Erhöhung der h3-Auflösung erhöht werden und hängt auch von der Eingabe- und Resampling-Technik ab

Aufräumen

Entfernen Sie die Tische, die wir nicht brauchen

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

Optional – Vektorkacheln

Um die Kacheln zu visualisieren, können wir mit pg_tileserv schnell Vektorkacheln erstellen

• Laden Sie pg_tileserv herunter Laden Sie es über den oben angegebenen Link herunter oder verwenden Sie Docker
• Anmeldeinformationen exportieren

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

• Erteilen Sie unserer Tabellenauswahlberechtigung

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

• Erstellen wir Geometrie auf H3-Indizes zur Visualisierung (Sie können dies direkt aus der Abfrage heraus tun, wenn Sie Kacheln manuell aus st_asmvt erstellen)

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

• Erstellen Sie einen Index für dieses H3-Geom, um es effektiv zu visualisieren

CREATE INDEX idx_esri_landcover_geometry 
ON esri_landcover 
USING GIST (geometry);

• Laufen

  ./pg_tileserv

• Jetzt können Sie diese MVT-Kacheln basierend auf dem Kachelwert oder auf jede gewünschte Weise visualisieren, z. B. mit Maplibre! Sie können das verarbeitete Ergebnis auch visualisieren oder mit anderen Datensätzen kombinieren. Dies ist die Visualisierung, die ich für h3-Indizes basierend auf ihrem Zellwert in QGIS erstellt habe

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

Referenzen:

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