How to plot continental shelf/slope floor next to the land

Q&A GMT.jl Q&A
1

As we all know, due to the existence of continental shelves and continental slopes, the ocean areas get smaller when it gets deeper. Below is my program to plot some contours in the ocean. It is mainly for the surface. In the surface, I use clip to draw the boundary of the land and ocean, but for the subsurface plots, the boundary of the land and ocean will move seaward correspondingly.

In my example, how do I allow the program to do the below?
(1) to continue to plot the land the same way as is,
(2) to plot the continental shelf/slope using a predefined color (For example, if I’m plotting a contour of my variable at 500m, this should be the region between the depth of 0 and -500m)
(3) to clip the contour using the new land/sea boundary at the deeper depth.

Below is my current code:

# Plot the land map, clipping water leaving land:
    coast(region=(180,320,10,80), clip=:land)
    grdimage!("@earth_relief_10m", shade=true)
    coast!(clip=:end)

# Plot the contour of my variable, clipping land and leaving water
    coast!(clip=:water)
    cpt = makecpt(cmap=:jet, range=(0,5,0.02) );
    grdimage!(G2, coast=true, color=cpt, frame=(annot=45,grid=45,title="test") )
    coast!(clip=:end, show=1)

G2 above is the gridded data of my variable using mat2grid.

For testing purposes, my G2 can be generated through:
lon = 20:0.25:379.75;
lat = -90:0.25:89.75;
Z = rand(0:5, 720, 1440);
G2 = mat2grid(Z, x=lon, y=lat);

Many thanks!

2

What can I say? You are getting what you are asking. We have no idea what is the G2 contents. Have you plotted it? Are you sure it’s correct?

3

Thanks for the reply, Joaquim! I’m able to plot the map as well as the contour. It is exactly what I want for the surface. Where I need help right now is for the subsurface plots, so that I could draw the continental shelf/slop area in a predefined color as well. That way, I do not want my contour to extend to the now land area (would be considered as part of the ocean in the surface plot) at a certain depth level.

You said, I’m already getting what I’m asking? Were you referring to the continental shelf/slope area? But I did not even find a place to specify the depth anywhere, am I right?

For testing purposes, my G2 can be generated through:
lon = 20:0.25:379.75;
lat = -90:0.25:89.75;
Z = rand(0:5, 720, 1440);
G2 = mat2grid(Z, x=lon, y=lat);

4

Don’t think grids with random numbers are useful in the context of plotting the continental shelf, slope and raise.
You have two options to do that

  • 1 Create a suitable color map (a cpt) with colors that show well those feature. Use grdimage to do the plot.
  • 2 Create (or get them) polygons that delimit those provinces and plot them colorfully using plot.

There are no clipping operations involved in this task.

1 Like
5

Thank you! Sounds like I would be able to extract what I need out of the same Earth_relief_10m data, and do something like the below?
grdimage!("@earth_relief_10m", ...)

  1. How could I define the polygon based on the depth range of 0m and -500m?
  2. You mentioned no “clipping” is needed. This is interesting, I remember you mentioned previously that unless we poke holes in these images before overlaying them, they will block each other.
  3. So I would need to use the new polygon to block part of the contour that is extended to the continental shelf/slope (now part of the land)?
6

Maybe you should read this
https://www.sciencedirect.com/science/article/abs/pii/S0025322714000310?via%3Dihub

1 Like
7

Is there a way to address this issue through the “range” as below, i.e., defining the cpt based on the depth range of 0m to -500m?

cpt = makecpt(range=(-500,0));
grdimage!("@earth_relief_10m", cmap=cpt)

Anyone could show me how to fit this into my code so that (a) the rest of the contour beyond 0m and -500m depth will be clipped, and (b) this new layer will cover part of my variable contour?

Also, is it possible to set the color to one solid color, like a shade of gray? Thanks!

8

If you don’t want anything below -500 m use grdclip … or plain Julia.

You said you’re going to read the docs :slight_smile:

9

Many thanks!

Sorry but reading the docs will not happen overnight. It will take me some time, unfortunately. Not to mention that I find the GMT.jl documentation really really hard to follow. In comparison, the Matlab documentation seems to be 100 times easier.

I also find it interesting that what I’m trying to achieve, i.e., making a contours in the subsurface, the need to draw the land between the coastline and the other side of the continental shelf/slope, is not something people have already done so many times previously.

10

When I say read the docs, that includes the GMT (not GMT.jl) docs too. I know GMT.jl docs are incomplete but that’s when one have to resort to the help of the other ones.

And haven’t I showed already how to plot contours on the ocean side only?

11

Hi @leon6j, if the Matlab documentation was on the same level as GMT documentation, perhaps you would complain that the $1000+ subscription for Matlab was a bit high compared to the 0$ model we operate under.
Anyway, I think your 100x is a bit exaggerated but I agree it is much better - amazing what a commercial company can do given a few hundred people working full-time on something!

12

He’s (rightfully aside the price) complaining on GMT.jl docs.

13

My apologies if my comments on the documentation were not considered constructive. To make that up, here are some real recommendations for you to consider:

  1. Instead of having a function like the below, could we have something more real like G2=grdclip(G, above=[n1 n2], low=[n3 n4], between="n5/n6/n7"), where G is the input matrix, G2 is the output grid, all values above n1 will be replaced by n2, etc.

grdclip(cmd0::String="", arg1=nothing, kwargs...). I have no idea what these things are.

  1. Can we make sure that for each function, instead of just showing an example like the below without much explanation, you would provide a fully independently workable example that people can run on their computer and actually see the map? Could we add documentation explaining what each of the argument does, e.g., “above=[5 6]”, replace=ing all values below 5 with 6, or vice versa.
    G=gmt(“grdmath”, “-R0/10/0/10 -I1 X”);
    G2=grdclip(G, above=[5 6], low=[2 2], between=“3/4/3.5”)

  2. BTW, I’m one of the few who much prefer to pay you thousands of dollars in exchange for dedicated assistance.

Above all, I really admire this wonderful tool you have built for the community, and the fact that you are offering it to all of us for free.

15

Hi Joaquim,

Thank you for the recommendation of using “grdclip”. Unfortunately, the colormap of my new map has changed from the normal “geo” to some shade of blue.

# Here was my original code:
    coast(region=(180,320,10,80), clip=:land)
    grdimage!("@earth_relief_10m", shade=true, show=true)
    coast!(clip=:end)

# Here is my new code, after using "grdclip" to cut off the data with a water depth of -500m or deeper:
    GG = grdclip("@earth_relief_10m", low=[-500 NaN])
    coast(region=(180,320,10,80))
    grdimage!(GG, shade=true, coast=true, show=true)

Do you know what I did wrong?

16

@leon6j I know what you mean and I should at least try to replace those generic module (cmd0::String="", arg1=nothing, kwargs...) by something more similar to what you are suggesting (though I need to least *all options in the synopsis). But the point I tried to explain earlier is that many modules (grdclip is one of them) do not yet have a Julia syntax manual page. When they do, they look like this, which has also many examples


but when they don’t one have to consult the GMT doc man page and adapt for the Julia syntax as briefly explained in the module doc strings. Or not adapt at all and use the GMT terse syntax given in text strings.
1 Like
17

Now you have two different figures (you call show twice) but since they have the same default name (you didn’t save them with different names by using the savefig option) you are only seeing the second.
@earth_relief grids know they are topo-bathymetric grids and automatically use the topo color map. Others do not know that and use GMT’s default color palette, which is turbo.
Answer: use cmap=:topo in your seconds grdimage command.

Note also that a grid at 10 arc minutes resolution is too coarse to show the continental shelf in most parts of the world.

1 Like
18

Many thanks, Joaquim!

That was exactly the problem. I’ll try to figure our how to mark this issue resolved. As to the resolution issue, I plan to use 1 arc minutes resolution for my final product.

19

Users are not supposed to have this as their first encounter with the docs. What is supposed is that the better documented modules teach users to recognize that, otherwise cryptic, instructions.

But I think we have made our points though I’ll add still a typical extra answer that I hate to receive. If you make those suggestions in GMT.jl issues, even if only for a more comprehensible synopsis, the documentation effort will advance faster.

1 Like
20

I have moved my recommendations about the GMT.jl documentation to the GMT.jl page. Sorry for posting them in the wrong place, Joaquim.