GMT5 has lost the ability to read netcdf-3 (older format) .grd files

Hi all,

I am having a problem with some netcdf .grd files that I created with MATLAB using Kelsey Jordahl’s grdwrite2 script. I’ve been using Kelsey’s MATLAB-GMT interface scripts for literally decades, and they’ve always worked fine. But now GMT5 won’t recognize them as legit netcdf files. Note the following output from gmt5 grdinfo:

MP-McGovern-14:LPSC_frantic mcgovern$ gmt5 grdinfo svv_out.grd

grdinfo (api_import_grid): NetCDF: Attempting netcdf-4 operation on netcdf-3 file [svv_out.grd]

[Session gmt5 (0)]: Error returned from GMT API: GMT_GRID_READ_ERROR (18)

[Session gmt5 (0)]: Error returned from GMT API: GMT_GRID_READ_ERROR (18)

[Session gmt5 (0)]: Error returned from GMT API: GMT_GRID_READ_ERROR (18)

This has never happened before, and only started happening after I implemented the joint install of GMT 5 and GMT6 from MacPorts. That changed the “prefix” for using GMT 5 from simply “gmt” to “gmt5”, but apparently some actual functionality changed as well.

I will note that GMT 6 can read the files:

MP-McGovern-14:LPSC_frantic mcgovern$ gmt6 grdinfo svv_out.grd

grdinfo [WARNING]: The y-coordinates and range attribute are in conflict; must rely on coordinates only

svv_out.grd: Title: svv_out.grd

svv_out.grd: Command: File written by MATLAB function grdwrite2.m

svv_out.grd: Remark: Created 29-Jul-2020 14:16:06

svv_out.grd: Gridline node registration used [Cartesian grid]

svv_out.grd: Grid file format: nd = GMT netCDF format (64-bit float), CF-1.7

svv_out.grd: x_min: 0 x_max: 1866 x_inc: 2 name: x n_columns: 934

svv_out.grd: y_min: -50 y_max: -0 y_inc: 1 name: y n_rows: 51

svv_out.grd: z_min: -30.5628037562 z_max: 0.00118786475074 name: z

svv_out.grd: scale_factor: 1 add_offset: 0

svv_out.grd: format: classic

But I’m not out of the woods yet, because using such files in programs like grdimage gives me problems I’ve never seen before:

grdimage [WARNING]: The y-coordinates and range attribute are in conflict; must rely on coordinates only

…and I end up not getting useful output from it; the .ps file won’t display on my Mac because Preview thinks it’s not a proper .ps file.

So, I’m wondering if any of you can give me some insight into what’s going on here. I will note that I’ve already seen responses of “Just use Python widget XLBLORB (or whatever) to convert to netcdf4”, and that’s a good last resort, but I’ve got literally hundreds of such files that I need to use, and laboriously finding and converting them all individually is something I am hoping (and reasonably, I think, expecting) to avoid. In other words, I’m hoping that there’s a solution that simply restores GMT 5’s previous ability to handle these files.

Thanks in advance for your help,
Pat

Hi Pat-

The netcdf messages were suddenly added by NetCDF developers earlier this year and threw us as well [in fact they also broke something else as well for us for a while but that got fixed]. They decided to act on a function return code for a call they had not done before and poof suddenly make it a fatal error. Very irritating, So we had to change the code - GMT 6 has this change. GMT 5 obviously is not self-changing so it will be caught unless you build with an older netcdf library. And we clearly have no capacity or interested in issues patches for old GMT versions.
The warnings you see are harmless - it just means the netcdf attributes weren’t done well or complete enough and we have to snoop around to sort things out.
The moral is that even if you wish to stay with an older (and you may believe, more stable, GMT version) the reality is that GMT depends on many other libraries, like netcdf, hdf, curl, etc., etc., and those guys are changing too, and when things break or get crazy it is better to be at GMT 6.x.y with the rest of us trying to sort this out than to be stuck at GMT 5.x.y and have far less help since none of the developers have run GMT 5 since we started on 6 many years ago! It is very difficult to stay fixed in time with all the packages that are needed for a particular version. Best, Paul

Hi Paul,

Thanks for your response. I agree with you about not expending effort on older versions. I’d be quite happy to take this as a signal to move on to GMT 6 (I only started using 5 a year or so ago!). But the problem is that GMT 6 isn’t properly handling the netcdf-3 files either. As I noted, there is the “The y-coordinates and range attribute are in conflict; must rely on coordinates only” which, yes, is only a warning, but unfortunately the final product is still corrupted (the file produced is not viewable). I will note that this line is interesting:

svv_out.grd: y_min: -50 y_max: -0 y_inc: 1 name: y n_rows: 51

There’s no reason for the y_max 0 value to have a minus sign. I think that this may be an indication of what the problem is. But if so it is probably a netcdf problem rather than a GMT one.

So at that point, I would ask the forum if they know of any successor to the “grdwrite2” MATLAB script or other modality to write netcdf-4 grds out of MATLAB.

Thanks,
Pat

Hi Pat-

  1. Send me one of the netcdf3 files so I can have a look.
  2. Well, there is the GMT/MEX toolbox for a full GMT/MATLAB interaction, but certainly more work than a short mex file.

Hi Paul,

Can I attach a file in this forum? Or should I use email?

98% of my graphical output is some combination of GMT and MATLAB, so I am MASSIVELY interested in that GMT/MEX toolbox. But I took a brief look at it quite a while ago, and I couldn’t figure out how to get started. And yes, I am looking at quicker solutions in the short term.

You can post here unless unreasonably large file. THe reason the GMT/MEX toolbox is hard for us to deliver and maintain is that Matlab ships with all the same libraries we use (netcdf, etc) but much older, leading to endless library conflicts… Oh, and they dont ship include files, so we cannot even build against their obsolete versions. It is a mess.

I have a 390 KB file. Do I just drag it into this window?

Re: the GMT/MEX toolbox, wow, that sounds like a lot of headaches. I may just have to have GMT making ALL of the grd files.

Yes, I think you just drop it in one of these messages.

OK, here it is. This one is more or less lithostatic stress in a lithosphere, so it’s monotonic increasing in magnitude (negative, because compressive) with depth.

Well, I tried to drag it, but it has a limited list of extensions (graphics types), and it won’t take a grd file. But it will take a zip file. I gzip-ed it then changed the extension to .zip. That sounds like something weird to do to an already binary file, but here goes…

svv_out.grd.zip (320.6 KB)

Nothing wrong with this netcdf grid. It is GMT 4 so 2-D. The ymax = -0 is in the netcdf file you created (so probably just a tad short of an exact 0; here is the header dump (ncdump -h svv_out.grd):

netcdf svv_out {
dimensions:
x = 934 ;
y = 51 ;
variables:
double x(x) ;
x:long_name = “x” ;
x:actual_range = 0., 1866. ;
double y(y) ;
y:long_name = “y” ;
y:actual_range = -50., -0. ;
double z(y, x) ;
z:long_name = “z” ;
z:_FillValue = NaN ;
z:actual_range = -30.5628037561865, 0.00118786475074117 ;

// global attributes:
:Conventions = “COARDS/CF-1.0” ;
:title = “svv_out.grd” ;
:history = “File written by MATLAB function grdwrite2.m” ;
:description = “Created 29-Jul-2020 14:16:06” ;
:GMT_version = “4.x” ;
}

The warning is harmless but looks like it should not happen in this case - so I will look at that. Plots fine for me.

OK, thanks for the readout! There are a couple of leads here for me to follow, things related to how I create the file in MATLAB, that I can check on later tonight. But first, a proposal is due so I have to do some writing (ugh).

OK, I corrected the “feature” that provided “-0” to the vector of y values in the grds. So the minus sign goes away, but I still get

grdinfo [WARNING]: The y-coordinates and range attribute are in conflict; must rely on coordinates only

at the top of the output when calling “gmt6 grdinfo” on the files. And they are still not converging to a plotable .ps file when I use grdimage.

The warning is annoying and will go away, but it does not cause any issues otherwise. I have just submitted a pull request [https://github.com/GenericMappingTools/gmt/pull/4251] to fix that part. The reason it happens is that the grid says y goes from -50 to 0 but the vector stored goes from 0 to -50. I can try to reproduce your non-plot with, what GMT version? 6.0? 6.1.1?

Thanks! I just this afternoon upgraded on MacPorts to their latest, 6.1.1_0 (it didn’t fix the problem, apparently)

Do you have a one-liner that should work but fails for you? For me, this plots right away:

gmt grdimage svv_out.grd -B -pdf map

Let me make one and I’ll get back to you with the result.

Ooo, plots up nice with that dirt-simple one-liner!

Part of my problem may be that I’m still using somewhat complex GMT 4-style bash scripts for all my GMT stuff, and there may be baggage from that that is causing the problem. For example, the plots I’m trying to make have 3 sub-panels with lots of hairy -Bs and so on. And being a late convert to GMTs 5 and 6, I’m still rooting out the atavistic 4-style syntax (my scripts are “deprecable”!)

So, I will try to simplify one of my scripts to see if I can build up the capacity to make my complex plot by putting together more simple ones.

Thanks for the guidance on this!

Pat

No worries, so it is something else tripping up since the grid itself clearly works. GMT 6 introduced new things like subplots and planes and insets and stuff which is brave new world but so much easier. Post a failing script using that grid if you get there.

OK I’ll work on that tomorrow.

But one correction: I use csh scripts, not bash like I said. I think you guys converted to bash in the late 90s or early 00s or something like that and I never followed because of inertia and the csh tricks I constantly used. And I think they are sufficiently different to cause issues (I would be delighted to be wrong).

And yes, I’m reading the GMT 6 documentation and it does seem very “brave new world”. We’ll see how it goes tomorrow. Thanks again!

Update: I’ve managed to re-create all three individual sub-plots of my figures using the new (to me) GMT 6 syntax. Wow, Smooth! I for one am very glad to see -K and -O banished to the fiery hell in which they belong!

I am now working on grokking the “subplot” module, so that I can re-create the integrated figure. I hope to accomplish this tomorrow.

I think I’m going to like this particular brave new world!

1 Like