I am trying to plot / manipulate a Bouguer gravity grid file from the GNS New Zealand database. I have tried both the original ascii file (attempting to convert it to a grid file) and a grd file of the same grid already processed by a colleague, but when I when grdimage I get the following error:
grdimage: Error: (x_max-x_min) must equal (NX + eps) * x_inc), where NX is an integer and |eps| <= 0.0001.
grdimage: Error: (y_max-y_min) must equal (NY + eps) * y_inc), where NY is an integer and |eps| <= 0.0001.
grdimage (gmt_grdio.c:1309(gmtlib_read_grd_info)): Use grdedit -A on your grid file to make region and increments compatible [Maps_and_Grids/NZbouguer.grd]
I am still able to plot a map, and everything looks okay, but I’m wondering what exactly this error means and how to go about fixing it. I believe it may have something to do with the fact that the original grid extent crosses the dateline and thus extends to longitudes of 190 degrees, which is throwing off GMT (though I am not actually plotting that part of the grid). I have tried the grdedit -A suggestion given in the error and it does not work. Any suggestions?
What it means is that the limits of your grid are not divisible into an integer number of grid increments (plus a small tolerance). For example, one cannot go from 1 to 5 at steps of 1.5. Here is the same x_max != x_min + N*x_inc But grdedit -A should had fixed it.
grdedit -A doesn’t do anything, even though it should. I continue to get the same error. Other than the -A flag is there some other flag I should be using? Or, are there specific arguments to the -A flag that need to be specified? From the documentation, it didn’t seem like it. I might be able to ignore it, but I think it’s causing some artefacts and aliasing in the grid when plotted, which is error that could potentially propagate in my gravity inversion, so I really would like to fix it. (I’m using gmt 5.4 btw).
Do we know for sure it is an equidistant grid? Or does it have variable spacing - in which case GMT generally does not support it.
The grids are from the New Zealand GNS database and were originally in Arc-ascii format. I converted them to netCDF format using grdconvert and specifying the ef format as the input format. I then cut them to my region of interest, which eliminates the points that pass the dateline (i.e. the ~190 degree lon points). Is there a way to tell whether it’s equidistant and/or make it equidistant if not? I would assume that such a widely used dataset for NZ is compatible with GMT.
Perhaps you can copy/past in the commands you ran. It is hard to know what went wrong from here. Perhaps your -R -I -r settings are not compatible or does not match the data file. And also how you cut the grid re dateline. All those steps could cause trouble.
Okay, so I’ve tried redoing things and discovered some stuff. When I convert the grid with grdconvert and don’t specify any format,
gmt grdconvert gns2016_NZ_freeair.asc gns2016_FA.nc -V
grdconvert: File gns2016_NZ_freeair.asc reads with GDAL driver AAIGrid
grdconvert: Translating file gns2016_NZ_freeair.asc (format gd = Import/export through GDAL) to file gns2016_FA.nc (format nf = GMT netCDF format (32-bit float), COARDS, CF-1.5)
grdconvert: File gns2016_NZ_freeair.asc reads with GDAL driver AAIGrid
gmt just uses gdal and it does not give an error. grdinfo on the resulting grid gives:
gns2016_FA.nc: Title: Produced by grdconvert
gns2016_FA.nc: Command: grdconvert gns2016_NZ_freeair.asc gns2016_FA.nc -V
gns2016_FA.nc: Remark:
gns2016_FA.nc: Gridline node registration used [Cartesian grid]
gns2016_FA.nc: Grid file format: nf = GMT netCDF format (32-bit float), COARDS, CF-1.5
gns2016_FA.nc: x_min: 160.0083335 x_max: 190.0089335 x_inc: 0.016667 name: x n_columns: 1801
gns2016_FA.nc: y_min: -59.9916665 y_max: -24.9909665 y_inc: 0.016667 name: y n_rows: 2101
gns2016_FA.nc: z_min: -252.645004272 z_max: 314.221008301 name: z
gns2016_FA.nc: scale_factor: 1 add_offset: 0
gns2016_FA.nc: format: netCDF-4 chunk_size: 129,132 shuffle: on deflation_level: 3
It says that the grid is cartesian, but it’s in lat lon, so it should say it’s in geographic (right?). Previously, I told grdconvert that the grid was geographic and it produced the error I showed in my original post. Now, grdimage still gives this error when I use -Jm, but not when I use -JX, which makes sense if it thinks the grid is cartesian, but when using -JX, everything gets stretched out all weird because the grid is actually geographic. Is there a way to get gmt to recognize it’s actually geographic without it throwing the error?
How can it know that those numbers represent lon, lat? Add -fg to your grdconvert command to inform GMT that the grid is geographic.
And looking at these numbers 160.0083335 x_max: 190.0089335 x_inc: 0.016667, I also fear that some truncation had occurred due to the use of ascii input data.
I can use -fg to tell gmt its geographic, but only then does it give the error about grdinfo: Error: (x_max-x_min) must equal (NX + eps) * x_inc), where NX is an integer and |eps| <= 0.0001. As you said, grdedit should fix this, but it does not. It seems like so long as the grid is left as cartesian, it doesn’t give the error, but if I tell it it’s geographic, then it does. I’ve tried grdedit on both the original esri ascii file and the converted grd file. Seems like there is one to many columns and rows in the grid; is there a way to manually cut the last column and last row of the grid? Could you please elaborate more on what you mean by “some truncation had occurred due to the use of ascii input data”?
You can try to convert with GMT instead of letting it go through GDAL with
grdconvert gns2016_NZ_freeair.asc=ef gns2016_FA.nc
but it shouldn’t make any difference.
I don’t understand why you want to strip out one row & column but see grdsample -T
What I mean by truncation is this. It’s not unusual that ASCII data does not contain enough decimal places to represent real numbers. For example your grid is at 1 min resolution, which is 1/60 = 0.016666666666666666 and x_min should be equal to
160 + 0.0166666666666666/2 = 160.00833333333333
but according to grdinfo it is 160.0083335
Can you make one those ascii files available? Preferably a small one that you can zip and post here.
I see; thanks for the description. I’ve uploaded two of the ascii files in the zip folder. Both are gravity grids.
Okay, it looks like the upload failed because even the smaller of the two files is pretty large. I’ll try again later, and keep working with the grids for now. Using =ef unfortunately does not work either.
Post a link to one of them
Those files have an error. The header says
xllcorner 160.000
yllcorner -60.000
but should say
xllcenter 160.000
yllcenter -60.000
because the grids are grid registered, not pixel registered.
and the cellsize should be more generous in 666. Make it
cellsize 0.01666666667
It worked! No more error. Thanks so much!
Maybe warn the data providers about the mistake?