Analisis Raster Menggunakan Uber hndexes dan PostgreSQL-Tutorial Python-php.cn

Hai , Dalam blog ini kita akan bercakap tentang bagaimana kita boleh melakukan analisis Raster dengan mudah menggunakan indeks h3.

Objektif

Untuk pembelajaran, Kami akan membuat pengiraan untuk mengetahui berapa banyak bangunan yang ada di kawasan penempatan yang ditentukan oleh ESRI Land Cover. Mari kita sasarkan data peringkat kebangsaan untuk kedua-dua vektor dan raster.

Mari kita cari data dahulu

Muat turun Data Raster

Saya telah memuat turun kawasan penempatan daripada Esri Land Cover .

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

Mari muat turun tahun 2023, bersaiz lebih kurang 362MB

Raster Analysis Using Uber hndexes and PostgreSQL

Muat turun Bangunan OSM Nepal

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

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

Salin selepas log masuk

Praproses data

Mari gunakan beberapa prapemprosesan pada data sebelum pengiraan sel h3 sebenar
Kami akan menggunakan program baris perintah gdal untuk langkah ini. Pasang gdal dalam mesin anda

Penukaran kepada Geotif Dioptimumkan Awan

Jika anda tidak mengetahui tentang cog , Semak di sini : https://www.cogeo.org/

Semak sama ada gdal_translate tersedia

gdal_translate --version

Salin selepas log masuk

Ia sepatutnya mencetak versi gdal yang anda gunakan

Unjurkan semula raster ke 4326

Raster anda mungkin mempunyai srs sumber yang berbeza, ubahnya dengan sewajarnya

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

Salin selepas log masuk

Sekarang mari tukar tif kepada geotif yang dioptimumkan awan

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

Salin selepas log masuk

Ia mengambil masa lebih kurang seminit untuk menukar tiff yang diunjurkan semula kepada geotif

Masukkan data osm ke jadual postgresql

Kami menggunakan osm2pgsql untuk memasukkan data osm ke jadual kami

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

Salin selepas log masuk

osm2pgsql mengambil masa 274s (4m 34s) secara keseluruhan.

Anda boleh menggunakan fail geojson juga jika anda ada menggunakan ogr2ogr

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

Salin selepas log masuk

ogro2gr mempunyai pelbagai sokongan untuk pemandu supaya anda agak fleksibel tentang input anda. Output ialah jadual Postgresql

Kira indeks h3

Postgresql

Pasang

pip install pgxnclient cmake pgxn install h3

Salin selepas log masuk

Buat sambungan dalam db anda

create extension h3; create extension h3_postgis CASCADE;

Salin selepas log masuk

Sekarang mari buat jadual bangunan

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

Salin selepas log masuk

Masukkan data ke jadual

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

Salin selepas log masuk

Log dan masa :

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

Salin selepas log masuk

Sekarang mari mengira indeks h3 untuk bangunan tersebut menggunakan centroid. Di sini 8 ialah resolusi h3 yang sedang saya usahakan. Ketahui lebih lanjut tentang resolusi di sini

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

Salin selepas log masuk

Operasi raster

Pasang

pip install h3 h3ronpy rasterio asyncio asyncpg aiohttp

Salin selepas log masuk

Pastikan cog yang diunjurkan semula berada dalam keadaan statik/

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

Salin selepas log masuk

Jalankan skrip yang disediakan dalam repo untuk mencipta sel h3 daripada raster . Saya sampling semula mengikut kaedah mod : Ini bergantung pada jenis data yang anda ada . Untuk mod data kategori lebih sesuai. Ketahui lebih lanjut tentang kaedah pensampelan semula di sini

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

Salin selepas log masuk

Log :

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

Salin selepas log masuk

Analisis

Mari buat fungsi untuk mendapatkan_h3_indexes dalam poligon

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;

Salin selepas log masuk

Mari dapatkan semua bangunan yang dikenal pasti sebagai kawasan binaan di kawasan minat kami

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;

Salin selepas log masuk

Pelan Pertanyaan :

Raster Analysis Using Uber hndexes and PostgreSQL

Ini boleh dipertingkatkan lagi jika ditambah indeks pada lajur h3_ix bangunan

create index on buildings(h3_ix);

Salin selepas log masuk

Apabila kiraan menembak : terdapat 24416 bangunan di kawasan saya dengan kelas binaan diklasifikasikan sebagai dari ESRI

Pengesahan

Mari sahkan sama ada output adalah benar : Mari dapatkan bangunan sebagai 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;

Salin selepas log masuk

Jom dapatkan sel h3 juga

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

Salin selepas log masuk

Raster Analysis Using Uber hndexes and PostgreSQL

Ketepatan boleh ditingkatkan selepas meningkatkan resolusi h3 dan juga bergantung pada input dan teknik pensampelan semula

Pembersihan

Jatuhkan meja yang kita tidak perlukan

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

Salin selepas log masuk

Pilihan - Jubin Vektor

Untuk menggambarkan jubin, mari bina jubin vektor dengan cepat menggunakan pg_tileserv

Muat turun pg_tileserv Muat turun dari pautan yang disediakan di atas atau gunakan docker
Eksport kelayakan

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

Salin selepas log masuk

Berikan kebenaran memilih meja kami

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

Salin selepas log masuk

Mari buat geometri pada indeks h3 untuk visualisasi (Anda boleh melakukan ini terus dari pertanyaan jika anda membina jubin dari st_asmvt secara manual)

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

Salin selepas log masuk

Buat indeks pada geom h3 itu untuk menggambarkannya dengan berkesan

CREATE INDEX idx_esri_landcover_geometry ON esri_landcover USING GIST (geometry);

Salin selepas log masuk

Lari

./pg_tileserv

Salin selepas log masuk

Raster Analysis Using Uber hndexes and PostgreSQL

Kini anda boleh memvisualisasikan jubin MVT tersebut berdasarkan nilai jubin atau apa-apa cara yang anda mahu cth : maplibre ! Anda juga boleh menggambarkan hasil yang diproses atau digabungkan dengan set data lain. Ini adalah visualisasi yang saya lakukan pada indeks h3 berdasarkan cell_value mereka dalam QGIS

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

Rujukan:

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

Atas ialah kandungan terperinci Analisis Raster Menggunakan Uber hndexes dan PostgreSQL. Untuk maklumat lanjut, sila ikut artikel berkaitan lain di laman web China PHP!