Vector world map

In this example, we will show how to work with a vector map of the world: draw a line of coasts and plot several points in the geodetic coordinate system.

Loading and outputting vector maps

One of the most notable industry standards for storing vector geodetic information is the Shapefile standard. SHP files containing country outlines, roads and rivers, airport and port markings are easy to find on the Internet and download, including using a script. In this example, we will use a map obtained from the [Natural Earth] website (https://www.naturalearthdata.com/downloads/110m-physical-vectors /).

We will need packages Shapefile and Downloads (already available in the Engee supply):

Pkg.add("Shapefile")

In this example, the file will be available. If you need to automate the download of any other file, you can use the following set of commands.

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)

This code places the unpacked archive with maps in the current project folder (but instead of the macro @__DIR__ it is also convenient to use a temporary folder, for exampletemp_path = mktempdir()).

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

Since we don't always need to download the map again, we'll write a code for rendering a pre-downloaded map in SHP format.

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

From the entire file Shapefile we will be interested in a set of fields geometry. The command line or command helps to explore the contents of complex structures. dump(объект).

We uploaded a physical map of the world, and only the coast, in low resolution. Let's draw it:

# Создаем карту
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)

And add a few points and a route.:

# Добавляем точки (пример: Лондон, Нью-Йорк, Москва)
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!()

Map output in other projections

In order to display vector information in another projection, it is enough to create an object that will transform coordinates and apply it to each pair of coordinates. We will perform this operation using the library. Proj and we will build a map in several projections, also drawing the equator line and displaying the location of the points of interest on the map.

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)

Conclusion

In many projects (for example, those dedicated to navigation), information output performed in such a simple way may be useful. And having mastered the language of the Proj library, you can flexibly transform the coordinates of points and output them using the usual library tools. Plots.