[Gdal-dev] GMT & GDAL Geotiff grid registration
Rich Signell
rsignell at usgs.gov
Fri Nov 18 14:44:08 EST 2005
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
More information about the Gdal-dev
mailing list