Slow grdcontour for seafloor age with NetCDF + .cpt

I’m probably doing something wrong, but I’m having some trouble with grdcontour:

  • it’s very slow, compared to gmt(grdcontour ...)
  • I don’t know how to reproduce the default behaviour of gmt(grdcontour ...)

First some version info:

julia> versioninfo()
Julia Version 1.6.2
Commit 1b93d53fc4* (2021-07-14 15:36 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: Intel(R) Core(TM) i7-8705G CPU @ 3.10GHz
  LIBM: libopenlibm
  LLVM: libLLVM-12.0.1 (ORCJIT, skylake)

(@v1.6) pkg> st GMT
      Status `~/.julia/environments/v1.6/Project.toml`
  [5752ebe1] GMT v0.36.1 ``

julia> GMT.GMTbyConda

shell> gmt --version

The files ( grid and 1x.cpt color map) were too large (is there a way to compress NetCDF files better than 1% from zip?) but there are similar ones: grid file, colormap. I used grdcut to cut the grid file to match the region.

The code:

using GMT

function t1()
    gmt("grdcontour -JM14c -R110/210/-40/40 -Nseafloor_age.cpt")

function t2()
    seafloor_age = gmtread("")
    cpt = gmtread("seafloor_age.cpt")
        region = (110, 210, -40, 40), 
        proj = :Mercator, 
        fill = true,
        contours = cpt,  # MethodError when trying to pass the filename here directly.
        annot = 5,  # Also tried without this, no change to output.

Using t2 is very slow and only a couple of labels are drawn. I don’t think annot is correct, I also tried without that argument. The first issue is that I’m unclear how to reproduce the plain gmt -N defaults. The other difference visually is that the map produced by t2 looks more “smoothed” in color gradients, maybe there is some change in defaults for the smoothing option?

I also tried without the modern mode gmt{begin,end} but the issues seem to persist.

Output from t2:

OK, the t2 function is slower because you read the entire grid first and only later apply the region, whilst t1 reads a sub-region directly.

I’m not sure why one cannot give a C=ctpt_name.cpt. @pwessel, is it because the KEY says CC( and thus always expects a cpt type?

A grid at 1 min (the link you provided is broken). I have one at 2 min but seems to have issues with the scale/offset values so I cannot try it. However, with an older 6 m one it works fine

grdcontour("C:/j/cd_grids/", region=(110, 210, -40, 40), proj=:merc, 
contours=gmtread("age_2020.cpt"), fill=true, annot=5, show=true)


  1. grdcontour fill=true is an heavy op. It’s much faster if you use a grdimage+grdcontour sequence.
  2. Not sure why your annot case did not work well and I’ll need to check it with that big grid, but I’m short of time right now.
  3. It would be good if you could update to GMT6.2. There were issues with grdcontour -N -C before.

My result

grscontour -C command line can deal with different types of input (CPT, list of contours, etc) but not sure if you pass a memory address - I think then it must be a CPT structure, no?

Probably. I think I’ll have to read the cpt file under the hood to workaround this situation.

Thanks for the fast response as always, sorry about the link (should be fixed now).

the t2 function is slower because you read the entire grid first

Ah, is there a way to only read regions from the grid file?

It’s much faster if you use grdimage+grdcontour

Thanks, I’ll read more about grdimage in that case.

I’ll see what happens when I fiddle with the annotations some more. I’ll look into compiling GMT 6.2 as well (my Linux distribution is still on 6.1).

A quick follow up with a simpler comparison:

julia> using BenchmarkTools

julia> function seafloor_contours_monolithic()
           gmt("begin seafloor_contours_monolithic pdf")
           gmt("grdcontour -JM14c -Baf -BWSen -R110/210/-40/40")

julia> function seafloor_contours_modern()
           grdcontour(gmtread(""), proj = :Mercator, region = (110, 210, -40, 40))

julia> function test_seafloor_contours()
           @btime seafloor_contours_modern()
           @btime seafloor_contours_monolithic()

julia> test_seafloor_contours()
  110.530 s (1399874805 allocations: 21.73 GiB)
  48.320 s (66 allocations: 3.69 KiB)


As mentioned in the issue the modern mode doesn’t plot axes by default (I forgot to add them manually…). So that’s a different story, but I’m more concerned about the performance. It makes sense that the second case uses way less Julia memory (because the data is not passed through Julia, right?) and I guess the difference in time comes from the issue with selecting the region as you said before…

Yes, use grdcut

Regarding the performance. The best is to always to use the file names directly when we are dealing with data on disk. When we pass in memory data the array has normally to be transposed, which implies a copy of that array. The transposition is due to the fact that C and Julia use different memory layouts (C is row-major). I try to avoid that in some known cases by keeping the row-major order when receiving the arrays from C, but this also creates some risky situations. And the pad is another copy-forcing-factor. So, in seafloor_contours_modern() just use grdcontour("", ...)

1 Like

Thanks, I thought that might be the case, it makes sense based on our discussions. I will use grdcut and then plot from the smaller grid file directly. I have been able to improve the speed a lot by using grdimage as well, as you suggested (which also works well with the “.cpt” file). There is one last thing: when using the file path directly, the region is not inherited from previous plotting commands:

function doesntwork()
    coast(region = (110, 210, -40, 40), proj = :Mercator, frame = :auto, shore = :thinner)
    grdimage!("", cmap = "seafloor_age.cpt")
    grdcontour!("", annot = 5)

function works()
    coast(region = (110, 210, -40, 40), proj = :Mercator, frame = :auto, shore = :thinner)
    grdimage!("", region = (110, 210, -40, 40), cmap = "seafloor_age.cpt")
    grdcontour!("", region = (110, 210, -40, 40), annot = 5)

It’s not a big deal: I can save the region to a variable and not have to type it twice. But, I’m just letting you know in case it’s a bug or could be fixed. It seems like the projection is inherited, for example.

(Also if you have time: any hints on how to stop the contour labels from overlapping?)

In modern mode one do not use the ! form, so use grdimage instead of grdimage!. See if makes the difference (it should inherit the region).

See the option G or label in grdcontour manual (and check the plain GMT version when details are missing in the jl version), but wont promise it will solve this sometimes difficult problem.

Unfortunately it’s not working:

julia> function seafloor_age_fullmodern()
           coast(region = (110, 210, -40, 40), proj = :Mercator, frame = :auto, land = :dimgray, shore = :thinner)
           grdimage("", cmap = "seafloor_age.cpt")
           grdcontour("", annot = 5)
seafloor_age_fullmodern (generic function with 1 method)

julia> seafloor_age_fullmodern()
grdimage [ERROR]: Option -R:  Cannot include south/north poles with Mercator projection!
grdimage [ERROR]: Internal Failure = error
grdimage (GMT_grdimage): (null)
ERROR: Something went wrong when calling the module. GMT error number = 73
 [1] gmt(::String, ::Nothing, ::Vararg{Nothing, N} where N)
   @ GMT ~/.julia/packages/GMT/lcweZ/src/gmt_main.jl:283
 [2] finish_PS_module(::Dict{Symbol, Any}, ::Vector{String}, ::String, ::Bool, ::Bool, ::Bool, ::Nothing, ::Vararg{Nothing, N} where N)
   @ GMT ~/.julia/packages/GMT/lcweZ/src/common_options.jl:3399
 [3] grdimage(cmd0::String, arg1::Nothing, arg2::Nothing, arg3::Nothing; first::Bool, kwargs::Base.Iterators.Pairs{Symbol, String, Tuple{Symbol}, NamedTuple{(:cmap,), Tuple{String}}})
   @ GMT ~/.julia/packages/GMT/lcweZ/src/grdimage.jl:120
 [4] seafloor_age_fullmodern()
   @ Main ./REPL[7]:4
 [5] top-level scope
   @ REPL[8]:1

I tried it first but I thought maybe I need the ! form. Anyway, I’ll use explicit regions for now, maybe I’ll get some time to look into this as well and open an issue.

Hmm, bringing @htyeim here. You are a modern mode partizan and this should certainly be a basic breaking thing if it happens all the times. Did you never get this error?

The other thing that I’ve now noticed is that the coastlines are missing. They appear only if the call to coast is made after grdimage.

The region error also happens for grdcontour, but not for e.g. coast:

julia> gmtbegin()

julia> basemap(region = (110, 210, -40, 40), proj = :Mercator, frame = :auto)

julia> grdcontour("", annot = 5)
grdcontour [ERROR]: Option -R:  Cannot include south/north poles with Mercator projection!
grdcontour [ERROR]: Internal Failure = error
ERROR: Something went wrong when calling the module. GMT error number = 73
 [1] gmt(::String, ::Nothing, ::Vararg{Nothing, N} where N)
   @ GMT ~/.julia/packages/GMT/lcweZ/src/gmt_main.jl:283
 [2] finish_PS_module(::Dict{Symbol, Any}, ::Vector{String}, ::String, ::Bool, ::Bool, ::Bool, ::Nothing, ::Vararg{Nothing, N} where N)
   @ GMT ~/.julia/packages/GMT/lcweZ/src/common_options.jl:3399
 [3] grdcontour(cmd0::String, arg1::Nothing; first::Bool, kwargs::Base.Iterators.Pairs{Symbol, Int64, Tuple{Symbol}, NamedTuple{(:annot,), Tuple{Int64}}})
   @ GMT ~/.julia/packages/GMT/lcweZ/src/grdcontour.jl:147
 [4] top-level scope
   @ REPL[43]:1

julia> gmtend()

julia> gmtbegin()

julia> basemap(region = (110, 210, -40, 40), proj = :Mercator, frame = :auto)

julia> coast(land = :dimgray, shore = :thinner)

julia> gmtend(show = true)

Things are working fine in classic mode.

Absolutely expected. coasts must be plotted on top of images, otherwise they’ll go under de blanket.

OK for completeness, I’m posting a version that works using the current GMT.jl master, and also features:

  • modern mode
  • specifying the region only once
  • sane label spacing
  • fast (using grdimage)

The files are a seafloor age grid (NetCDF) and a color map (CPT), see original question for download links. The code:

function seafloor_age()
        cmap = "seafloor_age.cpt",
        region = (110, 210, -40, 40),
        proj = :Mercator,
        frame = :auto,
    # This also works, but there were too many labels for my taste:
    # grdcontour("", contours = "seafloor_age.cpt")
    grdcontour("", contours = 5, annot = 25)
    coast(land = :dimgray, shore = :thinner)

The output:

1 Like