How to plot a grid with "filled contours"

Added by Juan Jose 7 months ago

Hi, I'd like to achieve a result that is rather common in most plotting programs, but I could not figure out how to achieve with GMT.

I have a netCDF file with some field. I can plot it with grdimage. But I get the raw matrix, which is pixeled and without any sort of "smoothing" to produce more "meaning" regions:

What I'd like to get is some sort of "countoured" colours. Something similar to what I'd get with grdcontour, but filling the region within contours with the same colour. Something similar to example 12, achieved with pscontour, but operating with netCDF fies.

I'm sure there must be a way to do that. Any idea?


Replies (18)

RE: How to plot a grid with "filled contours" - Added by Peter 7 months ago

Hi Juan

I agree the resulting image does not look "nice", however it is a conservative way of displaying data, after all you do not really know the data value except in those points where it has been measured (or calculated or ...). Any type of contouring of discrete data introduces information that is simply not known.

That said, what you can do is to interpolate the data onto a finer grid, have a look at grdsample to operate directly on the netCDF file, or alternatively (for a tunable interpolation algorithm) dump the data to an xyz file and use surface

RE: How to plot a grid with "filled contours" - Added by Juan Jose 7 months ago

Thanks for the fast reply. The surface approach is what I normally do, but I though there could be a more straightforward way to achieve the same results.

I agree that the pixels are the most accurate and conservative way of presenting the data. But the readability of the figure also matters, I believe. It is very difficult to read the values displayed (what is the value of correlation above the Alps, for instance?). With fewer colours in a number of discrete contours, one can more easily read the numbers. But this is a matter of taste, so I'm sure we could always find arguments pro and against . In this case, I'm pretty sure the "smooth" approach would benefit the overall result.

RE: How to plot a grid with "filled contours" - Added by Peter 7 months ago

So grdsample should be your friend then. BTW, I put in a feature request for grdsample to also have access to the interpolation algorithm used in surface in the future (at least from the man page it does not seem available currently)

RE: How to plot a grid with "filled contours" - Added by Juan Jose 7 months ago

Well, if we ask for new features, I'd ask an option to allows grdimage, or perhaps grdcontour, to do what pscontour already does, but working with netCDF instead of ascii as input files.

RE: How to plot a grid with "filled contours" - Added by Peter 7 months ago

I guess you mean filling of the contours, no? You have my support here, but I leave it to you to formulate the feature request

RE: How to plot a grid with "filled contours" - Added by Paul 7 months ago

Such plots are usually done with grdview -Qs.

RE: How to plot a grid with "filled contours" - Added by Juan Jose 6 months ago

Thanks Paul, that is exactly the option/approach I was missing.

To provide a constructive critic, perhaps it's worth to note that this is a bit hidden, or at least I find it counterintuitive. To me, grdview is "the tool to produce eye-candy 3D figures". I would expect that grdimage, or perhaps grdcontour would allow achieving this behaviour, as it is what you get with pscontour. I find it logical to expect grdcontour to be the "input binary counterpart" of pscontour.

But perhaps I'm misunderstanding the underlying rationale of the expected usage of the different tools...

RE: How to plot a grid with "filled contours" - Added by Paul 6 months ago

No, it is clearly a feature that belongs to grdcontour, and I will see if I could isolate a call to grdview from grdcontour for something like this. At the very least I can add some comments to the doc page.

RE: How to plot a grid with "filled contours" - Added by Paul 6 months ago

Implemented via a new -N[<cpt>] option in trunk.

RE: How to plot a grid with "filled contours" - Added by Juan Jose 6 months ago

Hi, I have a new issue that is very related this topic.

I'm trying to plot a field using the "grdimage -Qs" approach. It seems to work, but for some reason it produces an empty layer without any actual plot, despite the fact that the output of the command produces some postscript code. grdcontour and grdimage produce both the expected output, so it seems that the netCDF file is not the problem. I'm becoming crazy because I can not see what's wrong with this, and I'm starting to think it could be a bug.

Could you please double check?

I'm using gmt 5.4.3. I attach the script and netCDF files I use.

plot_WINDSLP_12km_area3.sh Magnifier - script (852 Bytes)

MV10ER-000001.nc - input file (152 KB)

RE: How to plot a grid with "filled contours" - Added by Joaquim 6 months ago

It's grdview -Qs, not grdimage -Qs

RE: How to plot a grid with "filled contours" - Added by Juan Jose 6 months ago

Sorry, I made a mistake writing the post. But in the attached script it's right, I'm using grdview -Qs (I also tried -Qm) without success.

RE: How to plot a grid with "filled contours" - Added by Joaquim 6 months ago

OK, it took me a little while to figure it out because your example grid raised a couple of issues. The main issue is that the z_min/z_max are NaN.
Tried to convert it with grdconvert but is crashed because the comment line is too long, so I had to use something like this

grd2xyz MV10ER-000001.nc?U10 | xyz2grd -Glixo.grd -RMV10ER-000001.nc

and than the grdview -Qs worked.

RE: How to plot a grid with "filled contours" - Added by Juan Jose 6 months ago

Thanks Joaquim. But I don't fully understand the origin of the issue. Where are those values NaN? The netCDF file as all values properly defined, as far as I can check with my usual tools...

RE: How to plot a grid with "filled contours" - Added by Joaquim 6 months ago

The problem is the the cdo command did not do a complete job and grdview is not checking for this eventuality. See

grdinfo MV10ER-000001.nc
MV10ER-000001.nc: Title: Longitude
MV10ER-000001.nc: Command: Wed Mar 14 13:19:37 2018: cdo splitrec seldate,2011-11-06T14:00:00,2011-11-08T23:00:00 MV10ER.nc MV10ER
Wed Mar 14 13:19:37 2018: cdo -sellonlatbox,0,13,34,46 ../DATOS/LIN/STORM-14-NOSALT/MV10ER.nc MV10ER.nc
Wed Mar 14 12:27:11 2018: cdo setgrid,LIN/STORM-14-NOSALT/SLP.nc LIN/STORM-14-NOSALT/MV10ER.nc kk.nc
Wed Mar
MV10ER-000001.nc: Remark:
MV10ER-000001.nc: Gridline node registration used [Cartesian grid]
MV10ER-000001.nc: Grid file format: nf = GMT netCDF format (32-bit float), CF-1.7
MV10ER-000001.nc: x_min: -1476 x_max: -252 x_inc: 12 name: x coordinate of projection [km] n_columns: 103
MV10ER-000001.nc: y_min: -612 y_max: 816 y_inc: 12 name: y coordinate of projection [km] n_rows: 120
MV10ER-000001.nc: z_min: NaN z_max: NaN name: Longitude [degrees_east]
MV10ER-000001.nc: scale_factor: 1 add_offset: 0
MV10ER-000001.nc: format: classic

RE: How to plot a grid with "filled contours" - Added by Paul 6 months ago

Placing

rm -f gmt*

in a script you share with us is a bit risky since I just deleted all my gmt C programs. Good thing we have version control.

RE: How to plot a grid with "filled contours" - Added by Paul 6 months ago

Addressed in r19934. If the grid is not following COARDS so zmin/xmax == NaN I reset those to 0/0 to prevent damage to z-scalings before data is even read.

RE: How to plot a grid with "filled contours" - Added by Juan Jose 6 months ago

Paul wrote:

Placing
[...]
in a script you share with us is a bit risky since I just deleted all my gmt C programs. Good thing we have version control.

Ups, I'm very sorry. I never thought that could be an issue.

Thank you for having that issue into account for future versions.

(1-18/18)