Як використовувати PostgreSQL для (військових) геоаналітичних задач

Підписуйтеся на Telegram-канал «DOU #tech», щоб не пропустити нові технічні статті

Привіт усім. Мене звати Тарас Кльоба, і я працюю в SoftServe на посаді Associate Director, Big Data & Analytics. Також я є співзасновником спільноти PostgreSQL Ukraine та military-tech волонтерського проєкту «Corvus Intelligence», який використовує технології для підвищення обороноздатності нашої країни, зокрема в розвідувальних операціях.

Також наша команда перемогла на Національному оборонному хакатоні, організованому РНБО, та в одному з челенджів NATO TIDE Hackathon 2023 у Варшаві. Моє залучення до різноманітних проєктів та ініціатив дозволили отримати практичний досвід у геоаналітиці, яким я поділюсь у цій публікації.

Геоаналітика є надзвичайно важливою у військовій справі, оскільки велика частина військових даних містить геоатрибути. У цій статті я розкажу про способи використання PostgreSQL для обробки геоданих і вирішення типових геоаналітичних завдань. Тут буде інформація про методи пошуку найближчих об’єктів, розрахунку відстаней і використанню геопросторових індексів для покращення цих процесів. Також розглянемо техніки визначення точки в полігоні та геопросторової агрегації. Метою статті є навести корисні приклади та поради, які можуть допомогти поліпшити роботу з геоданими та сприяти розвитку нових рішень.

Матеріали та дані у статті взяті з відкритих джерел та узгоджені з військовими.

Перше джерело даних: як імпортувати дані військових полігонів рф у PostgreSQL

Для початку аналізу та демонстрації можливостей PostgreSQL у геоаналітиці мені потрібні якісь набори даних. Я вирішив почати з даних про військові об’єкти рф, які доступні на OpenStreetMap (OSM). Першим кроком є завантаження цих даних до PostgreSQL, після чого ми зможемо використати інструменти для оптимізації запитів та підвищення їх ефективності.

Для імпорту даних про військові об’єкти рф з OpenStreetMap ми використаємо інструмент osm2pgsql. Цей open-source інструмент дозволяє ефективно переносити дані з OSM у PostgreSQL. Ми завантажимо файл russia-latest.osm.pbf (об’ємом 3,4 ГБ), який містить інформацію про точки, лінії, дороги та полігони з OSM. Після завантаження файл буде використаний для наповнення відповідних таблиць у PostgreSQL, де ми можемо почати аналіз та обробку даних.

Скрипт, який ми використовуємо, включає команди для завантаження OSM даних, створення нової бази даних PostgreSQL та імпорту даних за допомогою osm2pgsql:

#!/bin/bash
 
# Variable for the database and OSM file name
DB_NAME="russia_latest"  
 
# Construct the OSM data URL using the DB_NAME
OSM_DATA_URL="http://download.geofabrik.de/russia-latest.osm.pbf"
 
# Download the latest OSM data for Russia
echo "Downloading OSM data for Russia..."
wget $OSM_DATA_URL -O $DB_NAME.osm.pbf
 
# Create a new PostgreSQL database with the corrected name
echo "Creating new PostgreSQL database named $DB_NAME..."
sudo -u postgres psql -c "CREATE DATABASE $DB_NAME;"
sudo -u postgres psql -d $DB_NAME -c "CREATE EXTENSION postgis;"
 
# Import the OSM data into the PostgreSQL database using osm2pgsql
echo "Importing OSM data into PostgreSQL database named $DB_NAME..."
osm2pgsql --create \
          --verbose \
          --database $DB_NAME \
          --latlong \
          --number-processes 4 \
          --username postgres \
          --host localhost \
          --port 5432 \
          --password \
          $DB_NAME.osm.pbf
 
echo "Data import complete." 

Після виконання скрипта у нашій базі даних з’являться п’ять основних таблиць:

  • osm2pgsql_properties  зберігає налаштування та властивості, використані під час імпортування даних;
  • planet_osm_line  містить лінійні елементи, такі як дороги та річки;
  • planet_osm_point — включає точкові об’єкти, такі як будівлі (не всі будівлі позначені як географічні полігони, тому доведеться щось придумати для роботи із цими точками);
  • planet_osm_polygon — зберігає полігони, що представляють області, як-от військові бази;
  • planet_osm_roads — зберігає транспортні маршрути.

Для спрощення аналізу військових об’єктів створимо таблицю military_geometries. SQL-скрипт вибере дані з таблиць planet_osm_line, planet_osm_point, planet_osm_polygon та planet_osm_roads, фільтруючи військові об’єкти. Буферизація в 100 метрів застосовуватиметься до ліній, точок та доріг за допомогою ST_Buffer. Це дозволить також створити полігони на основі точок та ліній, забезпечуючи можливість аналізу, до прикладу, при розв’язанні задачі—чи входить точка у вказані полігони.

DROP TABLE IF EXISTS military_geometries;
 
CREATE TABLE military_geometries AS
SELECT osm_id, 'line' AS geom_type, landuse, military, building, name, operator,
       ST_Buffer(way, 0.0009)::geometry(Polygon, 4326) AS geom
FROM public.planet_osm_line
WHERE military IS NOT NULL OR building = 'military' OR landuse = 'military'
 
UNION ALL
 
SELECT osm_id, 'point' AS geom_type, landuse, military, building, name, operator,
       ST_Buffer(way, 0.0009)::geometry(Polygon, 4326) AS geom
FROM public.planet_osm_point
WHERE military IS NOT NULL OR building = 'military' OR landuse = 'military'
 
UNION ALL
 
SELECT osm_id, 'polygon' AS geom_type, landuse, military, building, name, operator,
       way::geometry(Polygon, 4326) AS geom
FROM public.planet_osm_polygon
WHERE military IS NOT NULL OR building = 'military' OR landuse = 'military'
 
UNION ALL
 
SELECT osm_id, 'road' AS geom_type, landuse, military, building, name, operator,
       ST_Buffer(way, 0.0009)::geometry(Polygon, 4326) AS geom
FROM public.planet_osm_roads
WHERE military IS NOT NULL OR building = 'military' OR landuse = 'military';
 
--
SELECT 9252 Query returned successfully in 12 secs 151 ms. 

Виконання наведеного SQL-скрипта дозволить нам сформувати таблицю military_geometries, яка міститиме полігони для 9,252 військових об’єктів, ідентифікованих на OSM:

9,252 військових об’єктів на росії та в тимчасово окупованому АР Крим, візуалізовано як точки за допомогою QGIS

В OSM, як і в інших відкритих джерелах, інформація може змінюватися. Наприклад, за період від початку 2022 року було видалено 2,995 військових об’єктів на росії.

Але, як то кажуть, «скриншоти не горять». Такі видалення часто призводять до ефекту Барбари Стрейзанд, коли спроби сховати інформацію лише привертають до неї ще більше уваги. Якщо ви хочете покопатися в історичних даних OSM і допомогти виявити такі аномалії, можна використовувати ресурси, як-от GeoFabrik.de. І хоча це не зовсім стосується нашого аналізу, хочу показати, як ці видалені об’єкти виглядають на карті, що відображує спроби росіян приховати важливі дані.

Видалені після 01.01.2022 (сині) та наявні (червоні) географічні полігони військових об’єктів у москві

Друге джерело даних: дані про пожежі із супутників NASA

Як наступне джерело даних ми використаємо інформацію з системи FIRMS (Fire Information for Resource Management System), розробленої у Мерілендському університеті за підтримки NASA та ООН у 2007 році. FIRMS забезпечує можливість спостерігати майже у режимі реального часу за активними пожежами по всьому світу, завдяки даним з супутників Aqua і Terra, оснащених спектрорадіометрами MODIS, а також VIIRS на супутниках S-NPP і NOAA 20. Інформація оновлюється кожні три години, а для США та Канади навіть частіше.

Ми ж використаємо дані FIRMS для ідентифікації пожеж на території військових об’єктів росії з 2022 року.

Щоб завантажити дані про пожежі з системи FIRMS, ми використаємо наступний скрипт, який витягує всі записи про пожежі на території росії з 1 січня 2022 року до сьогоднішнього дня. Ці дані потім імпортуються до нової таблиці viirs_fire_events в базі даних PostgreSQL.

from datetime import datetime, timedelta
import time
import requests
import psycopg2
import csv
from psycopg2.extras import execute_values
 
# Your NASA FIRMS API key
api_key = 'your_api_key'
 
# Start and end dates for the data retrieval
start_date = datetime(2022, 1, 1)
end_date = datetime(2023, 12, 24)
 
# URL for the NASA FIRMS API call
base_url = "https://firms.modaps.eosdis.nasa.gov/api/country/csv/{}/VIIRS_NOAA20_NRT/RUS/1/{}"
 
# Database connection details
db_name = 'your_db_name'
db_user = 'your_db_user'
db_password = 'your_db_password'
db_host = 'your_db_host'
db_port = 'your_db_port'
 
# Function to insert data into PostgreSQL
def insert_data_into_db(data, cursor):
    insert_query = """INSERT INTO viirs_fire_events (country_id, latitude, longitude, bright_ti4, scan, track, acq_date, acq_time, satellite, instrument, confidence, version, bright_ti5, frp, daynight) VALUES %s;"""
    execute_values(cursor, insert_query, data)
 
# Function to fetch and insert data into the database
def fetch_and_insert_data(date, cursor):
    formatted_date = date.strftime("%Y-%m-%d")
    try:
        response = requests.get(base_url.format(api_key, formatted_date))
        
        if "Invalid MAP_KEY" in response.text:
            print(f"Invalid MAP_KEY detected for {formatted_date}. Waiting for 10 minutes before retrying.")
            time.sleep(600)  # Wait for 10 minutes
            return fetch_and_insert_data(date, cursor)  # Retry
 
        decoded_content = response.content.decode('utf-8')
        csv_reader = csv.reader(decoded_content.splitlines(), delimiter=',')
        next(csv_reader, None)  # Skip the header
        rows = [tuple(row) for row in csv_reader if row]
 
        if rows:
            insert_data_into_db(rows, cursor)
            print(f"Data for {formatted_date} inserted successfully (rows = {len(rows)}).")
        else:
            print(f"No valid data to insert for {formatted_date}.")
 
    except Exception as e:
        # General exception handler for any other unexpected errors
        print(f"Error fetching data for {formatted_date}: {e}")
 
try:
    conn = psycopg2.connect(dbname=db_name, user=db_user, password=db_password, host=db_host, port=db_port)
    cur = conn.cursor()
 
    current_date = start_date
    while current_date <= end_date:
        fetch_and_insert_data(current_date, cur)
        conn.commit()
        current_date += timedelta(days=1)
    
except Exception as e:
    print(f"Database error: {e}")
finally:
    if cur is not None:
        cur.close()
    if conn is not None:
        conn.close()
 
print("Data fetching and insertion complete.") 

Таким способом ми наповнимо таблицю viirs_fire_events, яка буде містити 1,711,475 записів про пожежі на території росії, і фактично ці пожежі виглядають ось так:

Візуалізація пожеж на росії від 1 січня 2022 року (1,711,475 пожеж)

Таблиця viirs_fire_events у базі даних PostgreSQL буде використовуватися для зберігання деталізованих даних про пожежі, із полями для координат, супутникових параметрів, дати та часу зйомки, та інших критичних метаданих. Нова колонка з типом даних GEOMETRY(POINT, 4326) буде автоматично наповнюватись на основі даних із колонок longitude та latitude.

CREATE TABLE viirs_fire_events (
    country_id TEXT,
    latitude FLOAT,
    longitude FLOAT,
    bright_ti4 FLOAT,
    scan FLOAT,
    track FLOAT,
    acq_date DATE,
    acq_time INT,
    satellite TEXT,
    instrument TEXT,
    confidence TEXT,
    version TEXT,
    bright_ti5 FLOAT,
    frp FLOAT,
    daynight TEXT,
    geom GEOMETRY(POINT, 4326) GENERATED ALWAYS AS (ST_MakePoint(longitude, latitude)) STORED
); 

Якщо вам цікаво працювати з цими даними про військові об’єкти та пожежі, але описаний процес витягування наборів даних здається занадто довгим, я підготував для вас експортовані таблиці в форматі CSV. Ви можете завантажити їх за вказаними посиланнями: military_geometries, viirs_fire_events.

Пошук військових об’єктів, де були пожежі: точки у полігоні

Отже, фактично тепер у нас є дві таблиці military_geometries та viirs_fire_events. Спробуємо знайти ті військові об’єкти, на яких були пожежі (від початку 2022 року), або знайти ті, на яких ще не було 🙂.

Тож давайте використаємо SQL-запит з використанням функції ST_Contains, щоб знайти ті військові об’єкти, на яких були зафіксовані пожежі з супутників NASA:

SELECT mg.osm_id, 
    mg.geom_type, 
    mg.landuse, 
    mg.military, 
    mg.building, 
    mg.name, 
    mg.operator, 
    mg.geom
FROM public.military_geometries AS mg
WHERE EXISTS (
    SELECT 1
    FROM public.viirs_fire_events AS vfe
    WHERE ST_Contains(mg.geom, vfe.geom)
);
--
Successfully run. Total query runtime: 54 min 15 secs. 129 rows affected. 

Як ви могли зауважити, ми виявили 129 військових об’єктів, де були пожежі з початку 2022 року. Цікаво, що на деяких об’єктах пожежі могли виникати неодноразово:

Військові об’єкти, на яких відбувались пожежі від початку 2022 року. Прозорість об’єктів вказує на частоту виникання пожеж

Друге, що ви могли зауважити, це те, що зазначений запит виконувався 54 хвилини 15 секунд, що є досить довго для такої простої операції. Щоб зрозуміти причини такої тривалості, корисно застосувати команду EXPLAIN ANALYZE. Ця команда дозволить проаналізувати процес виконання запиту, виявити можливі «вузькі місця» та надалі оптимізувати запит для покращення продуктивності.

Nested Loop Semi Join  (cost=0.00..395827843592.52 rows=8 width=524) (actual time=4628.471..4005790.130rows=129 loops=1)
  Join Filter: st_contains(mg.geom, vfe.geom)
  Rows Removed by Join Filter: 15695214562
  ->  Seq Scan on military_geometries mg  (cost=0.00..649.52 rows=9252 width=524) (actual time=0.009..22.516rows=9252 loops=1)
  ->  Materialize  (cost=0.00..70285.12 rows=1711475 width=32) (actual time=0.015..121.310 rows=1696413loops=9252)
        ->  Seq Scan on viirs_fire_events vfe  (cost=0.00..50027.75 rows=1711475 width=32) (actual time=124.339..351.922 rows=1711475 loops=1)
Planning Time: 0.923 ms
Execution Time: 4005798.854 ms 

Як ви бачите, у цьому випадку, при використанні оператора Nested Loop Semi Join ми отримали задачу складності O(n*m), де n — це 9,252 рядків у таблиці military_geometries, а m — це 1,711,475 рядків у viirs_fire_events. Це означає, що кожен рядок з першої таблиці порівнюється з кожним рядком другої таблиці, що призводить до дуже великої кількості операцій.

Отже, давайте розберемось, як ми можемо пришвидшити виконання такого запиту за допомогою використання індексів.

Підвищення продуктивності: використання індексів у геоаналітиці

PostgreSQL відомий своїми можливостями по розширюваності: у цій базі даних є багато методів доступу до геоданих. І щоб відшукати всі методи, які підходять для роботи із точками у двовимірному просторі, ми можемо виконати такий запит:

SELECT am.amname AS access_method, typ.typname, opc.opcname AS operator_class 
FROM pg_am am
INNER JOIN pg_opclass opc ON am.oid = opc.opcmethod
INNER JOIN pg_type typ ON typ.oid = opc.opcintype
WHERE (opc.opcname LIKE '%geometry%' OR opc.opcname LIKE '%point%')
AND (
    opc.opcname NOT LIKE '%3d%' AND
    opc.opcname NOT LIKE '%4d%' AND
    opc.opcname NOT LIKE '%nd%'
); 

В результаті ми побачимо як мінімум п’ять методів доступу, включаючи btree, hash, gist, brin та spgist. Пропоную провести дослідження, створивши індекси для кожного з цих методів та класів операторів. Після створення індексів ми перевіримо швидкість виконання нашого запиту про пожежі на військових об’єктах росії, щоб з’ясувати, які методи є найефективнішими для нашої задачі.

Тип індексу

Клас оператора індексу

Оператор фільтрації

Час на створення індексу

Розмір індексу

Час виконання запиту

Коротке пояснення

btree

btree_geometry_ops

Відсутній відповідний оператор — індекс ігнорується для цього запиту

1 sec 918 ms

81 MB

53 min 45 sec (129 rows affected)

Підтримує запити на рівність та по діапазону; швидке та упорядковане отримання даних

hash

hash_geometry_ops

Відсутній відповідний оператор — індекс ігнорується для цього запиту

3 secs 158 ms

59 MB

53 min 15 sec (129 rows affected)

Швидкий пошук за рівністю; не підходить для упорядкування або діапазонних запитів

brin

brin_geometry_inclusion_ops_2d

@(geometry,geometry)

536 ms

0.032 MB

28 min 3 sec

(129 rows affected)

Ефективний для великих наборів даних з природно упорядкованими даними; індексує діапазони блоків, а не окремі рядки.

gist

gist_geometry_ops_2d

@(geometry,geometry)

11 secs 659 ms

94 MB

493 ms (129 rows affected)

Підтримує широкий спектр запитів, включаючи просторові пошуки на перекриття та близькість

spgist

spgist_geometry_ops_2d

@(geometry,geometry)

6 secs 290 ms

78 MB

353 ms

(129 rows affected)

Підходить для даних з нерівномірним розподілом; підтримує різноманітні розділені деревоподібні структури

Наступні класи оператора не використовують geometry типи даних для пошуку, а працюють із: <@(point,polygon) та <@(point,box) — тому у результатах, кількості рядків не співпадають (наприклад, складний географічний полігон був спрощений до прямокутника)

gist

point_ops

<@(point,polygon)

1 secs 426 ms

81 MB

306 ms (повернуло 132 записи)

Ідеально підходить для точкових даних; підтримує запити на просторові відносини, такі як міститься в межах та перетин

spgist

quad_point_ops

<@(point,box)

4 secs 849 ms

77 MB

243 ms (173 rows affected)

Використовує квадродерева для індексації точкових даних; ефективний у певних сценаріях просторового аналізу

spgist

kd_point_ops

<@(point,box)

5 secs 204 ms

93 MB

199 ms (173 rows affected)

Використовує kd-дерева для багатовимірних точкових даних; відмінно підходить для пошуку найближчих сусідів

В таблиці з результатами ми бачимо, що найефективнішими для нашої задачі є GiST та SP-GiST індекси. Розгляньмо, як вони працюють.

Як працює GiST

GiST (Generalized Search Tree) індекси у PostgreSQL дозволяють ефективно впорядковувати та виконувати пошук по різноманітних типах даних, використовуючи концепцію збалансованих дерев. Вони дають можливість розробляти власні оператори для індексації, що робить GiST доволі універсальним та адаптивним до специфічних потреб.

Ієрархічна структура GiST індексу в PostgreSQL [1]

На зображені приклад GiST дерева: на верхньому рівні розміщені R1 та R2, які є bounding boxes для інших елементів. R1 містить R3, R4 та R5, а R3 в свою чергу охоплює R8, R9 і R10. Індекс GiST, має ієрархічну структуру, внаслідок чого можна значно пришвидшувати пошук. На відміну від B-дерев, GiST підтримує операції перекриття та визначення просторових відносин. Саме тому GiST добре підходить для індексування геометричних даних.

Як працює SP-GiST

SP-GiST (Space Partitioning Generalized Search Tree) — індекси в PostgreSQL, призначені для структур даних, які розділяють простір на області, що не перетинаються, як-от дерева квадрантів або префіксні дерева. Вони дозволяють рекурсивно розділяти дані на підмножини, формуючи незбалансовані дерева. Це робить SP-GiST індекси особливо ефективними для використання в оперативній пам’яті, де вони можуть швидко обробляти запити за рахунок меншої кількості рівнів і малих груп даних у кожному вузлі.

Однак SP-GiST індекси мають недоліки при зберіганні на диску через високу кількість дискових операцій, потрібних для їх функціонування, особливо у великих базах даних.

Враховуючи це, GiST індекси часто є кращим вибором, зокрема при роботі з полігонами та складними просторовими структурами.

Пошук найближчих сусідів: 10 пожеж біля заводу з виробництва «Шахедів»

Тепер спробуймо вирішити задачу «пошуку найближчих сусідів» за допомогою PostgreSQL. Використовуючи наші набори даних, ми спробуємо виявити 10 пожеж, що виникли поряд із заводом на росії, де виробляють іранські дрони «Шахед». За детальнішою інформацією про завод можна звернутися до дослідження команди «Molfar». Завод розташований у спеціальній економічній зоні «Алабуга» в Татарстані, де раніше виробляли корм для котів, автоскло та вирощували гриби. Однак після санкцій проти росії його пріоритети змінилися, і тепер він має ключову роль у планах росії щодо виробництва дронів.

Одним із методів розв’язання цієї задачі є створення буфера у формі круга навколо обраної цілі. Цей буфер рекурсивно збільшується до отримання потрібної кількості результатів. У PostgreSQL це може бути реалізовано за допомогою ось такого SQL-запиту, який формує буфер та визначає пожежі, що відбулися в межах вказаного радіусу від обраного об’єкта:

SELECT *
FROM public.viirs_fire_events
WHERE ST_DWithin(
    geom, 
    ST_SetSRID(ST_MakePoint(52.048470, 55.821688), 4326)::geography, 
    10000 -- can be changed to 1000, 5000, 10000, ...
)
LIMIT 10;
 
--
Successfully run. Total query runtime: 1 secs 157 ms.
10 rows affected. 

Така методика включає поетапне збільшення буфера та аналіз результатів, що може зайняти багато часу.

Завод в Татарстані з виробництва «Шахедів», із радіусами у 1,5 та 10 км та візуалізація пожеж

Для оптимізації геопросторових запитів можна скористатися різноманітними операторами, які підтримуються GiST індексами. Щоб отримати перелік доступних операторів для використання із GiST індексом, можна виконати SQL-запит, який сканує системні таблиці PostgreSQL і повертає інформацію про оператори, асоційовані з класом операторів «gist_geometry_ops_2d». Це дозволить знайти найбільш ефективні оператори для виконання конкретних геопросторових операцій у базі даних.

SELECT
    amopopr::regoperator, 
    oprcode::regproc, 
    obj_description(opr.oid, 'pg_operator') description
FROM pg_am am
INNER JOIN pg_opclass opc ON opcmethod = am.oid
INNER JOIN pg_amop amop ON amopfamily = opcfamily
INNER JOIN pg_operator opr ON opr.oid = amopopr
WHERE amname = 'gist'
    AND opcname = 'gist_geometry_ops_2d'
ORDER BY amopstrategy;
 
--
 amopopr                |    oprcode            | description
------------------------+-----------------------+-----------------------------------
 <<(geometry,geometry)  | geometry_left         | Determines if one geometry is left of another
 &<(geometry,geometry)  | geometry_overleft     | Checks if one geometry overlaps the left of another
 &&(geometry,geometry)  | geometry_overlaps     | Checks if two geometries overlap
 &>(geometry,geometry)  | geometry_overright    | Determines if one geometry overlaps the right of another
 >>(geometry,geometry)  | geometry_right        | Determines if one geometry is right of another
 ~=(geometry,geometry)  | geometry_same         | Checks if two geometries are the same
 ~(geometry,geometry)   | geometry_contains     | Determines if one geometry contains another
 @(geometry,geometry)   | geometry_within       | Determines if one geometry is within another
 &<|(geometry,geometry) | geometry_overbelow    | Determines if one geometry overlaps below another
 <<|(geometry,geometry) | geometry_below        | Determines if one geometry is below another
 |>>(geometry,geometry) | geometry_above        | Determines if one geometry is above another
 |&>(geometry,geometry) | geometry_overabove    | Determines if one geometry overlaps above another
 <->(geometry,geometry) | geometry_distance_centroid | Calculates centroid distance between geometries
 <#>(geometry,geometry) | geometry_distance_box | Calculates box distance between geometries 

Наш GiST індекс надає широкі можливості для роботи з геоданими, дозволяючи визначати просторове розташування об’єктів та вимірювати відстані. Використання оператора «<->» дозволяє впорядкувати об’єкти за близькістю до заданої точки. У цьому прикладі ми використовуємо цей оператор для визначення десяти найближчих пожеж до вказаної локації.

SELECT *
FROM public.viirs_fire_events
ORDER BY geom <-> ST_SetSRID(ST_MakePoint(52.048470, 55.821688), 4326)
LIMIT 10;
 
--
Successfully run. Total query runtime: 78 ms.
10 rows affected. 

Запит виявився значно швидшим — у 15 разів, порівняно з попередньою методикою, і це без повторних виконань зі зміненим радіусом. Щоб підтвердити, що швидкість зросла завдяки використанню індексу та оператора, ми можемо аналізувати план запиту. Так ми переконаємось, що індекс дійсно був залучений, а це і є ключем до збільшення продуктивності.

EXPLAIN ANALYZE
SELECT *
FROM public.viirs_fire_events
ORDER BY geom <-> ST_SetSRID(ST_MakePoint(52.048470, 55.821688), 4326)
LIMIT 10;
 
--
QUERY PLAN
Limit  (cost=0.41..14.22 rows=10 width=143) (actual time=0.230..0.254 rows=10 loops=1)
  ->  Index Scan using idx_viirs_fire_events_geom_gist on viirs_fire_events  (cost=0.41..2363046.98rows=1711475 width=143) (actual time=0.229..0.252 rows=10 loops=1)
        Order By: (geom <-> '0101000020E6100000276BD44334064A4001C287122DE94B40'::geometry)
Planning Time: 0.105 ms
Execution Time: 0.283 ms 

Як бачимо, індекси, подібні до GiST, розширюють аналітичні можливості за межі простих порівнянь, дозволяючи вирішувати складніші завдання. І як показано у цій статті, відкриті дані можна ефективно використовувати для швидкої оцінки та визначення цілей на глобальному рівні, у тому числі для оцінювання успішності ураження цілей.

H3 від Uber: один із підходів до геоаналітики та агрегації даних

H3, розроблений компанією Uber, є системою гексагональної (шестикутної) сітки, що забезпечує гнучке та ефективне розподілення геопросторових даних. Мені здається, що H3 має всі шанси стати спільним стандартом для роботи із геоданими у ЗСУ. Ми дослідимо, як цей інструмент можна використовувати для агрегації даних і вирішення складних геоаналітичних задач.

Ілюстрація гексагональної сітки H3 від Uber

Як видно на зображенні, кожен шестикутник є окремою географічною одиницею, що дозволяє спрощувати обробку складних геоформ до рівномірних сегментів.

Рівень

Загальна кількість обʼєктів

Кількість гексагонів

Кількість пентагонів

0

122

110

12

1

842

830

12

2

5882

5870

12

3

41 162

41 150

12

4

288 122

288 110

12

5

2 016 842

2 016 830

12

6

14 117 882

14 117 870

12

7

98 825 162

98 825 150

12

8

691 776 122

691 776 110

12

9

4 842 432 842

4 842 432 830

12

10

33 897 029 882

33 897 029 870

12

11

237 279 209 162

237 279 209 150

12

12

1 660 954 464 122

1 660 954 464 110

12

13

11 626 681 248 842

11 626 681 248 830

12

14

81 386 768 741 882

81 386 768 741 870

12

15

569 707 381 193 162

569 707 381 193 150

12

Це ієрархічна система, що включає 15 рівнів поділу земної поверхні на гексагони (шестикутники), нульовий рівень, поділений на 122 секції, з яких 12 є п’ятикутниками для точного відображення сферичної форми Землі. На найдрібнішому рівні ми маємо близько 569 трильйонів гексагонів, кожен із яких представляє окремий геопросторовий об’єкт. У прикладі, на відео можна побачити, як це працює на практиці.

PostgreSQL може інтегрувати функціонал H3 через додаткове розширення. Це розширення встановлюється командою CREATE EXTENSION h3; і доступне у сервісах хмарних обчислень, включно з AWS RDS (принагідно дякую AWS за їхню підтримку України). І разом із встановленням цього розширення у вас зʼявляються нові функції. Розглянемо ті, які можуть знадобитись у роботі початківців:

Функція

Вхідні дані

Вихідні дані

Опис

h3_lat_lng_to_cell

latitude: FLOAT, longitude: FLOAT, resolution: INT

H3 index: BIGINT

Перетворює координати широти та довготи в H3 індекс на заданому рівні розділення.

h3_cell_to_boundary

H3 index: BIGINT

Array of boundary coordinates: GEOMETRY(POLYGON, 4326)

Перетворює H3 індекс у геометричний полігон, що представляє границі гексагона.

h3_get_resolution

H3 index: BIGINT

Resolution level: INT

Повертає рівень розділення заданого H3 індексу.

h3_cell_to_parent

H3 index: BIGINT, desired resolution: INT

Parent H3 index: BIGINT

Перетворює H3 індекс на батьківський індекс на вищому рівні ієрархії.

h3_cell_to_children

H3 index: BIGINT, desired resolution: INT

Array of child H3 indexes: SETOF BIGINT

Перетворює H3 індекс на масив дочірніх індексів на нижчому рівні ієрархії.

h3_polygon_to_cells

geometry: GEOMETRY, resolution: INT

Array of H3 indexes: SETOF BIGINT

Перетворює полігон на множину H3 індексів, що повністю або частково покривають полігон.

h3_grid_disk

H3 index: BIGINT, range: INT

Array of H3 indexes: SETOF BIGINT

Генерує масив H3 індексів, що представляють сітку гексагонів (grid), навколо центрального H3 індексу, формуючи «диск» визначеного радіуса.

h3_compact_cells

Array of H3 indexes: SETOF BIGINT

Array of compact H3 indexes: SETOF BIGINT

Консолідує масив H3 індексів, зменшуючи кількість індексів, які покривають ту саму область.

Для ефективного вирішення першої задачі, ми можемо перетворити всі наші полігони на масиви H3 індексів (гексагонів) заданого рівня, а також аналогічно обробити центроїди пожеж, перетворюючи їх у H3 індекси. Отримавши тип даних BIGINT для цих H3 індексів, ми можемо застосувати звичайний індекс B-дерево, який є особливо ефективним у виконанні операцій порівняння на рівність. Це значно підвищить швидкість виконання запитів у складних задачах геоаналізу, забезпечуючи швидкі та точні результати.

Розглянемо декілька простих функцій H3, які допоможуть краще зрозуміти, як це працює на практиці:

h3_polygon_to_cells(geom, 8)

h3_grid_disk(h3_polygon_to_cells(geom, 8), 1)

h3_polygon_to_cells(geom, 9)

Ця функція конвертує геометрію військового полігону у набір H3 індексів на 8-му рівні розділення, який ефективно ділить полігон на шестикутники, кожен з яких покриває площу в 0.737327598 квадратних кілометрів, забезпечуючи детальний просторовий аналіз.

Якщо деякі зони полігону залишилися непокритими після застосування h3_polygon_to_cells, можна використати h3_grid_disk для створення додаткового кільця H3 індексів, що розширює покриття, додаючи шестикутники навколо наявних індексів та забезпечуючи повне покриття визначного географічного полігону.

Використання функції h3_polygon_to_cells з рівнем 9 підвищує роздільну здатність сітки до більш дрібної шкали, де кожен шестикутник відображає площу в 0.105332513 квадратних кілометрів. Це дозволяє досягти більшої точності у відтворенні геометрії географічного полігону для деталізованого просторового аналізу, але разом із тим ми отримуємо більше число гексагонів, що може негативно вплинути на подальшу швидкодію виконання запитів.

Під час мого виступу на PGConf.2023 — найбільшій конференції в Європі, присвяченій PostgreSQL — я мав нагоду презентувати ряд складніших задач, які можна вирішити через агрегацію геоданих із застосуванням H3. Одним із прикладів був пошук інших дронів, що знаходились у тій самій локації та часі, а також аналіз маршрутів дронів, які подорожували разом, ідентифікованих у різних місцях протягом певного періоду часу (Companion Analysis). Детальніше ознайомитися з цією темою можна буде у продовженні цієї статті. До цього часу, ви можете переглянути матеріали моєї презентації за посиланням.

У межах наших наборів даних ми можемо проаналізувати військові обʼєкти і за допомогою агрегації у H3 порахувати густоту військових обʼєктів на росії. Ця візуалізація виглядає так:

Візуалізація густоти військових обʼєктів на території росії та у тимчасово окупованому АР Крим за допомогою H3 шестикутників.

Використання H3 для агрегації геоданих значно покращує аналітичні можливості, дозволяючи глибше інтерпретувати та візуалізувати складні просторові відносини.

На завершення

Якщо ви є представником Збройних Сил України і шукаєте кваліфіковану безкоштовну підтримку у сфері даних, Big Data або геоаналітики, запрошую у приватні повідомлення, я і моя команда волонтерів з радістю підтримаємо вас своїми знаннями та ресурсами.

Додаткові ресурси, які я використав для підготовки цієї статті

  1. Mastering PostgreSQL by Hans-Jürgen Schönig
  2. Inside the russian effort to build 6,000 attack drones with Iran’s help
  3. Uber H3. Tables of Cell Statistics Across Resolutions
  4. H3-PG Extension. API Reference
  5. PostGIS documentation
  6. PostGIS in Action, Third Edition by Leo S. Hsu and Regina Obe

P.S. Для підготовки цієї статті я використовував Ubuntu сервер із такими характеристиками та налаштування бази даних PostgreSQL:


# DB Version: 16
# OS Type: linux
# DB Type: dw
# Total Memory (RAM): 8 GB
# CPUs num: 4
# Connections num: 20
# Data Storage: ssd

max_connections = 20

shared_buffers = 2 GB
effective_cache_size = 6 GB
maintenance_work_mem = 1 GB
checkpoint_completion_target = 0.9
wal_buffers = 16 MB
default_statistics_target = 500
random_page_cost = 1.1
effective_io_concurrency = 200
work_mem = 26214 kB
huge_pages = off
min_wal_size = 4 GB
max_wal_size = 16 GB
max_worker_processes = 4
max_parallel_workers_per_gather = 2
max_parallel_workers = 4
max_parallel_maintenance_workers = 2

# Recommendations are generated by pgtune.leopard.in.ua

👍ПодобаєтьсяСподобалось60
До обраногоВ обраному25
LinkedIn
Дозволені теги: blockquote, a, pre, code, ul, ol, li, b, i, del.
Ctrl + Enter
Дозволені теги: blockquote, a, pre, code, ul, ol, li, b, i, del.
Ctrl + Enter

чудова робота :)

Мєм ще в тому, що багато росіян докладало і докладає руки до розробки постгресу. Ті ж постгрес про. xd

Мєм ще в тому, що багато росіян докладало і докладає руки до розробки постгресу. Ті ж постгрес про. xd

Крупній открітій проект, кого там только нету.

Не самого постгреса, а його форка

Крута технічна стаття. Дякую!

Невже на DOU повернулися технічні статті? 😂
Тарасе, дякую. Приємно було читати

як запустити скрипти підкажіть ?

встановіть osm2pgsql
Тоді запустіть баш скрипт на початку статті. А все наступне — то звичайний SQL, запускайте де вам зручно — psql, pgadmin, dbeaver, datagrip etc.
Інфо про пожежі — то Python, відповідно треба встановити Python)

Так, я вже зрозумів, Тарас підказав... Я працюю у Windows. Зроблю через термінал. а потім через pgadmin останні дії. Головне щоб SQL запит спрацював.
Дякю за приділену увагу !

Отличная работа. Кстати в своё время устав от GeoMesa перешёл на прототип на базе PostGIS + Tile38 — интересный зверь, рекомендую

Дякую, цікаво. Все ніяк не вдається попрацювати на проекті\продукті з геоданими

Підписатись на коментарі