مرحبًا، سنتحدث في هذه المدونة عن كيفية إجراء تحليل البيانات النقطية بسهولة باستخدام فهارس h3.

موضوعي

للتعلم، سنقوم بإجراء العمليات الحسابية لمعرفة عدد المباني الموجودة في منطقة المستوطنات التي تحددها ESRI Land Cover. دعونا نستهدف البيانات على المستوى الوطني لكل من المتجهات والنقطية.

دعونا أولا العثور على البيانات

تنزيل البيانات النقطية

لقد قمت بتنزيل منطقة التسوية من Esri Land Cover .

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

لنقوم بتنزيل العام 2023 بحجم 362 ميجا تقريبا

تحميل مباني OSM في نيبال

المصدر: http://download.geofabrik.de/asia/nepal.html

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

المعالجة المسبقة للبيانات

لنطبق بعض المعالجة المسبقة على البيانات قبل إجراء حسابات خلايا h3 الفعلية
سوف نستخدم برنامج سطر الأوامر gdal لهذه الخطوة. تثبيت gdal في جهازك

التحويل إلى Geotiff المحسنة

إذا لم تكن على علم بالعتاد، قم بالخروج هنا: https://www.cogeo.org/

• تحقق من توفر gdal_translate

gdal_translate --version

يجب طباعة إصدار gdal الذي تستخدمه

• إعادة عرض البيانات النقطية إلى 4326

قد تحتوي بياناتك النقطية على مصدر مختلف، قم بتغييره وفقًا لذلك

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

• الآن يتيح تحويل tif إلى Geotif المحسّن على السحابة

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

استغرق الأمر دقيقة تقريبًا لتحويل tiff المعاد إسقاطه إلى Geotiff

أدخل بيانات osm في جدول postgresql

نحن نستخدم osm2pgsql لإدراج بيانات osm في جدولنا

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

استغرق osm2pgsql 274 ثانية (4 دقائق و34 ثانية) بشكل عام.

يمكنك استخدام ملفات Geojson أيضًا إذا كان لديك أي ملفات تستخدم ogr2ogr

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

يتمتع ogro2gr بمجموعة واسعة من الدعم لبرامج التشغيل، لذا فأنت مرن جدًا بشأن مدخلاتك. الإخراج هو جدول Postgresql

حساب فهارس h3

Postgresql

ثَبَّتَ

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 مبنى في منطقتي مع فئة البناء المصنفة من ESRI

تَحَقّق

لنتحقق مما إذا كان الإخراج صحيحًا: لنحصل على المباني كـ 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;

لنحصل على خلايا 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;

• قم بإنشاء فهرس على تلك Geom h3 لتصورها بشكل فعال

CREATE INDEX idx_esri_landcover_geometry 
ON esri_landcover 
USING GIST (geometry);

• يجري

  ./pg_tileserv

• الآن يمكنك تصور مربعات MVT تلك بناءً على قيمة المربعات أو بأي طريقة تريدها، على سبيل المثال: Maplibre! يمكنك أيضًا تصور النتيجة المعالجة أو دمجها مع مجموعات البيانات الأخرى. هذا هو التصور الذي قمت به على فهارس h3 استنادًا إلى قيمة الخلية الخاصة بها في QGIS

المصدر الريبو: https://github.com/kshitijrajsharma/raster-ana Analysis-using-h3

مراجع :

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