How to use "grdclip" or other functions to clip my contour based on the build in bathymetric data?

Hi All,

Sorry for the many questions. I get excited about GMT and really want to plot something I’ve been dreaming about.

My final goal is a subsurface filled contour of an oceanographic variable, along with the land as well as the continental shelf/slope between the land and my variable contour. Below are what I’m trying to do specifically:
(a) Plot a map showing the land, as well as the continental shelf/slop, e.g., at 500m water depth. I’m able to do that using grdclip, thanks to help from Joaquim.
(b) What I’m trying to figure out next (in this thread) is how to clip off part of my variable contour along the 500m-isobath line?

My variable is on a 1440 x 720 longitude and latitude grid. Will it be possible to clip off the part of the variable contour that falls into the region landward of the 500m-isobath.

Below is my current code to clip off the land (leaving water) along the coastline. So basically, instead of clipping off my variable contour along the coastline, I want the clipping to be done along the 500m-isobath line.

# 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"), show=true )
    coast!(clip=:end)

G2 is my gridded variable created using mat2grid.

Any advice is greatly appreciated!

Use grdclip to make a mask from the bathymetry grid that has z < -500 = NaN and z >= -500 = 1. Then multiply this mask with your G2 grid. You may have to use grdsample to make the mask grid compatible with G2 since you cannot multiply two grids if they are not exactly compatible in terms of origin, increments a grid registration.

1 Like

Thank you, Joaquim.

Below is what I did:

# put the Earth relief data onto a new grid defined by my variable grid, G2:
G_ocean = grdsample("@earth_relief_10m", G2);

# Create a mask based on your advice:
mask1 = grdclip(G_ocean, low=[-500 NaN], high=[-500 1]);

# Multiply the mask with my variable grid to derive the new grid:
G2_clip = G2 * mask1;

There is an error about the last line of code:
Grids have different sizes, so they cannot be multiplied.

It seems that I wrote my grdsample code wrong? Unfortunately, I’m unable to find an example code for grdsample to verify.

Please let’s not start with this again. From grdsample — GMT 6.6.0 documentation


which clearly tells that one has to use the limits (R or region) and the grid increment (I or ``increment). What it doesn’t say in this case is if we have to also change the grid registration (I mentioned in a Julia post that understanding the different registrations is fundamental

So your if your G2 and the @earth_relief_10m do not have the same origin, increment and registration (which the probably do) the grdsample should resample the mask to have exactly the same parameters of your G2 grid. That is the region the increment and, eventually, the registration (registration)

1 Like

My apologies for not knowing that I can use the same command in Julia as if I’m using GMT in a terminal? I thought Julia would need something like this instead?
G_ocean = grdsample("@earth_relief_10m", region=(180,320,10,80), increment=0.25, G2);

My G2 is on a 0.25 decimal degrees grid: longitude: 20:0.25:379.75, and latitude:-90:0.25:89.75. How do I specify the increment in this case then?

That is explained in the GMT.jl documentation. See


and here
grdsample("@earth_relief_10m". ...)

No. You want to resample the mask to be consistent with your G2 grid or vice-versa so that you can multiply the two.

Thanks, Joaquim! I agree with you in the proposed steps. Unfortunately, what I’m struggling with is trying to figure out how to write a proper grdsample syntax. I have decided to give up on this route and get around that issue by using my own elevation grid data E. Now I have to move on to the next issue :grinning:

# A grid data of the global ocean:
E_grid = mat2grid(Array(transpose(E)), x=lon, y=lat);

# Create a mask:
mask1 = grdclip(E_grid, low=[-500 NaN], high=[-500 1]);

G2_clip = G2 * mask1;