GMT DRIVER out-by-1 problem

[petermorgan@…]

The out-by-1 problem  with the GMT driver was again recently reported. It appears that this is quite an old problem as the initial postings to the GDAL community were made in 2005. The first appears to be  #796 on 13 March 2004 by jluis.  The second was on 15 April 2005  also by jluis. Another on November 18 2005, Rich Signal  reports the same problem.

We have discussed the problem in the GMT community.  GMT is able to support grid referenced data AND pixel registered data. We have routines to transfer between these if there is a need to do so. The use of grid referenced data in GMT is mainly historical so resampling can be very lossy especially at scales near the Nyquist frequency of the grid. Thus the solution to the problem may be to output the region of interest as a pixel registered grid and to make the header information consistent with this option. What I am attempting to say here is that the GMT community is not locked into a grid reference model when the underlying data set is pixel registered. grdinfo in GMT will return the type of grid that is being used so the onus can be placed on the GMT community to understand the imported grid.



For these tests I cut from the USGS Seamless server the following region using the direct entry of the region of interest.

                                    |                                          |
                                    |                                          |
                     ----------------------------------------------------------  05 01 00.000    (5.01666667)
                                    |                                          |
                                    |                                          |
                                    |                                          |                                      |                                         |
                 ------------------------------------------------------------   03 59 00.000   (3.983333333)
                                    |                                          |
                                    |                                          |
                            113 59 00.000                     115 01 00.000
                         (113.98333333)                     (115.01666667)


Using gdalinfo on the returned file  we see the following
[peterm@currawong 16733119]$ gdalinfo w001001.adf
Driver: AIG/Arc/Info Binary Grid
Size is 1240, 1240
Coordinate System is:
GEOGCS["WGS 84",
   DATUM["WGS_1984",
       SPHEROID["WGS 84",6378137,298.257223563,
           AUTHORITY["EPSG","7030"]],
       TOWGS84[0,0,0,0,0,0,0],
       AUTHORITY["EPSG","6326"]],
   PRIMEM["Greenwich",0,
       AUTHORITY["EPSG","8901"]],
   UNIT["degree",0.0174532925199433,
       AUTHORITY["EPSG","9108"]],
   AXIS["Lat",NORTH],
   AXIS["Long",EAST],
   AUTHORITY["EPSG","4326"]]
Origin = (113.983333328038498,5.016666667786583)
Pixel Size = (0.000833333333300,-0.000833333333300)
Corner Coordinates:
Upper Left  ( 113.9833333,   5.0166667)
Lower Left  ( 113.9833333,   3.9833333)
Upper Right ( 115.0166667,   5.0166667)
Lower Right ( 115.0166667,   3.9833333)
Center      ( 114.5000000,   4.5000000)
Band 1 Block=512x4 Type=Float32, ColorInterp=Undefined
 Min=-7.000 Max=2371.000
 NoData Value=-3.40282346638529e+38

Note  that the subtle difference between the definition of the origin and the Upper Left corner definition is a numeric
problem. There are indeed 1240 cells, (columns) from 113.9833333 to 115.0166667 similarly there are 1240 cells, (rows), from
3.9833333 to 5.0166667. These cells are 0.000833333333300 degrees which is equivalent to 3 arc seconds.
This indicates that we do indeed have a pixel registered data set.

Running the  routine gdal_translate [peterm@currawong 16733119]$ gdal_translate w001001.adf -of GMT w001001.gmt
Input file size is 1240, 1240
[peterm@currawong 16733119]$
says that the input file size was  1240-by-1240
Now we can get information on the output file by two methods:
 1. using gdalinfo

   [peterm@currawong 16733119]$ gdalinfo w001001.gmt
Driver: GMT/GMT NetCDF Grid Format
Size is 1240, 1240
Coordinate System is `'
Origin = (113.982916325079145,5.017083670745943)
Pixel Size = (0.000834005918718,-0.000834005918718)
Corner Coordinates:
Upper Left  ( 113.9829163,   5.0170837)
Lower Left  ( 113.9829163,   3.9829163)
Upper Right ( 115.0170837,   5.0170837)
Lower Right ( 115.0170837,   3.9829163)
Center      ( 114.5000000,   4.5000000)
Band 1 Block=1240x1 Type=Float32, ColorInterp=Undefined

Notice that the cell size is no longer 0.00083333333 but .000834005918718.
This change in size is coupled with changes in the origin which is now
(113.982916325079145,5.017083670745943) which is
113 58 58.498    5 01 01.5012
The lower right corner is ( 115.0170837,   3.9829163)  or 115 01 01.501,  03 58 58.49868
The longitude difference of  1.034167 corresponds to 1241 cells when the size is 0.00083333
or                                                                                  1240 cells when the size is 0.00083400591871
A similar computation results for latitude.

This is a bastard system as it appears that the origins have been adjusted to move from pixel to grid registration but that is all. To make things consistent the cell width has be modified.

2. getting information with grdinfo a GMT command.
[peterm@currawong 16733119]$ grdinfo w001001.gmt
w001001.gmt: Title:
w001001.gmt: Command:
w001001.gmt: Remark:
w001001.gmt: Normal node registration used
w001001.gmt: grdfile format: cf (# 10)
w001001.gmt: x_min: 1.139833e+02 x_max: 1.150167e+02 x_inc: 8.333333e-04 name: meters nx: 1240
w001001.gmt: y_min: 3.983333e+00 y_max: 5.016667e+00 y_inc: 8.333333e-04 name: meters ny: 1240
w001001.gmt: z_min: -7.000000e+00 z_max: 2.371000e+03 name: meters
w001001.gmt: scale_factor: 1.000000e+00 add_offset: 0.000000e+00

Here the grid size is 0.0008333 and the number of grid cells is 1240-by-1240.
The origin of the grid is consistent with the pixel grid yet the grid is labelled as a normal default grid rather than a pixel grid.

I am a little puzzled that these two info routines return different information. I presume that this is because they interrogated different parts of the header.


It is hypothesised that the GDAL routines have failed to output the additional row and column that result from a conversion from a pixel registered grid to a grid reference grid. The correct process of transforming between these two grid types involves some form of resampling which is why the initial grid should be oversized if the process is pixed to reference grid.