Анализ геометрических полигонов для оценки территориальных изменений

В данной статье рассматривается реализация алгоритма для анализа временны́х изменений геометрических полигонов. Код загружает актуальные геоданные, вычисляет площади территорий на разные даты и визуализирует результаты, демонстрируя возможности Engee для задач обработки геопространственной информации.

В данном примере используются следующие пакеты:

Downloads: Для простой загрузки файлов по URL.
JSON3: Высокопроизводительный пакет для парсинга JSON-данных.
CodecZlib: Обеспечивает работу с сжатыми данными, в частности, распаковку GZIP.
Dates: Работа с датами для фильтрации и форматирования.
Plots: Мощная и гибкая библиотека для построения графиков и визуализации данных, используется для отрисовки карты.
GeoJSON (хоть и добавлен, но не используется напрямую в коде): Пакет для работы с геоданными. В данном скрипте парсинг осуществляется через JSON3, что также является валидным подходом.

Pkg.activate("MyGeoProject", shared=true)
Pkg.add(["JSON3", "CodecZlib", "GeoJSON"])
using Downloads, JSON3, CodecZlib, Dates, GeoJSON

Код начинается с определения констант:

DATE_START — фиксированная дата начала анализа.
DATE_CURRENT — текущая дата, вычисляемая в момент запуска скрипта. Это делает анализ всегда актуальным.
DSM_URL — URL-адрес, указывающий на источник данных.
REGIONS_DATA - содержит справочную информацию по областям Украины: название, площадь и координаты для подписи на карте.
UKRAINE_TOTAL_AREA — общая площадь страны.

DATE_START = Date("2022-02-24")
DATE_CURRENT = Date(Dates.now())

2025-09-23

DSM_URL = "https://raw.githubusercontent.com/cyterat/deepstate-map-data/main/deepstate-map-data.geojson.gz"

"https://raw.githubusercontent.com/cyterat/deepstate-map-data/main/deepstate-map-data.geojson.gz"

REGIONS_DATA = [
    ("Винницкая", 26513, [28.47, 49.23]),
    ("Волынская", 20144, [25.33, 50.75]),
    ("Днепропетровская", 31914, [35.02, 48.45]),
    ("Донецкая", 26517, [37.80, 48.00]),
    ("Житомирская", 29832, [28.65, 50.25]),
    ("Закарпатская", 12777, [22.29, 48.62]),
    ("Запорожская", 27183, [35.18, 47.84]),
    ("Ивано-Франковская", 13928, [24.71, 48.92]),
    ("Киевская", 28131, [30.52, 50.45]),
    ("Кировоградская", 24588, [32.26, 48.51]),
    ("Луганская", 26684, [39.30, 48.57]),
    ("Львовская", 21833, [24.03, 49.84]),
    ("Николаевская", 24598, [31.98, 46.97]),
    ("Одесская", 33310, [30.73, 46.48]),
    ("Полтавская", 28748, [34.55, 49.59]),
    ("Ровненская", 20047, [26.25, 50.62]),
    ("Сумская", 23834, [34.80, 50.91]),
    ("Тернопольская", 13823, [25.60, 49.55]),
    ("Харьковская", 31415, [36.25, 50.00]),
    ("Херсонская", 28461, [32.62, 46.63]),
    ("Хмельницкая", 20629, [27.00, 49.42]),
    ("Черкасская", 20900, [32.06, 49.44]),
    ("Черниговская", 31865, [31.30, 51.50]),
    ("Черновицкая", 8097, [25.94, 48.29]),
    ("Крым", 26081, [34.10, 45.35]),
    ("Севастополь", 864, [33.52, 44.61])
]
UKRAINE_TOTAL_AREA = 603628

Функция get_gzip_data:

Downloads.download(url) загружает файл по URL во временный файл и возвращает путь к нему.
open(temp_file) do file открывает этот временный файл.
GzipDecompressorStream(file) создает поток декомпрессии, который "оборачивает" файл. При чтении из этого потока данные автоматически распаковываются.
read(stream, String) читает весь распакованный поток в строку, которая и возвращается функцией.

Это решение для работы с сжатыми данными прямо в памяти, без необходимости управлять временными файлами вручную.

function get_gzip_data(url::String)
    temp_file = Downloads.download(url)
    open(temp_file) do file
        stream = GzipDecompressorStream(file)
        return read(stream, String)
    end
end

get_gzip_data (generic function with 1 method)

Функция get_polygons_for_date — ядро фильтрации:

JSON3.read парсит строку с GeoJSON.
Цикл for проходит по всем объектам в массиве features.
Для каждого объекта проверяется наличие свойства "date".
Дата преобразуется в объект Date и сравнивается с target_date.
Если даты совпали, извлекается геометрия. В зависимости от ее типа (Polygon или MultiPolygon) координаты добавляются в результирующий массив. Важно, что MultiPolygon разбивается, и его составляющие добавляются по отдельности.

function get_polygons_for_date(geojson_data::String, target_date::Date)
    data = JSON3.read(geojson_data)
    features = get(data, "features", [])
    result_polygons = []
    
    for feature in features
        properties = get(feature, "properties", nothing)
        if properties !== nothing && haskey(properties, "date")
            feature_date = Date(properties["date"])
            if feature_date == target_date
                geometry = get(feature, "geometry", nothing)
                if geometry !== nothing
                    geom_type = get(geometry, "type", "")
                    coordinates = get(geometry, "coordinates", [])
                    if geom_type == "Polygon"
                        push!(result_polygons, coordinates)
                    elseif geom_type == "MultiPolygon"
                        append!(result_polygons, coordinates)
                    end
                end
            end
        end
    end
    return result_polygons
end

get_polygons_for_date (generic function with 1 method)

Функция calculate_polygon_area реализует алгоритм «шнуровки» Гаусса (или формулу площадей барицентрических координат).

Формула проходит по всем вершинам полигона и суммирует векторные произведения координат. Результат дает удвоенную площадь полигона.
Адаптация для сферы: Это плоскостной алгоритм. Чтобы учесть сферичность Земли, координаты долготы (ring[i][1]) и широты (ring[i][2]) умножаются на приблизительные длины градусов в километрах (85.0 и 111.0 соответственно). Это грубое, но работающее приближение.

function calculate_polygon_area(coordinates)
    total_area = 0.0
    for polygon in coordinates
        for ring in polygon
            if length(ring) >= 3
                area = 0.0
                n = length(ring)
                for i in 1:n
                    j = (i % n) + 1
                    area += (ring[i][1] * 85.0) * (ring[j][2] * 111.0) - 
                           (ring[j][1] * 85.0) * (ring[i][2] * 111.0)
                end
                total_area += abs(area) / 2.0
            end
        end
    end
    return total_area
end

calculate_polygon_area (generic function with 1 method)

Функция draw_ukraine_map используется для отрисовки геометрии.

seriestype=:shape — ключевой параметр, указывающий, что мы рисуем замкнутую фигуру (полигон) по заданным координатам.
Использование plot! и scatter! с восклицательным знаком означает, что мы добавляем новые элементы к существующему графику.
annotate! добавляет текстовые подписи в заданные координаты.
Параметр aspect_ratio=:equal гарантирует, что масштаб по осям X и Y одинаков, и карта не будет искажена.

function draw_ukraine_map(start_polygons, current_polygons)
    plt = plot()
    for polygon in start_polygons
        for ring in polygon
            x_coords = [point[1] for point in ring]
            y_coords = [point[2] for point in ring]
            plot!(x_coords, y_coords, 
                  seriestype=:shape, 
                  color=:lightblue, 
                  alpha=0.4, 
                  linecolor=:blue,
                  linewidth=1,
                  label="")
        end
    end 
    for polygon in current_polygons
        for ring in polygon
            x_coords = [point[1] for point in ring]
            y_coords = [point[2] for point in ring]
            plot!(x_coords, y_coords, 
                  seriestype=:shape, 
                  color=:red, 
                  alpha=0.4, 
                  linecolor=:darkred,
                  linewidth=1,
                  label="")
        end
    end
    
    for (name, area, center) in REGIONS_DATA
        scatter!([center[1]], [center[2]], 
                 markersize=4, 
                 markercolor=:black, 
                 markerstrokewidth=0.5,
                 label="")
        
        annotate!(center[1], center[2], text(name, 6, :left))
    end
    
    plot!(plt,
        legend=:bottomright,
        aspect_ratio=:equal, 
        size=(1000, 800),
        ylims=(44, 52),
        title="Территориальные изменения с $(Dates.format(DATE_START, "dd.mm.yyyy")) по $(Dates.format(DATE_CURRENT, "dd.mm.yyyy"))",
        xlabel="Долгота",
        ylabel="Широта",
        grid=true
    )
    return plt
end

draw_ukraine_map (generic function with 1 method)

Теперь выполним запуск всего алгоритма целиком, в котором данные загружаются, фильтруются, анализируются, и результаты выводятся в виде текста и изображения.

geo_data = get_gzip_data(DSM_URL)
start_polygons = get_polygons_for_date(geo_data, DATE_START)
current_polygons = get_polygons_for_date(geo_data, DATE_CURRENT)
start_polygon_area = calculate_polygon_area(start_polygons)
current_polygon_area = calculate_polygon_area(current_polygons)

println("Период: с $(Dates.format(DATE_START, "dd.mm.yyyy")) по $(Dates.format(DATE_CURRENT, "dd.mm.yyyy"))")
println()
println("  Общая площадь Украины: $UKRAINE_TOTAL_AREA км²")
println("  Начальная площадь: $(round(start_polygon_area, digits=0)) км² ($(round((start_polygon_area/UKRAINE_TOTAL_AREA)*100, digits=2))%)")
println("  Текущая площадь:   $(round(current_polygon_area, digits=0)) км² ($(round((current_polygon_area/UKRAINE_TOTAL_AREA)*100, digits=2))%)")

draw_ukraine_map(start_polygons, current_polygons)

Период: с 24.02.2022 по 23.09.2025

  Общая площадь Украины: 603628 км²
  Начальная площадь: 0.0 км² (0.0%)
  Текущая площадь:   128781.0 км² (21.33%)

Вывод

Представленный код является отличным примером того, как на современном языке высокого уровня можно быстро создать эффективный инструмент для анализа сложных геопространственных данных.

Этот код служит прочной основой для реализации этих и многих других функций, демонстрируя гибкость и производительность Engee в области Data Science.