GMT.jl : Coloring lakes, keeping oceans bathymetry

[Balinus]

Hello! I’m trying to plot some points over a map. I would like to keep the default topography/bathymetry visuals, but the lakes are “green”, being over 0m at sea level (my guess). Here’s my current code. I tried to use the coast function, but it overrides the ocean visuals.

using GMT

grdimage("@earth_relief_01m", proj=(name=:lambertConic, center=[-97.5 48], parallels=[40 50]), region=(-140,-55,32,64), coast=true, show=false)
GMT.coast!(water="skyblue", shore=true)

GMT.plot!(station, marker=:circ, ms=0.23, mc=:white, markeredgecolor=:red, fmt=:png, figname=joinpath(repfig, "stations_NA.png"))

This gives the following figure:

While, if I omit the coast function, it gives a nice ocean visuals:

grdimage("@earth_relief_01m", proj=(name=:lambertConic, center=[-97.5 48], parallels=[40 50]), region=(-140,-55,32,64), coast=true, show=false)

GMT.plot!(station, marker=:circ, ms=0.23, mc=:white, markeredgecolor=:red, fmt=:png, figname=joinpath(repfig, "stations_NA.png"))

Which gives the following (lakes are green)

I tried inverting the call of coast and grdimage without success.

Thanks!

[joa-quim]

The key here is to understand that coastlines in GMT are a set of polygons with different hierarchical levels (see brief here) and option area. Together, we can split it in to calls. The first plots only the coastlines and the second plots only the lakes and hence we paint/pen style them differently. See

grdimage("@earth_relief_01m", proj=(name=:lambertConic, center=[-97.5 48], parallels=[40 50]), region=(-140,-55,32,64))

# Plot only the ocean-land interface (polygon level 1)
coast!(shore=true, area=(0,1,1))

# Plot only the coast-lakes interface (lakes are polygons level 2)
# The the 200 means that polygons with area < 200 km2 are nor drawn) 
coast!(shore=true, water=:red, area=(200,2,2), show=true)