GMT.jl : Coloring lakes, keeping oceans bathymetry

Specific Domains Geo
1

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:

stations_NA
stations_NA1865×1023 539 KB

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)

stations_NA2
stations_NA21865×1023 558 KB

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

Thanks!

2

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)

GMTjl_j
GMTjl_j1865×1023 465 KB

4 Likes
3

Thanks! I works very well:

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

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

stations_NA3
stations_NA31865×1023 555 KB

4

Just two quick notes.

One do not normally need to set show=false because that is the default for almost all functions (so that we can overlay like in current example). The exceptions are the viz function that is meant to be a quick visualizer ad hence has the show=true by default. Other exceptions are the demo functions like seismicity and alike.

Adding shade=true to the grdimage command makes an even nicer image (IMO).

1 Like
5

Thanks, I will try!

6

I agree that it is much better with the shade=true option!

image
image1863×1018 750 KB

7

Nice. And actually I over complicated this issue. There is an option in coast exactly for this (river_fill) but its long names parsing is broken, so one must use condensed syntax. The two coast! calls above can be replaced by a simpler:

coast!(shore=true, area=200, river_fill="skyblue+l")

and when a fix is committed, that can be expressed as:

coast!(shore=true, area=200, river_fill=(lake=true, fill=:skyblue))
1 Like