नमस्कार, इस ब्लॉग में हम इस बारे में बात करेंगे कि हम h3 इंडेक्स का उपयोग करके आसानी से रैस्टर विश्लेषण कैसे कर सकते हैं।

उद्देश्य

सीखने के लिए, हम यह पता लगाने पर गणना करेंगे कि ईएसआरआई लैंड कवर द्वारा निर्धारित निपटान क्षेत्र में कितनी इमारतें हैं। आइए वेक्टर और रैस्टर दोनों के लिए राष्ट्रीय स्तर के डेटा का लक्ष्य रखें।

आइए सबसे पहले डेटा ढूंढें

रैस्टर डेटा डाउनलोड करें

मैंने ईएसआरआई लैंड कवर से निपटान क्षेत्र डाउनलोड किया है।

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

आइए 2023 वर्ष डाउनलोड करें, आकार लगभग 362एमबी

नेपाल की ओएसएम बिल्डिंग डाउनलोड करें

स्रोत: http://download.geofabrik.de/asia/nepal.html

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

डेटा को प्रीप्रोसेस करें

आइए वास्तविक h3 सेल गणना से पहले डेटा पर कुछ प्रीप्रोसेसिंग लागू करें
हम इस चरण के लिए जीडीएएल कमांडलाइन प्रोग्राम का उपयोग करेंगे। अपनी मशीन में गदल स्थापित करें

क्लाउड अनुकूलित जियोटीफ़ में रूपांतरण

यदि आप कोग से अनजान हैं, तो यहां चेकआउट करें: https://www.cogeo.org/

• जांचें कि क्या gdal_translate उपलब्ध है

gdal_translate --version

इसे आपके द्वारा उपयोग किए जा रहे जीडीएल संस्करण को प्रिंट करना चाहिए

• रैस्टर को 4326 पर पुन: प्रोजेक्ट करें

आपके रैस्टर का स्रोत भिन्न हो सकता है, इसे तदनुसार बदलें

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

• अब टीआईएफ को क्लाउड ऑप्टिमाइज्ड जियोटीफ में परिवर्तित करते हैं

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

पुन: प्रक्षेपित टिफ़ को जियोटीफ़ में बदलने में लगभग एक मिनट का समय लगा

Postgresql तालिका में osm डेटा डालें

हम अपनी तालिका में ओएसएम डेटा डालने के लिए osm2pgsql का उपयोग कर रहे हैं

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

osm2pgsql ने कुल मिलाकर 274 सेकेंड (4 मिनट 34 सेकेंड) का समय लिया।

यदि आपके पास ogr2ogr का उपयोग करने वाली कोई फ़ाइल है तो आप जियोजोन फ़ाइलों का भी उपयोग कर सकते हैं

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

ogro2gr के पास ड्राइवरों के लिए व्यापक समर्थन है, इसलिए आप अपने इनपुट के बारे में काफी लचीले हैं। आउटपुट Postgresql तालिका है

h3 सूचकांक की गणना करें

पोस्टग्रेस्क्ल

स्थापित करना

pip install pgxnclient cmake
pgxn install h3

अपने डेटाबेस में एक्सटेंशन बनाएं

create extension h3;
create extension h3_postgis CASCADE;

अब भवन तालिका बनाएं

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

तालिका में डेटा डालें

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

लॉग और समय :

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

अब सेंट्रोइड का उपयोग करके उन इमारतों के लिए h3 इंडेक्स की गणना करते हैं। यहां 8 h3 रिज़ॉल्यूशन है जिस पर मैं काम कर रहा हूं। यहां संकल्पों के बारे में और जानें

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

रेखापुंज संचालन

स्थापित करना

pip install h3 h3ronpy rasterio asyncio asyncpg aiohttp

सुनिश्चित करें कि पुन: प्रक्षेपित कोग स्थिर है/

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

रैस्टर से h3 सेल बनाने के लिए रेपो में दी गई स्क्रिप्ट चलाएँ। मैं मोड विधि द्वारा पुन: नमूनाकरण कर रहा हूं: यह आपके पास मौजूद डेटा के प्रकार पर निर्भर करता है। श्रेणीबद्ध डेटा के लिए मोड बेहतर फिट बैठता है। पुन: नमूनाकरण विधियों के बारे में यहां अधिक जानें

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

लकड़ी का लट्ठा :

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

विश्लेषण

आइए बहुभुज में get_h3_indexes के लिए एक फ़ंक्शन बनाएं

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;

आइए उन सभी इमारतों को प्राप्त करें जिन्हें हमारे रुचि के क्षेत्र में निर्मित क्षेत्र के रूप में पहचाना गया है

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;

क्वेरी योजना :

इमारतों के h3_ix कॉलम पर इंडेक्स जोड़ने पर इसे और बढ़ाया जा सकता है

create index on buildings(h3_ix);

शूटिंग के समय गिनती: मेरे क्षेत्र में 24416 इमारतें थीं, जिनका निर्माण ईएसआरआई के अनुसार वर्गीकृत किया गया था

सत्यापन

आइए सत्यापित करें कि क्या आउटपुट सही है: आइए इमारतों को जियोजोन के रूप में प्राप्त करें

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;

आइए 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

H3 रिज़ॉल्यूशन बढ़ाने के बाद सटीकता बढ़ाई जा सकती है और यह इनपुट और पुनः नमूनाकरण तकनीक पर भी निर्भर करेगा

साफ - सफाई

उन तालिकाओं को हटा दें जिनकी हमें आवश्यकता नहीं है

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

वैकल्पिक - वेक्टर टाइलें

टाइल्स को विज़ुअलाइज़ करने के लिए pg_tileserv का उपयोग करके शीघ्रता से वेक्टर टाइल्स बनाएं

• pg_tileserv डाउनलोड करें ऊपर दिए गए लिंक से डाउनलोड करें या डॉकर का उपयोग करें
• निर्यात क्रेडेंशियल

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

• हमारी तालिका चयन की अनुमति दें

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

• आइए विज़ुअलाइज़ेशन के लिए h3 इंडेक्स पर ज्यामिति बनाएं (यदि आप मैन्युअल रूप से st_asmvt से टाइल्स बना रहे हैं तो आप इसे सीधे क्वेरी से कर सकते हैं)

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

• उस h3 जियोम को प्रभावी ढंग से देखने के लिए उस पर इंडेक्स बनाएं

CREATE INDEX idx_esri_landcover_geometry 
ON esri_landcover 
USING GIST (geometry);

• दौड़ना

  ./pg_tileserv

• अब आप उन एमवीटी टाइलों को टाइल्स के मूल्य के आधार पर या किसी भी तरह से देख सकते हैं, उदाहरण के लिए: मैपलिब्रे! आप संसाधित परिणाम की कल्पना भी कर सकते हैं या अन्य डेटासेट के साथ संयोजन कर सकते हैं। यह वह विज़ुअलाइज़ेशन है जो मैंने QGIS में उनके सेल_वैल्यू के आधार पर h3 इंडेक्स पर किया था

स्रोत रेपो: https://github.com/क्षितिजराजशर्मा/रास्टर-एनालिसिस-यूजिंग-एच3

सन्दर्भ :

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