几何多边形分析以评估领土变化

本文讨论了用于分析几何多边形的时间变化的算法的实现。该代码加载最新的地理数据，计算不同日期的领土面积，并可视化结果，展示了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在数据科学领域的灵活性和性能。