Slow grdcontour for seafloor age with NetCDF + .cpt

adigitoleo · September 21, 2021, 3:00pm

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
  WORD_SIZE: 64
  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 `https://github.com/GenericMappingTools/GMT.jl.git#master`

julia> GMT.GMTbyConda
false

shell> gmt --version
6.1.1

The files (1x.nc 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()
    gmtbegin("seafloor_age.pdf")
    gmt("grdcontour -JM14c -R110/210/-40/40 seafloor_age_sample.nc -Nseafloor_age.cpt")
    gmtend()
end

function t2()
    seafloor_age = gmtread("seafloor_age_sample.nc")
    cpt = gmtread("seafloor_age.cpt")
    gmtbegin("seafloor_age.pdf")
    grdcontour(
        seafloor_age, 
        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.
    )
    gmtend()
end

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:

Joaquim · September 21, 2021, 4:44pm

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/age.2.6.nc", region=(110, 210, -40, 40), proj=:merc, 
contours=gmtread("age_2020.cpt"), fill=true, annot=5, show=true)

Notes:

grdcontour fill=true is an heavy op. It’s much faster if you use a grdimage+grdcontour sequence.
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.
It would be good if you could update to GMT6.2. There were issues with grdcontour -N -C before.

My result

pwessel · September 21, 2021, 6:53pm

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?

Joaquim · September 21, 2021, 7:41pm

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

adigitoleo · September 22, 2021, 1:56am

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).

adigitoleo · September 22, 2021, 9:05am

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 seafloor_age_sample.nc")
           gmt("end")
       end

julia> function seafloor_contours_modern()
           gmtbegin("seafloor_contours.pdf")
           grdcontour(gmtread("seafloor_age_sample.nc"), proj = :Mercator, region = (110, 210, -40, 40))
           gmtend()
       end

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

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

Outputs:

seafloor_contours.pdf (105.5 KB)
seafloor_contours_monolithic.pdf (205.3 KB)

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…

Joaquim · September 22, 2021, 9:55am

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("seafloor_age_sample.nc", ...)

adigitoleo · September 22, 2021, 11:29am

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()
    gmtbegin("seafloor_age.pdf")
    coast(region = (110, 210, -40, 40), proj = :Mercator, frame = :auto, shore = :thinner)
    grdimage!("seafloor_age_sample.nc", cmap = "seafloor_age.cpt")
    grdcontour!("seafloor_age_sample.nc", annot = 5)
    gmtend()
end

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

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?)

Joaquim · September 22, 2021, 11:53am

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.

adigitoleo · September 22, 2021, 12:34pm

Unfortunately it’s not working:

julia> function seafloor_age_fullmodern()
           gmtbegin("seafloor_age.pdf")
           coast(region = (110, 210, -40, 40), proj = :Mercator, frame = :auto, land = :dimgray, shore = :thinner)
           grdimage("seafloor_age_sample.nc", cmap = "seafloor_age.cpt")
           grdcontour("seafloor_age_sample.nc", annot = 5)
           gmtend()
       end
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
Stacktrace:
 [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.

Joaquim · September 22, 2021, 1:23pm

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?

adigitoleo · September 22, 2021, 1:35pm

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.

adigitoleo · September 22, 2021, 2:08pm

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("seafloor_age_sample.nc", 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
Stacktrace:
 [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.

Joaquim · September 22, 2021, 2:25pm

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

adigitoleo · September 22, 2021, 3:27pm

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()
    gmtbegin("seafloor_age.pdf")
    grdimage(
        "seafloor_age_sample.nc",
        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("seafloor_age_sample.nc", contours = "seafloor_age.cpt")
    grdcontour("seafloor_age_sample.nc", contours = 5, annot = 25)
    coast(land = :dimgray, shore = :thinner)
    gmtend()
end

The output: