Combine Makie with GMT.jl

filchristou · January 10, 2022, 5:55pm

Short story:
I want to plot graphs on top of a map. Can I use GMT.jl to draw the region of interest and plot over it using GraphMakie.jl ?

Long story:
GraphMakie.jl
Let’s say I want to plot the Abilene network.
I have my coordinates data taken from the SNDlib (longtitude, latitude):

Dict{Tuple{Float64, Float64}, String} with 12 entries:
  (-96.5967, 38.9617) => "KSCYng"
  (-122.3, 47.6)      => "STTLng"
  (-105.0, 40.75)     => "DNVRng"
  (-122.026, 37.3858) => "SNVAng"
  (-87.6167, 41.8333) => "CHINng"
  (-85.5, 34.5)       => "ATLAng"
  (-77.0268, 38.8973) => "WASHng"
  (-73.9667, 40.7833) => "NYCMng"
  (-86.1595, 39.7806) => "IPLSng"
  (-95.5174, 29.77)   => "HSTNng"
  (-118.25, 34.05)    => "LOSAng"
  (-84.3833, 33.75)   => "ATLAM5"

Using a simple layout function I can get a realistic graph visualization:

GraphMakie.graphplot(gr; layout=mycoordlayout)

But I would like to have as a background the USA map.

GMT.jl
Inspired from the code here Map projections · GMT, I do:

using GMT
coast(region=[-130 -70 24 52], proj=(name=:lambertConic, center=[-100 35], parallels=[33 45]),
                    frame=:ag, borders=((type=1, pen=("thin","red")), (type=2, pen=("thinner",))),
                    land=:tan, water=:darkblue)
scatter!(coordslong, coordslat, markersize=0.2, mc=:darkgreen, show=true)

So the scatter does overwrite the image and I get my nodes at the right coordinations.
I guess I can continue going and create the edges between the nodes as lines, but I am hoping for a better solution.

Desirable
I would like to be able to do something like:

GMT.coast(region=[-130 -70 24 52], proj=(name=:lambertConic, center=[-100 35], parallels=[33 45]),
                    frame=:ag, borders=((type=1, pen=("thin","red")), (type=2, pen=("thinner",))),
                    land=:tan, water=:darkblue)
GraphMakie.graphplot!(gr; layout=mycoordlayout)

and get my graph on top of the map.

I gues that means (correct me if I am wrong) combing 2 different plotting backends on the same plot; would that be possible ?

joa-quim · January 10, 2022, 11:22pm

Hi @filchristou I’m afraid the short answer is no. GMT.jl creates PS fig files that are later converted to raster formats. I have no idea what Makie intermediate steps do but they certainly cannot append to the PS file generated by GMT.jl.

However, CairoMakie uses Cairo that can produce PostScript so it’s not inconceivable that CairoMakie could one day append to an existing PostScript file, but … not today.

filchristou · January 11, 2022, 6:42am

Hi! Not surprised. I was kind of hoping for the Julia out-of-the-box inter-package magic.
So I guess the most straightforward way to plot a graph on top of a GMT.jl map, is to indeed indeed use scatter() and arrows().
Wouldn’t it be cool to multiple dispatch graphplot() in order to do it automatically and add it to the toolset ?

joa-quim · January 11, 2022, 1:55pm

Maybe you want to explore the connection option (-F in GMT syntax) of the plot module. It’s not an easy option to describe but you can do things like the attached image. This image comes from one of the GMT tests. That test is a bash shell but it can be useful to help translate into the GMT.jl syntax

sdanisch · January 11, 2022, 3:00pm

You could also try replicate the plot with GeoMakie:
https://github.com/JuliaPlots/GeoMakie.jl

GeoMakie isn’t in the best shape, but @lazarusA has done some amazing plots with it already:
https://lazarusa.github.io/BeautifulMakie/GeoPlots/

There is also a release pending, which should make GeoMakie much more user friendly

lazarusA · January 11, 2022, 7:27pm

Ok. It’s not that hard to do with GeoMakie.jl (well I did it during dinner ), it will be better soon, after the pending release is out. (xticks, yticks could be also added manually.)

using GeoMakie, GLMakie
using GraphMakie, Graphs
using Proj4
using Downloads
using GeoJSON, GeoInterface

states = Downloads.download("https://raw.githubusercontent.com/PublicaMundi/MappingAPI/master/data/geojson/us-states.json")
states_geo = GeoJSON.read(read(states, String))
n = length(GeoInterface.features(states_geo))
#g = wheel_graph(10)
gpos = Dict(
    (-96.5967, 38.9617) => "KSCYng",
    (-122.3, 47.6) => "STTLng",
    (-105.0, 40.75) => "DNVRng",
    (-122.026, 37.3858) => "SNVAng",
    (-87.6167, 41.8333) => "CHINng",
    (-85.5, 34.5) => "ATLAng",
    (-77.0268, 38.8973) => "WASHng",
    (-73.9667, 40.7833) => "NYCMng",
    (-86.1595, 39.7806) => "IPLSng",
    (-95.5174, 29.77) => "HSTNng",
    (-118.25, 34.05) => "LOSAng",
    (-84.3833, 33.75) => "ATLAM5")

g = complete_graph(length(keys(gpos)))
positions = Point2f.(collect(keys(gpos)))


fig = Figure(resolution = (1200, 800))
ax = Axis(fig[1, 1], aspect = DataAspect())
trans = Proj4.Transformation("+proj=longlat +datum=WGS84",
    "+proj=lcc +lon_0=-100 +lat_1=33 +lat_2=45", always_xy = true)
ptrans = Makie.PointTrans{2}(trans)
ax.scene.transformation.transform_func[] = ptrans
lons = -130:-70
lats = 24:52
points = [Point2f0(lon, lat) for lon in lons, lat in lats]
rectLimits = FRect2D(Makie.apply_transform(ptrans, points))
limits!(ax, rectLimits)
lines!(ax, GeoMakie.coastlines(), color = :black)
poly!(ax, states_geo, color = 1:n, colormap = (:viridis, 0.25), strokecolor = :black,
    strokewidth = 1)
graphplot!(ax, g; layout = _ -> positions, node_size = 20,
    edge_color = to_colormap(:plasma, 66),
    node_color = to_colormap(:plasma, length(keys(gpos))))
# drawing the correct lines and ticks is still a little bit cumbersome and manual
# you could ignore the rest and just hide everything
lonrange = range(-130, -70, length = 5)
latrange = range(24, 52, length = 5)

lonlines = [Point2f0(j, i) for i in lats, j in lonrange]
latlines = [Point2f0(j, i) for j in lons, i in latrange]

[lines!(ax, lonlines[:, i], color = (:black, 0.25),
    linestyle = :dash, overdraw = true) for i in 1:size(lonlines)[2]]
[lines!(ax, latlines[:, i], color = (:black, 0.25), linestyle = :dash,
    overdraw = true) for i in 1:size(latlines)[2]]

hidedecorations!(ax)
hidespines!(ax)
save("graphMap.png", fig)
fig

(proj4Geo) pkg> st
      Status `~/Desktop/proj4Geo/Project.toml`
  [13f3f980] CairoMakie v0.6.6
  [e9467ef8] GLMakie v0.4.7
  [ddc7317b] GeoDatasets v0.1.6
  [cf35fbd7] GeoInterface v0.5.6
  [61d90e0f] GeoJSON v0.5.1
  [db073c08] GeoMakie v0.2.2
  [5c1252a2] GeometryBasics v0.4.1
  [1ecd5474] GraphMakie v0.3.0
  [86223c79] Graphs v1.5.1
  [ee78f7c6] Makie v0.15.3
  [85f8d34a] NCDatasets v0.11.7
  [9a7e659c] Proj4 v0.7.6

Jeff_Emanuel · January 11, 2022, 8:14pm

FYI, your Denver is 1 degree too far north.

lazarusA · January 11, 2022, 8:28pm

it’s not my Denver, I just copy/paste the info from above , but good to know.

filchristou · January 18, 2022, 11:43am

While trying to simplify this script for GeoMakie v0.3.0 not only I couln’t get it working but also I realized that the same script you provided deesn’t work for GeoMakie v0.3.0.

For example the following extraction from your script:

let
	states = Downloads.download("https://raw.githubusercontent.com/PublicaMundi/MappingAPI/master/data/geojson/us-states.json")
	states_geo = GeoJSON.read(read(states, String))
	n = length(GeoInterface.features(states_geo))
	
	fig = Figure(resolution = (1200, 800))
	ax = Axis(fig[1, 1], aspect = DataAspect())
	trans = Proj4.Transformation("+proj=longlat +datum=WGS84",
	    "+proj=lcc +lon_0=-100 +lat_1=33 +lat_2=45", always_xy = true)
	ptrans = Makie.PointTrans{2}(trans)
	ax.scene.transformation.transform_func[] = ptrans
	lons = -130:-70
	lats = 24:52
	points = [Point2f0(lon, lat) for lon in lons, lat in lats]
	rectLimits = FRect2D(Makie.apply_transform(ptrans, points))
	limits!(ax, rectLimits)
	lines!(ax, GeoMakie.coastlines(), color = :black)
	poly!(ax, states_geo, color = 1:n, colormap = (:viridis, 0.25), strokecolor = :black, strokewidth = 1)
	fig
end

works for GeoMakie v0.2.2:

  [13f3f980] CairoMakie v0.6.6
  [e9467ef8] GLMakie v0.4.7
  [cf35fbd7] GeoInterface v0.5.7
  [61d90e0f] GeoJSON v0.5.1
  [db073c08] GeoMakie v0.2.2
  [1ecd5474] GraphMakie v0.3.0
  [86223c79] Graphs v1.5.1
  [682c06a0] JSON v0.21.2
  [ee78f7c6] Makie v0.15.3
  [9a7e659c] Proj4 v0.7.6

but not for v0.3.0:

 [13f3f980] CairoMakie v0.7.0
  [e9467ef8] GLMakie v0.5.0
  [cf35fbd7] GeoInterface v0.5.7
  [61d90e0f] GeoJSON v0.5.1
  [db073c08] GeoMakie v0.3.0
  [1ecd5474] GraphMakie v0.3.0 `https://github.com/JuliaPlots/GraphMakie.jl.git#compathelper/new_version/2022-01-08-00-39-37-422-00161005115`
  [86223c79] Graphs v1.5.1
  [682c06a0] JSON v0.21.2
  [ee78f7c6] Makie v0.16.1
  [9a7e659c] Proj4 v0.7.6

Do you mind posting the v0.3.0 equivalent ?
Don’t mind the graphplot.

lazarusA · January 18, 2022, 6:45pm

Hi,
unfortunately, the new projection method in GeoMakie is not working as intended. Hence, the output doesn’t always look good. Hence, I personally still use the old version.

The idea was to have a simpler interface, as the following:

using GeoMakie, Proj4, GLMakie
using Downloads
using GeoJSON, GeoInterface

states = Downloads.download("https://raw.githubusercontent.com/PublicaMundi/MappingAPI/master/data/geojson/us-states.json")
states_geo = GeoJSON.read(read(states, String))
n = length(GeoInterface.features(states_geo))

fig = Figure(resolution = (1200, 800), fontsize = 22)
ax = GeoAxis(fig[1, 1], source = "+proj=longlat +datum=WGS84",
    dest = "+proj=lcc +lon_0=-100 +lat_1=33 +lat_2=45",
    title = "Projection: lcc +lon_0=-100 +lat_1=33 +lat_2=45", coastlines = true,
    lonticks = -130:5:-70, latticks = 20:5:50)
poly!(ax, states_geo, color = 1:n, colormap = (:viridis, 0.25),
    strokecolor = :black, strokewidth = 1)
fig

Note that the ticks don’t work well, plus the right limit is not correct (compare to the previous answer).