Векторная карта мира

В этом примере мы покажем, как работать с векторной картой мира: отрисуем линию побережий и нанесем несколько точек в геодезической системе координат.

Загрузка и вывод векторных карт

Одним из наиболее заметных промышленных стандартов для хранения векторной геодезической информации является стандарт Shapefile. Файлы SHP, содержащие контуры стран, дороги и реки, отметки аэропортов и портов, легко найти в интернете и скачать, в том числе при помощи скрипта. В этом примере мы будем пользоваться картой, полученной с сайта Natural Earth.

Нам потребуются пакеты Shapefile и Downloads (уже есть в поставке Engee):

Pkg.add("Shapefile")

В этом примере файл будет доступен. Если вам нужно автоматизировать скачивание какого-либо другого файла, можно воспользоваться следующим рядом команд.

using Shapefile
using Downloads

# Скачиваем и загружаем данные границ стран (Natural Earth)
function download_world_shapefile()
    shapefile_url = "https://naciscdn.org/naturalearth/110m/physical/ne_110m_coastline.zip"
    zip_path = joinpath(@__DIR__, "countries.zip")
    Downloads.download(shapefile_url, zip_path)
    run(`unzip -o $zip_path -d $(@__DIR__)/countries`)
    shp_path = joinpath(@__DIR__, "countries", "ne_110m_coastline.shp")
    return shp_path
end

download_world_shapefile (generic function with 1 method)

Этот код помещает распакованный архив с картами в текущую папку проекта (но вместо макроса @__DIR__ также удобно воспользоваться временной папкой, напримерtemp_path = mktempdir()).

# Скачиваем и распаковываем карту
#shp_path = download_world_shapefile()

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

shp_path = "$(@__DIR__)/countries/ne_110m_coastline.shp"
table = Shapefile.Table(shp_path)

Shapefile.Table{Union{Missing, Shapefile.Polyline}} with 134 rows and the following 4 columns:
	
geometry, scalerank, featurecla, min_zoom

Из всего файла Shapefile нас будет интересовать набор полей geometry. Исследовать содержимое сложных структур помогает командная строка или команда dump(объект).

Мы загрузили физическую карту мира, причем только побережья, в низком разрешении. Отрисуем ее:

# Создаем карту
plot([-180,180], [-90,-90], fillrange=[90,90], color=:lightblue)

# Нарисуем контур побережья
for shape in table
    x = [point.x for point in shape.geometry.points]
    y = [point.y for point in shape.geometry.points]
    plot!(x, y, color=:white, linewidth=0.5)
end

plot!(legend=false, grid=false, axis=false, framestyle=:none, 
        size=(1200, 600), xlim=(-180, 180), ylim=(-90, 90), aspect_ratio=:equal)

И добавим несколько точек и маршрут:

# Добавляем точки (пример: Лондон, Нью-Йорк, Москва)
cities_lon = [-0.1276, -74.0060, 37.6176]
cities_lat = [51.5072, 40.7128, 55.7558]
city_names = ["Лондон", "Нью Йорк", "Москва"]

scatter!( cities_lon, cities_lat, color=:red, markersize=5, label="Cities")

# Добавляем линию (пример: маршрут)
route_lon = [-0.1276, 37.6176]  # Лондон -> Москва
route_lat = [51.5072, 55.7558]
plot!(route_lon, route_lat, color=:blue, linewidth=2, label="Route")

# Подписываем названия городов
for c in zip(cities_lon, cities_lat, city_names)
    annotate!(c[1], c[2], text(c[3], 8, :bottom))
end

plot!()

Вывод карты в других проекциях

Для того, чтобы вывести векторную информацию в другой проекции, достаточно создать объект, который будет осуществлять трансформацию координат, и применить его к каждой паре координат. Мы выполним эту операцию при помощи библиотеки Proj и построим карту в нескольких проекциях, нарисовав также линию экватора и выведя на карту местонахождение интересующих нас точек.

Pkg.add("Proj")

   Resolving package versions...
  No Changes to `~/.project/Project.toml`
  No Changes to `~/.project/Manifest.toml`

using Proj

# Определяем проекции
projections = [
    ("Mercator", "+proj=merc +lon_0=0 +lat_ts=0 +x_0=0 +y_0=0"),
    ("Robinson", "+proj=robin +lon_0=0 +x_0=0 +y_0=0"),
    ("Mollweide", "+proj=moll +lon_0=0 +x_0=0 +y_0=0"),
    ("Orthographic (сферическая)", "+proj=ortho +lat_0=30 +lon_0=0"),
    ("Eckert IV (псевдоцилиндрическая)", "+proj=eck4 +lon_0=0"),
    ("Winkel Tripel (компромиссная)", "+proj=wintri +lon_0=0"),
    ("Sinusoidal (равновеликая)", "+proj=sinu +lon_0=0"),
    ("Azimuthal Equidistant (азимутальная)", "+proj=aeqd +lat_0=0 +lon_0=0")
]

# Создаем массив подграфиков
plt = plot(layout=(ceil(Int32, length(projections) / 2),2), size=(800, 1200))

for (i, (name, proj)) in enumerate(projections)
    trans = Proj.Transformation("+proj=longlat", proj)
    
    # Указываем текущий подграфик
    subplot = plt[i]
    
    for shape in table
        x = [point.x for point in shape.geometry.points]
        y = [point.y for point in shape.geometry.points]
        xy_transf = trans.(x, y)
        xx,yy = getindex.(xy_transf, 1), getindex.(xy_transf, 2)
        plot!(subplot, xx, yy, c=1, lw=0.5, fill=:lightblue, title=name, legend=false, grid=false, framestyle=:none)
    end
    
    # Добавляем экватор
    lons = range(-180, 180, length=100)
    lats = zeros(100)
    line_transf = trans.(lons, lats)
    x_line, y_line = getindex.(line_transf, 1), getindex.(line_transf, 2)
    plot!(subplot, x_line, y_line, c=:red, lw=2, label="Equator")

    # Добавляем отметки
    cities_transf = trans.(cities_lon, cities_lat)
    cities_x = getindex.(cities_transf, 1)
    cities_y = getindex.(cities_transf, 2)
    scatter!(subplot, cities_x, cities_y, color=:red, markersize=5, label="Cities")
    for c in zip( cities_x, cities_y, city_names)
        annotate!(subplot, c[1], c[2], text(c[3], 8, :bottom))
    end
    
end

display(plt)

Заключение

Во многих проектах (например, посвященных навигации) может пригодиться вывод информации, выполненный таким несложным способом. А освоив язык библиотеки Proj можно вполне гибко преобразовывать координаты точек и выводить их привычными средствами библиотеки Plots.