Conversion of GMT grd files to AAIGrid files reverses north and south

[hellyj]

I'm doing some grid format conversions and ran into an odd result. I think it's a problem with gdal_translate but want to mention it here in case anyone else has seen it before and can verify the problem but also to warn anyone who is doing GMT->AAIGrid transformations. Looks to me like gdal_translate is reading the *.grd file in the wrong order and inverting north and south. This doesn't happen when I use gdal_translate with geotiffs so I think it's localized to the GMT driver that it is using. When I have time I will check the code and see if I can find the problem. I've confirmed that the *.grd file is fine. I can convert it to an Arc grid using the gmt2arc script from Dawn Wright at OSU.

[warmerdam]

Hi,

Could you provide a sample input file, and the command you used to translate it?

Ivan,

Could you look into this? I'm presuming it is an orientation issue in the GMT driver, and it should show in a gdalinfo report on the input file.

[ilucena]

Yes, but I am going to need data samples.

[ilucena]

I download some GMT/grd sample files and I could find any problem with the conversion.

% gdal_translate -of aaigrid tbase15m.grd out.asc
Input file size is 1441, 721
0...10...20...30...40...50...60...70...80...90...100 - done.

from osgeo import gdal
a = gdal.Open("tbase15m.grd")
b = gdal.Open("out.asc")
a.GetGeoTransform()
(-180.125, 0.25, 0.0, 90.125, 0.0, -0.25)
b.GetGeoTransform()
(-180.125, 0.25, 0.0, 90.125, 0.0, -0.25)

I load it on QGIS and everything fits well. See screenshot on attachment.

[jluis]

I'm sorry but I have a file that reproduces this behavior

gdal_translate -of GTiff mola.grd lixo.tiff

When I open both "mola.grd" and "lixo.tiff" in Mirone they are upside down one relative to the other.

Since the file is to big to attach, I placed it temporarily at w3.ualg.pt/~jluis/ftp/tmp/mola.grd.zip

Joaquim

[ilucena]

Joaquim,

Thank you for providing that file. See what I got:

% gdalinfo mola.grd -checksum Driver: netCDF/Network Common Data Format Files: mola.grd Size is 1281, 1281 Coordinate System is `' Origin = (-0.007807530000000,-15.000000000000000) Pixel Size = (0.007812500000000,-0.007812500000000) Metadata:

NC_GLOBAL#Conventions=COARDS/CF-1.0 NC_GLOBAL#title=Grid imported trhough GDAL NC_GLOBAL#history=grdmath mola.nc=gd FLIPUD = mola.grd=ns NC_GLOBAL#GMT_version=4.5.3 (CVS Mar 20 2010 04:10:17) z#long_name=z z#_FillValue=-32768 z#actual_range=-1862, 2379 x#long_name=x x#actual_range=-0.00390128, 9.9961 y#long_name=y y#actual_range=-25.0039, -15.0039

Corner Coordinates: Upper Left ( -0.0078075, -15.0000000) Lower Left ( -0.0078075, -25.0078125) Upper Right ( 10.0000050, -15.0000000) Lower Right ( 10.0000050, -25.0078125) Center ( 4.9960987, -20.0039062) Band 1 Block=1281x1 Type=Int16, ColorInterp?=Undefined

Checksum=48852 NoData? Value=-32768 Metadata:

NETCDF_VARNAME=z

% gdal_translate -of aaigrid mola.grd out.asc
Input file size is 1281, 1281
0...10...20...30...40...50...60...70...80...90...100 - done.

% gdalinfo out.asc -checksum
Driver: AAIGrid/Arc/Info ASCII Grid
Files: out.asc
       out.asc.aux.xml
Size is 1281, 1281
Coordinate System is `'
Origin = (-0.007807530000000,-15.000000000000000)
Pixel Size = (0.007812500000000,-0.007812500000000)
Metadata:
  NC_GLOBAL#Conventions=COARDS/CF-1.0
  NC_GLOBAL#title=Grid imported trhough GDAL
  NC_GLOBAL#history=grdmath mola.nc=gd FLIPUD = mola.grd=ns
  NC_GLOBAL#GMT_version=4.5.3 (CVS Mar 20 2010 04:10:17)
  z#long_name=z
  z#_FillValue=-32768
  z#actual_range=-1862, 2379
  x#long_name=x
  x#actual_range=-0.00390128, 9.9961
  y#long_name=y
  y#actual_range=-25.0039, -15.0039
Corner Coordinates:
Upper Left  (  -0.0078075, -15.0000000)
Lower Left  (  -0.0078075, -25.0078125)
Upper Right (  10.0000050, -15.0000000)
Lower Right (  10.0000050, -25.0078125)
Center      (   4.9960987, -20.0039062)
Band 1 Block=1281x1 Type=Int32, ColorInterp=Undefined
  Checksum=48852
  NoData Value=-32768
  Metadata:
    NETCDF_VARNAME=z

The image looks the same, cheksum are the same.

Can you please discribe what versions/platform you are using?

Obrigado.

[jluis]

Olá Ivan,

I'm using a (recent) trunk version on Win 7 and Mirone, but you can see the UD flipping issue in GMT as well

gdal_translate -of aaigrid mola.grd out.asc

# Create a GMT palette grd2cpt mola.grd > pal.cpt

# Make an image of the netCDF grid grdimage mola.grd -Cpal.cpt -JX12c -P > grd.ps

# Convert the aai grid to netCDF with GMT (GMT5 will be able to do this directly) xyz2grd out.asc -E -Gaai.grd

# Make an image out of it grdimage aai.grd -Cpal.cpt -JX12c -P > aai.ps

Now grd.ps and aai.ps are flipped one relatively to the other

(Much easier to use Mirone an see the same)

Joaquim

[kyle]

Can someone point me to or post a small dataset that is being flipped?

[ilucena]

I can't. I tried with the sample provided by Joaquim but on my computer I can't see it flipped.

w3.ualg.pt/~jluis/ftp/tmp/mola.grd.zip

There is a small one here too:

http://trac.osgeo.org/gdal/ticket/796

[jluis]

Ivan

That "small one" is from a 5 years old ticket. This problem has lasted for many years but was apparently fixed recently (1.7.1 or close) and I thought it was solved for good, but I than found this case where it still shows up. I don't understand why it happens but the fact is that I'm able to reproduce with either Mirone or GMT. When you say that you are not able to see it on your computer you are referring to the gdalinfo reports only or you did display them in another program? Did you try my GMT recipe?

Joaquim

[ilucena]

Joaquim, Yes, I can load that image on QGIS. But if there is something wrong with GDAL that is not going to be conclusive, right. So I am checking out GMT now from CVS and I am going to build/install it and I will get back to you as soon as I can.

[ilucena]

Joaquim, I finally got GMT to run on Linux. So, yes, I can see the resulting images flipped after been converted to Postscript. I didn't see that problem with tbase15m.grd tough.

[ilucena]

I have some clues that I want to share so maybe you guys can make some of it.

On my system, the file mola.grd is been supported by the netCDF driver, not the GMT driver.

The code on netCDF driver cannot figure out that the image is bottom-up:

http://trac.osgeo.org/gdal/browser/trunk/gdal/frmts/netcdf/netcdfdataset.cpp#L1426

So, it reads it top-down:

http://trac.osgeo.org/gdal/browser/trunk/gdal/frmts/netcdf/netcdfdataset.cpp#L501

I tried to skip the netCDF driver to see if the GMT driver would take it but it didn't. See the before-and-after gdalinfo report:

% gdalinfo mola.grd
Driver: netCDF/Network Common Data Format
...
% export GDAL_SKIP=netCDF
% gdalinfo mola.grd
ERROR 4: `mola.grd' not recognised as a supported file format.

So now, looking at the GMT driver code, I noticed that the Open method recognize the file at first but then it give up (desiste[pt]) because it cannot find the dimension field:

http://trac.osgeo.org/gdal/browser/trunk/gdal/frmts/netcdf/gmtdataset.cpp#L226

If I change the ifdef to ifndef, just to turn that warming message on, here is what I get:

% gdalinfo mola.grd
Warning 1: mola.grd is a GMT file, but not in GMT configuration.
gdalinfo failed - unable to open 'mola.grd'.

That leads me to believe that or the file is corrupted or the dimension field is stored elsewhere, just like it happens all the time with HDF files.

So, I try to use some of the GMT command line tools to get information about that file. That is what I've got so far:

% gmtinfo mola.grd
gmtinfo : Cannot find leg mola.grd
usage: gmtinfo leg(s)

% grdinfo mola.grd
mola.grd: Title: Grid imported trhough GDAL
mola.grd: Command: grdmath mola.nc=gd FLIPUD = mola.grd=ns
mola.grd: Remark:
mola.grd: Gridline node registration used
mola.grd: Grid file format: ns (# 16) GMT netCDF format (short)  (COARDS-compliant)
mola.grd: x_min: -0.00390128 x_max: 9.99609872 x_inc: 0.0078125 name: x nx: 1281
mola.grd: y_min: -25.00390625 y_max: -15.00390625 y_inc: 0.0078125 name: y ny: 1281
mola.grd: z_min: -1862 z_max: 2379 name: z
mola.grd: scale_factor: 1 add_offset: 0

I don't know much about that stuff, but it seems to me that this file is not GMT and it was created as GRD with a FLIPUP parameter on a command named grdmath.

Does explain everything? Does it means that whatever is happen on GMT and netCDF drivers is coerrent with what is on that file, or not?

[jluis]

Hi Ivan,

I am on a mobile connection and won't be able to follow this closely until next week, but I have some information that could be useful for you right now so here it goes.

A bit of story. The GDAL netCDF driver had troubles during many years in not recognizing if a netcdf file was written top-down or bottom up. Since it was me who wrote the GMT GDAL wraper I hard-coded a test to always UD flip netCDF grids imported through GDAL. That is why you see that FLIPUD on argument to GMT's grdmath program. That test is now removed as things seamed to work well since GDAL 1.7.1(?)

GMT test on the way the netCDF file is written, so we don't have that sort of problems. Regarding the GMT driver on GDAL. We cannot use it for reading because it's for the old version of GMT netCDF grids (non COARDS) that used a different construct on the data array. Now, GMT does not have a netCDF "format" of its own. All it does is to write/read a COARDS compliant grid.

Since you compiled the GMT CVS version (otherwise the flipud test had not been removed yet) you can do it to have GDAL support (there is a compile switch to ling against GDAL). You will than access netCDF via GDAL by appending '=gd' to the grid's name. Like in

"grdmath mola.nc=gd FLIPUD = mola.grd"

that you saw in the comment's section of the mola.grd grid.

But to the point. If GDAL now recognizes when a grid is written top-down or bottom-up, why does it fail in current case?

[ilucena]

Joaquim,

If I understand you correctly the file mola.grd should not be read by the GMT driver.

Looking at the netCDF driver I noticed that there is bug related to that variable that decide about flipping the image, the szDimNameY:

1425	                /* for CF-1 conventions, assume bottom first */
1426	                if( EQUAL( szDimNameY, "lat" ) && pdfYCoord[0] < pdfYCoord[1] )
1427	                    poDS->bBottomUp = TRUE;

Since it was not initialized properly, the content was "y t". Originally it was "y type=int" and then it has terminated at position [3] by that code here:

783	    for( i = 0; (i < strlen( poDS->papszDimName[ poDS->nDimXid ] )  && 
784	                 i < 3 ); i++ ) {
785	        szDimNameX[i] = tolower( ( poDS->papszDimName[poDS->nDimXid] )[i] );
786	    }
787	    szDimNameX[3] = '\0';
788	    for( i = 0; (i < strlen( poDS->papszDimName[ poDS->nDimYid ] )  && 
789	                 i < 3 ); i++ ) {
790	        szDimNameY[i] = tolower( ( poDS->papszDimName[poDS->nDimYid] )[i] );
791	    }
792	    szDimNameY[3] = '\0';

Fixing that issue didn't help us, unfortunately.

In ReadBlock? "bBottomUp" is still FALSE at reading:

499	    start[nBandXPos] = 0;          // x dim can move arround in array
500	    // check y order
501	    if( ( ( netCDFDataset *) poDS )->bBottomUp ) {
502	        start[nBandYPos] = ( ( netCDFDataset * ) poDS )->ydim - 1 - nBlockYOff;
503	    } else {
504	        start[nBandYPos] = nBlockYOff; // y
505	    }

So it is reading the image from TOP down TO the BOTTON. Regards.

[ilucena]

Please try r20006

Here is the major change:

-          if( EQUAL( szDimNameY, "lat" ) && pdfYCoord[0] < pdfYCoord[1] )
+          if( ( EQUAL( szDimNameY, "lat" ) || EQUAL( szDimNameY, "y" ) )
+                  && pdfYCoord[0] < pdfYCoord[1] )
               poDS->bBottomUp = TRUE;

It assumes that if a maximum Y is smaller than a minimum Y the image is flipped.

Please let me know if there is any objection and reopen the ticket.