Support irregular grid

[htyeim]

Irregular grid is commonly used in many simulation programs.
I think it would be necessary to support plotting the results of these simulations.
Below is an example output of an ionosphere simulation program (SAMI2). Code and data are attached in the end.

The grid can be plotted using plot3, but I haven't found a way to plot data using GMT.

Matplotlib provides pcolor/pcolormesh which can plot the data in irregular grid easily.

I think it just add some polygons to the figure ref.

I found a post might be related to this. I think there is no need to modify grdimage.
A new function named pcolor/pcolormesh would be great.

data used:
alts.txt
glats.txt
glons.txt
ne.txt

code used above:

using Libdl
"/usr/local/Cellar/gmt/6.0.0_5/lib" in Libdl.DL_LOAD_PATH ||
push!(Libdl.DL_LOAD_PATH, "/usr/local/Cellar/gmt/6.0.0_5/lib")
using GMT
using DelimitedFiles

function build_AGMTgrid(alts, glats, glons)
    grd1 = Array{GMT.GMTdataset,1}()
    grd2 = Array{GMT.GMTdataset,1}()
    n1, n2 = size(glons)
    for i in 1:n1
        if i % 3 != 1 continue end
        push!(grd1, GMT.GMTdataset(hcat(glons[i,:], glats[i,:], alts[i,:])))
    end
    for i in 1:n2
        if i % 2 != 1 continue end
        push!(grd2, GMT.GMTdataset(hcat(glons[:,i], glats[:,i], alts[:,i])))
    end
    grd1, grd2
end



tmp_array = Array{Array{Float64,2},1}()
for fn in ["alts","glats","glons","ne", ]
    push!(tmp_array,
        DelimitedFiles.readdlm("$fn.txt", '\t', Float64, '\n'))
end

alts, glats, glons, ne = tmp_array

grd1, grd2 = build_AGMTgrid(alts, glats, glons)

R =  (65, 160, -80, 80, 0, 4000)
view = (150, 45)

proj = (name = :PlateCarree, center = (0, 20),)

gmtbegin()
gmtfig("test1",fmt = "png,jpg")
basemap(region = R, proj = proj,
            frame = (axes = "WSE"),
            xaxis = (annot = :auto, label = :lon),
            yaxis = (annot = :auto, label = :lat),
            view = view, 
            figsize = 20, zsize = 20,)

topo = makecpt(color = :globe, range = (-8800, 8800,), continuous = true)

grdimage("@earth_relief_06m.grd", region = R, view = view, 
        shade = (azim = view[1], norm = "e0.8"),color = topo)

plot3(grd1,region = R, view = view,
    zsize = 20,
    pen = (0.5, :red, :solid)
)
plot3(grd2,region = R, view = view,
    frame = (axes = "Z1234"),
    zsize = 20,
    zaxis = (annot = 1000, label = :altitude),
    pen = (0.3, :cyan, :solid),
)
gmtend()

using PyPlot

PyPlot.pcolormesh(glats,alts,ne,cmap = "jet")
PyPlot.savefig("test2.png", )

Are you willing to help implement and maintain this feature? Yes

[welcome]

👋 Thanks for opening your first issue here! Please make sure you filled out the template with as much detail as possible. We appreciate that you took the time to contribute!

Please make sure you read our Contributing Guide and abide by our Code of Conduct.

[joa-quim]

Hi @htyeim Thanks for this example and question raised. I'm in the middle of other things for today, but I think the GMT module contour (or pscontour) should be able to do that. I'm not certain how pcolormesh works but suspect it's not very different. contour will do a triangulation and lets us give a color palette to fill the contours.

[htyeim]

I tried to use contour and got one figure:

julia> GMT.contour(ds, region = R[1:4], show = true, fmt = :png,
               fill = true, 
               # cont = 1e6,# color = cb_jet, # cont = 1,
               nan_alpha = true, figsize = (12, 12),
               Vd = 1
               )
        pscontour  -R-27.0/44.0/84.0/4000.0 -JX12/12 -Baf -BWSen -I -C -P -K > /var/folders/ss/65r_tg917hg47rkrfjxk9ftc0000gn/T/GMTjl_tmp.ps

It might be better to plot the value at each grid point (along and perpendicular to the geomagnetic field lines). The SAMI2 can capture disturbances with very small scale, the contour will eliminate these useful information. But I agree with that a pcolormesh function should be very similar to contour.
In Matplotlib, the key sentences of pcolormesh might be:

       C = C.ravel()
        coords = np.column_stack((X, Y)).astype(float, copy=False)
        collection = mcoll.QuadMesh(Nx - 1, Ny - 1, coords,
                                    antialiased=antialiased, shading=shading,
                                    **kwargs)

I feel that the data structures in GMT are very different from those of Matplotlib so I don't know if it is hard to implement the function (or modify from contour).

[joa-quim]

Nice plot, and looking at it there are two things that come up to me

1. The GMT plot is way nicer than the Matplotlib one 😄
2. Your second matplotlib fig (the one zoomed) shows lozenges, which makes think that some gridding process took place under the hood. The GMT fig on the other hand has sharp clean edges.

Supporting irregular grids in plain GMT is an old request but hard to implement. Workarounds via the Julia wrapper seem more practical. To make it plot at grid points one would have to, under the hood, interpolate into a regular grid and plot it. That implies a time consuming step and a considerable amount of guessing. Isn't a a wrap-around contour a better option?

[htyeim]

.....So GMT already support irregular grid? contour with dense interval:

rg = (0, 2.1e6, 1e4)

cb_jet = makecpt(color = :jet, range = rg, 
                    background = (:blue, :red),
                    # continuous = false,
                    );
R = (-27.0, 44.0, 84.0, 4000.0)
GMT.contour(ds, region = R, show = true, fmt = :png,
    fill = true, 
    cont = cb_jet,
    range = rg,
    nan_alpha = true, figsize = (12, 12),
    Vd = 1
)

With grid:

Zoom:

I think the only thing needed is to update document.. When I want to plot some figures similar to above and if I check document, I will only try to use grdimage not contour. Thank you for this great project!

[htyeim]

And I have another question. As GMT can process this kind of data and plot it, GMT acts like a data processing program. Is there a way to get the data after triangulation? I think it is very useful.

[joa-quim]

.....So GMT already support irregular grid?

No, but a recent adding to GMT master was to allow grd2xyz to output the correct x,y,z coordinates of irregular grids so we can pipe them to contour and make it work. It's a bit of a wasteful way because it will have to compute a triangulation when simple rectangles would do but ok.

Is there a way to get the data after triangulation?

Yes, see the triangulate module.

think the only thing needed is to update document..

Yes, documentation is needed for lots of modules and so far the best docs are those of GMT itself.

[joa-quim]

Just a final comment. It's a pity that this is not in the Forum as it shows very nice usage examples and here, as an issue, they will be hidden soon when we close this.

[htyeim]

All data and code here can be used to build an example in the document if you think it is useful~ I am still learning how to use GMT and GMT.jl. The more I learned the more I found I like it.

[joa-quim]

Thanks. It's the time and slot in priorities that worry me. But I'll do it.