# Combine Makie with GMT.jl

**URL:** <https://discourse.julialang.org/t/combine-makie-with-gmt-jl/74347>\
**Category:** Visualization\
**Tags:** makie, graphs, gmt\
**Created:** [January 10, 2022, 5:55pm UTC](https://discourse.julialang.org/t/combine-makie-with-gmt-jl/74347 "2022-01-10T17:55:01Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![filchristou](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/filchristou/32/26760_2.png) [@filchristou](https://discourse.julialang.org/u/filchristou)\
**Post date:** [January 10, 2022, 5:55pm UTC](https://discourse.julialang.org/t/combine-makie-with-gmt-jl/74347/1 "2022-01-10T17:55:01Z")

</div>

**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](http://sndlib.zib.de/home.action) (longtitude, latitude):

```julia
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:

```julia
GraphMakie.graphplot(gr; layout=mycoordlayout)

```

 ![image](https://global.discourse-cdn.com/julialang/original/3X/7/9/79a8306f3be74cee32bf54846e69c118fb9a460c.png)

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

_GMT.jl_  
Inspired from the code here [Map projections · GMT](https://docs.juliahub.com/GMT/EoU0j/0.30.1/proj_examples/#Lambert-conic-conformal), I do:

```julia
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)

```

 ![GMTjl_tmp](https://global.discourse-cdn.com/julialang/original/3X/2/d/2df5fdb8fb12cdd9f7bfc008652d804cd5a64df1.png)

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:

```julia
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 ?

---

<div class="post-metadata">

**Author:** ![joa-quim](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joa-quim/32/227_2.png) [@joa-quim](https://discourse.julialang.org/u/joa-quim)\
**Post date:** [January 10, 2022, 11:22pm UTC](https://discourse.julialang.org/t/combine-makie-with-gmt-jl/74347/2 "2022-01-10T23:22:18Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![filchristou](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/filchristou/32/26760_2.png) [@filchristou](https://discourse.julialang.org/u/filchristou)\
**Post date:** [January 11, 2022, 6:42am UTC](https://discourse.julialang.org/t/combine-makie-with-gmt-jl/74347/3 "2022-01-11T06:42:07Z")

</div>

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 ?

---

<div class="post-metadata">

**Author:** ![joa-quim](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joa-quim/32/227_2.png) [@joa-quim](https://discourse.julialang.org/u/joa-quim)\
**Post date:** [January 11, 2022, 1:55pm UTC](https://discourse.julialang.org/t/combine-makie-with-gmt-jl/74347/4 "2022-01-11T13:55:19Z")

</div>

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](https://github.com/GenericMappingTools/gmt/blob/master/test/psxy/geosegmentize.sh). That test is a bash shell but it can be useful to help translate into the GMT.jl syntax

 ![geosegmentize](https://global.discourse-cdn.com/julialang/original/3X/4/6/468ac7a42f27bdd12311bc92e985151370138784.png)

---

<div class="post-metadata">

**Author:** ![sdanisch](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sdanisch/32/1406_2.png) [@sdanisch](https://discourse.julialang.org/u/sdanisch)\
**Post date:** [January 11, 2022, 3:00pm UTC](https://discourse.julialang.org/t/combine-makie-with-gmt-jl/74347/5 "2022-01-11T15:00:13Z")

</div>

You could also try replicate the plot with GeoMakie:  
[https://github.com/JuliaPlots/GeoMakie.jl](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/](https://lazarusa.github.io/BeautifulMakie/GeoPlots/)

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

---

<div class="post-metadata">

**Author:** ![lazarusA](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lazarusa/32/6571_2.png) [@lazarusA](https://discourse.julialang.org/u/lazarusA)\
**Post date:** [January 11, 2022, 7:27pm UTC](https://discourse.julialang.org/t/combine-makie-with-gmt-jl/74347/6 "2022-01-11T19:27:12Z")

</div>

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.)

 ![graphMap](https://global.discourse-cdn.com/julialang/original/3X/2/1/21188180d1d70643d67ea479a31646d378c30a0e.png)

```julia
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

```

```julia
(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

```

---

<div class="post-metadata">

**Author:** ![Jeff\_Emanuel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jeff_emanuel/32/15440_2.png) [@Jeff\_Emanuel](https://discourse.julialang.org/u/Jeff_Emanuel)\
**Post date:** [January 11, 2022, 8:14pm UTC](https://discourse.julialang.org/t/combine-makie-with-gmt-jl/74347/7 "2022-01-11T20:14:07Z")

</div>

> [@lazarusA](#):
>
> `(-105.0, 40.75) => "DNVRng",`

FYI, your Denver is 1 degree too far north.

---

<div class="post-metadata">

**Author:** ![lazarusA](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lazarusa/32/6571_2.png) [@lazarusA](https://discourse.julialang.org/u/lazarusA)\
**Post date:** [January 11, 2022, 8:28pm UTC](https://discourse.julialang.org/t/combine-makie-with-gmt-jl/74347/8 "2022-01-11T20:28:02Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![filchristou](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/filchristou/32/26760_2.png) [@filchristou](https://discourse.julialang.org/u/filchristou)\
**Post date:** [January 18, 2022, 11:43am UTC](https://discourse.julialang.org/t/combine-makie-with-gmt-jl/74347/9 "2022-01-18T11:43:01Z")

</div>

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:

```julia
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:

 ![image](https://global.discourse-cdn.com/julialang/original/3X/1/b/1b935095bd991d3d02bd83abacea0e323a531f81.jpeg)

```julia
  [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:

 ![image](https://global.discourse-cdn.com/julialang/original/3X/8/4/84e995cfef684a031a9036ad84683659b57fef24.png)

```julia
 [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.

---

<div class="post-metadata">

**Author:** ![lazarusA](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lazarusa/32/6571_2.png) [@lazarusA](https://discourse.julialang.org/u/lazarusA)\
**Post date:** [January 18, 2022, 6:45pm UTC](https://discourse.julialang.org/t/combine-makie-with-gmt-jl/74347/10 "2022-01-18T18:45:22Z")

</div>

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:

```plaintext
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

```

 ![lcc_proj](https://global.discourse-cdn.com/julialang/original/3X/d/6/d60001a9e2a799683a02adf271cdde1ae25f0080.png)

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