[Gdal-dev] GMT & GDAL Geotiff grid registration

[Rich Signell]

GDAL-folk,

I think I'm seeing an issue with grid registration when I convert 
grids from GMT to Geotiff using GDAL:  GMT can be either PIXEL or 
LINE registered, but GDAL always produces PIXEL-registered images 
from GMT files.

The result is that the geotiff created from my GMT grid via:

   gdal_translate -a_srs EPSG:32618 rich.grd rich.tif

has the wrong grid spacing.

GMT's grdinfo reports the correct 5 m spacing:

[rsignell at ricsigdtlx seth]$ grdinfo rich.grd
rich.grd: Title: surface
rich.grd: Command: surface rich_5m.xyz 
-R/683500/684500/4565000/4566000/ -Grich.grd -C0.1 -I5 -T0.2 -V
rich.grd: Remark:
rich.grd: Normal node registration used
rich.grd: grdfile format # 0
rich.grd: x_min: 683500.0 x_max: 684500.0 x_inc: 5.0 units: 
user_x_unit nx: 201
rich.grd: y_min: 4565000.0 y_max: 4566000.0 y_inc: 5.0 units: 
user_y_unit ny: 201
rich.grd: z_min: 3.1 z_max: 17.0 units: user_z_unit
rich.grd: scale_factor: 1.0 add_offset: 0.0


while gdalinfo reports incorrect grid spacing, because it 
incorrectly assumed PIXEL==AREA registration:

gdalinfo rich.tif
[rsignell at ricsigdtlx seth]$ gdalinfo rich.tif
Driver: GTiff/GeoTIFF
Size is 201, 201
Coordinate System is:
PROJCS["WGS 84 / UTM zone 18N",
...
Origin = (683500.000000,4566000.000000)
Pixel Size = (4.97512438,-4.97512438)
Metadata:
   AREA_OR_POINT=Area
Corner Coordinates:
Upper Left  (  683500.000, 4566000.000)
Lower Left  (  683500.000, 4565000.000)
Upper Right (  684500.000, 4566000.000)
Lower Right (  684500.000, 4565000.000)
Center      (  684000.000, 4565500.000)
Band 1 Block=201x10 Type=Float32, ColorInterp=Gray


Perhaps I've missed a GDAL option somewhere, but I did at least 
try to look!

In case this is really a missing feature, the PIXEL vs. LINE 
registration of a GMT file is really easy to check.  You just read 
the "node_offset" attribute of the grid variable.

If node_offset=0, it's LINE
If node_offset=1, it's PIXEL

[rsignell at ricsigdtlx seth]$ ncdump -h rich.grd
netcdf rich {
dimensions:
         side = 2 ;
         xysize = 40401 ;
variables:
         double x_range(side) ;
                 x_range:units = "user_x_unit" ;
         double y_range(side) ;
                 y_range:units = "user_y_unit" ;
         double z_range(side) ;
                 z_range:units = "user_z_unit" ;
         double spacing(side) ;
         int dimension(side) ;
         float z(xysize) ;
                 z:scale_factor = 1. ;
                 z:add_offset = 0. ;
                 z:node_offset = 0 ;

// global attributes:
                 :title = "surface" ;
                 :source = "surface rich_5m.xyz 
-R/683500/684500/4565000/4566000/ -Grich.grd -C0.1 -I5 -T0.2 -V" ;
}

I'm working around this for now using the "-a_ullr" option on 
gdal_translate to specify different bounds:

gdal_translate -a_ullr 683497.5 4566002.5 684502.5 4564997.5 
-a_srs EPSG:32618 rich.grd rich3.tif

[rsignell at ricsigdtlx seth]$ gdalinfo rich3.tif
Driver: GTiff/GeoTIFF
Size is 201, 201
Coordinate System is:
...
Origin = (683497.500000,4566002.500000)
Pixel Size = (5.00000000,-5.00000000)
Metadata:
   AREA_OR_POINT=Area
Corner Coordinates:
Upper Left  (  683497.500, 4566002.500)
Lower Left  (  683497.500, 4564997.500)
Upper Right (  684502.500, 4566002.500)
Lower Right (  684502.500, 4564997.500)
Center      (  684000.000, 4565500.000)
Band 1 Block=201x10 Type=Float32, ColorInterp=Gray

but obviously it would be nice to have a *slightly* more elegant 
solution.


Thanks,
Rich

P.S. I'm using GDAL 1.3.1.0, FWTools 1.0.0a7, released 2005/11/10
-- 
Richard P. Signell
U.S. Geological Survey       Phone: (508) 457-2229
384 Woods Hole Road          Fax:   (508) 457-2310
Woods Hole, MA  02543-1598