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

# Downloading and uploading country border data (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()).

# Download and unpack the map
# 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:

# Creating a map
plot([-180,180], [-90,-90], fillrange=[90,90], color=:lightblue)

# Let's draw the outline of the coast
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.:

# Adding points (example: London, New York, Moscow)
cities_lon = [-0.1276, -74.0060, 37.6176]
cities_lat = [51.5072, 40.7128, 55.7558]
city_names = ["London", "new york", "Moscow"]

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

# Adding a line (example: route)
route_lon = [-0.1276, 37.6176]  # London -> Moscow
route_lat = [51.5072, 55.7558]
plot!(route_lon, route_lat, color=:blue, linewidth=2, label="Route")

# We sign the names of cities
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

# Defining projections
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 (spherical)", "+proj=ortho +lat_0=30 +lon_0=0"),
    ("Eckert IV (pseudo-cylindrical)", "+proj=eck4 +lon_0=0"),
    ("Winkel Tripel (compromise)", "+proj=wintri +lon_0=0"),
    ("Sinusoidal (equilateral)", "+proj=sinu +lon_0=0"),
    ("Azimuthal Equidistant (azimuthal)", "+proj=aeqd +lat_0=0 +lon_0=0")
]

# Creating an array of subgraphs
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)
    
    # Specifying the current subgraph
    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
    
    # Adding the equator
    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")

    # Adding markers
    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.