Geometric polygon analysis to assess territorial changes
This article discusses the implementation of an algorithm for analyzing temporal changes in geometric polygons. The code loads up-to-date geodata, calculates the area of territories for different dates, and visualizes the results, demonstrating Engee's capabilities for geospatial information processing tasks.
The following packages are used in this example:
Downloads: For easy file upload by URL.JSON3: A high-performance package for parsing JSON data.CodecZlib: Provides work with compressed data, in particular, GZIP decompression.Dates: Working with dates for filtering and formatting.Plots: A powerful and flexible library for plotting and visualizing data, used for rendering maps.GeoJSON(although added, it is not directly used in the code): A package for working with geodata. In this script, parsing is performed viaJSON3, which is also a valid approach.
In [ ]:
Pkg.activate("MyGeoProject", shared=true)
Pkg.add(["JSON3", "CodecZlib", "GeoJSON"])
using Downloads, JSON3, CodecZlib, Dates, GeoJSON
The code starts with defining constants:
DATE_START— a fixed date for the start of the analysis.DATE_CURRENT— the current date calculated at the time of the script launch. This makes the analysis always relevant.DSM_URL— The URL indicating the data source.REGIONS_DATA- contains background information on the regions of Ukraine: name, area and coordinates for the caption on the map.UKRAINE_TOTAL_AREA— the total area of the country.
In [ ]:
DATE_START = Date("2022-02-24")
DATE_CURRENT = Date(Dates.now())
Out[0]:
In [ ]:
DSM_URL = "https://raw.githubusercontent.com/cyterat/deepstate-map-data/main/deepstate-map-data.geojson.gz"
Out[0]:
In [ ]:
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
Function get_gzip_data:
Downloads.download(url)uploads the file by URL to a temporary file and returns the path to it.open(temp_file) do fileopens this temporary file.GzipDecompressorStream(file)creates a decompression stream that "wraps" the file. When reading from this stream, the data is automatically decompressed.read(stream, String)reads the entire decompressed stream into a string, which is returned by the function.
It is a solution for working with compressed data directly in memory, without having to manage temporary files manually.
In [ ]:
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
Out[0]:
Function get_polygons_for_date — the filtration core:
JSON3.readparses a string with GeoJSON.- The cycle
forpasses through all objects in the arrayfeatures. - The presence of a property is checked for each object.
"date". - The date is converted to an object.
Dateand is compared withtarget_date. - If the dates match, the geometry is extracted. Depending on its type (
PolygonorMultiPolygon) coordinates are added to the resulting array. It is important thatMultiPolygonit is split, and its components are added separately.
In [ ]:
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
Out[0]:
Function calculate_polygon_area implements the Gauss lacing algorithm (or the formula for the areas of barycentric coordinates).
- The formula goes through all the vertices of the polygon and summarizes the vector products of the coordinates. The result is twice the area of the polygon.
- Adaptation for the sphere: This is a planar algorithm. To account for the sphericity of the Earth, the longitude coordinates (
ring[i][1]) and latitude (ring[i][2]) multiplied by the approximate length of degrees in kilometers (85.0and111.0respectively). This is a rough but working approximation.
In [ ]:
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
Out[0]:
Function draw_ukraine_map it is used for drawing geometry.
seriestype=:shape— a key parameter indicating that we are drawing a closed shape (polygon) at the specified coordinates.- Usage
plot!andscatter!with an exclamation mark, it means that we are adding new elements to the existing schedule. annotate!adds text captions to the specified coordinates.- Parameter
aspect_ratio=:equalensures that the scale on the X and Y axes is the same, and the map will not be distorted.
In [ ]:
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
Out[0]:
Now we will run the entire algorithm, in which the data is loaded, filtered, analyzed, and the results are displayed as text and images.
In [ ]:
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)
Out[0]:
Conclusion
The presented code is an excellent example of how an effective tool for analyzing complex geospatial data can be quickly created in a modern high-level language.
This code provides a solid foundation for implementing these and many other features, demonstrating Engee's flexibility and performance in the field of Data Science.