Opened 13 years ago

Closed 12 years ago

#1462 closed defect (fixed)

GMT DRIVER out-by-1 problem

Reported by: petermorgan@… Owned by: warmerdam
Priority: normal Milestone: 1.4.4
Component: GDAL_Raster Version: 1.4.0
Severity: normal Keywords: gmt
Cc: warmerdam

Description (last modified by Mateusz Łoskot)

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

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.

Change History (6)

comment:4 Changed 13 years ago by warmerdam

Cc: warmerdam added
Description: modified (diff)
Milestone: 1.4.2
Owner: changed from warmerdam to Mateusz Łoskot
Priority: highestnormal
Severity: blockernormal

comment:5 Changed 13 years ago by Mateusz Łoskot

Description: modified (diff)

comment:6 Changed 13 years ago by warmerdam

Component: defaultGDAL_Raster
Keywords: gmt added
Milestone: 1.4.21.4.3

Put off till 1.4.3.

comment:7 Changed 12 years ago by warmerdam

Owner: changed from Mateusz Łoskot to warmerdam

I'll take in order to try and polish off for 1.4.3.

comment:8 Changed 12 years ago by warmerdam

Milestone: 1.4.31.4.4

This will not be addressed in time for 1.4.3.

comment:9 Changed 12 years ago by warmerdam

Resolution: fixed
Status: newclosed

Peter,

I appologise for the lengthy delay in dealing with this issue. After reading through your notes carefully, and going through the driver I have discovered that the core of the problem is that the "CreateCopy?" implementation for gmtdataset.cpp was writing "node_offset=0" (which essentially means origin at pixel center) but putting out bounds values were based on the top left corner.

So any file generated by GDAL would end up with the wrong bounds and pixel size.

I have corrected this in trunk (r12890) and 1.4-branch (r12891). The fix will be in the GDAL 1.4.4 release, likely within 2-5 weeks.

I have also added a very modest test to the GDAL test suite (r12893).

Hopefully GDAL's GMT support will be smoother in the future.

Note: See TracTickets for help on using tickets.