Feature #254

GDAL driver imports gridline registered instead of pixel registered

Added by Florian over 5 years ago. Updated over 5 years ago.

Status:NewStart date:2013-04-25
Priority:NormalDue date:
Assignee:-% Done:

0%

Category:-
Target version:-
Platform:

Description

GDAL only handles pixel registered grids with GTiff being an exception. Hence, GMT has to convert to pixel registration when writing through GDAL. But when reading through GDAL, grids should stay pixel registered and not be converted to gridline registration. Since GDAL's Gtiff driver supports both modes (http://trac.osgeo.org/gdal/wiki/rfc33_gtiff_pixelispoint), GMT should set AREA_OR_POINT properly and not convert to pixel registration.

grdmath -R0/10/0/10 -I1 1 = foo.tif=gd:gtiff
grdmath -r -R0/10/0/10 -I1 1 = fooP.tif=gd:gtiff

grdinfo foo.tif fooP.tif
foo.tif: Title: Grid imported via GDAL
foo.tif: Command: 
foo.tif: Remark: 
foo.tif: Gridline node registration used [Cartesian grid]
foo.tif: Grid file format: gd = Import/export through GDAL
foo.tif: x_min: 0 x_max: 10 x_inc: 1 name: x nx: 11
foo.tif: y_min: 0 y_max: 10 y_inc: 1 name: y ny: 11
foo.tif: z_min: 1 z_max: 1 name: z
foo.tif: scale_factor: 1 add_offset: 0
fooP.tif: Title: Grid imported via GDAL
fooP.tif: Command: 
fooP.tif: Remark: 
fooP.tif: Gridline node registration used [Cartesian grid]
fooP.tif: Grid file format: gd = Import/export through GDAL
fooP.tif: x_min: 0.5 x_max: 9.5 x_inc: 1 name: x nx: 10
fooP.tif: y_min: 0.5 y_max: 9.5 y_inc: 1 name: y ny: 10
fooP.tif: z_min: 1 z_max: 1 name: z
fooP.tif: scale_factor: 1 add_offset: 0

gdalinfo foo.tif
Driver: GTiff/GeoTIFF
Files: foo.tif
       foo.tif.aux.xml
Size is 11, 11
Coordinate System is `'
Origin = (-0.500000000000000,10.500000000000000)
Pixel Size = (1.000000000000000,-1.000000000000000)
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  -0.5000000,  10.5000000) 
Lower Left  (  -0.5000000,  -0.5000000) 
Upper Right (  10.5000000,  10.5000000) 
Lower Right (  10.5000000,  -0.5000000) 
Center      (   5.0000000,   5.0000000) 
Band 1 Block=11x11 Type=Float32, ColorInterp=Gray
  NoData Value=nan

gdalinfo fooP.tif
Driver: GTiff/GeoTIFF
Files: fooP.tif
       fooP.tif.aux.xml
Size is 10, 10
Coordinate System is `'
Origin = (0.000000000000000,10.000000000000000)
Pixel Size = (1.000000000000000,-1.000000000000000)
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (   0.0000000,  10.0000000) 
Lower Left  (   0.0000000,   0.0000000) 
Upper Right (  10.0000000,  10.0000000) 
Lower Right (  10.0000000,   0.0000000) 
Center      (   5.0000000,   5.0000000) 
Band 1 Block=10x10 Type=Float32, ColorInterp=Gray
  NoData Value=nan

gdal.tif (642 Bytes) Florian, 2013-05-01 21:30

area.tif (807 Bytes) Florian, 2013-05-01 21:30

point.tif (807 Bytes) Florian, 2013-05-01 21:30


Related issues

Related to Bug #403: GDAL import of small grid fails Closed 2013-10-09

Associated revisions

Revision 13222
Added by Joaquim over 4 years ago

Check the Metadata content to search for the GeoTiff tag "AREA_OR_POINT". If found and equals to 'Area' for a grid, than accept it to mean that the grid was original in pixel registration.
This partially addresses #254 and #541

History

#1 Updated by Florian over 5 years ago

Note also that grdinfo lies: Both grids are of course pixel registered.

#2 Updated by Joaquim over 5 years ago

I agree that when writing geotiff grids the original registration should be preserved via the AREA_OR_POINT property. At the time of the driver writing that property was being ignored and in practice all geotiffs were assumed pixel registration.
The other way around is much less clear. Saying that grdinfo is lying is a question of perspective. It is telling the true from the perspective of how the data will be represented internally by GMT. And it is not lying anymore than gdalinfo does when reports the equivalent info. For example SURFER grids do not even support the concept of pixel registration, however gdalinfo reports as if they were.

grdmath -R0/10/0/10 -I1 1 = lixo.grd=gd:GSBG
gdalinfo lixo.grd
Driver: GSBG/Golden Software Binary Grid (.grd)
Files: lixo.grd
Size is 11, 11
Coordinate System is `'
Origin = (-0.500000000000000,10.500000000000000)
Pixel Size = (1.000000000000000,-1.000000000000000)
Corner Coordinates:
Upper Left ( -0.5000000, 10.5000000)
Lower Left ( -0.5000000, -0.5000000)
Upper Right ( 10.5000000, 10.5000000)
Lower Right ( 10.5000000, -0.5000000)
Center ( 5.0000000, 5.0000000)
Band 1 Block=11x1 Type=Float32, ColorInterp=Undefined
Min=1.000 Max=1.000
NoData Value=1.701410009187828e+038

while GMT

grdinfo lixo.grd
lixo.grd: Title: Grid originally in Surfer 6 format
lixo.grd: Command:
lixo.grd: Remark:
lixo.grd: Gridline node registration used [Cartesian grid]
lixo.grd: Grid file format: sf = Golden Software Surfer format 6 (32-bit float)
lixo.grd: x_min: 0 x_max: 10 x_inc: 1 name: x nx: 11
lixo.grd: y_min: 0 y_max: 10 y_inc: 1 name: y ny: 11
lixo.grd: z_min: 1 z_max: 1 name: z
lixo.grd: scale_factor: 1 add_offset: 0
The sad thing in GDAL is that it doesn't provide a mechanism to inform 3rd parties the original registration type of the file accessed via its lib, so one has to make a choice.
  1. 1) accept that everything (but geotiff and eventually netcdf) that comes from GDAL is pixel registered
  2. 2) convert all GRIDS to grid registration and for images let them maintain the pixel registration

We decided for 2) that also gives the liberty to adapt to cases when we clearly know that original file was pixel registered (e.g. using the AREA_OR_POINT mechanism)

#3 Updated by Joaquim over 5 years ago

  • Tracker changed from Bug to Feature

Tried to implement the AREA_OR_POINT when writing geotiff grids but, ... case closed

grdmath -R0/10/0/10 -I1 1 = lixo.tiff=gd:gtiff
Warning 6: Driver GTiff does not support AREA_OR_POINT creation option

BTW, I also implemented a case insensitivity for the driver name. Before 'gtiff' was ignored and and fall to the default driver, which is ... GTiff

There is still the reading part but maybe too much work for what it worth, well at least for now.

#4 Updated by Florian over 5 years ago

Tried to implement the AREA_OR_POINT when writing geotiff grids but, ... case closed

grdmath -R0/10/0/10 -I1 1 = lixo.tiff=gd:gtiff
Warning 6: Driver GTiff does not support AREA_OR_POINT creation option

Strange because this works:

grdmath -R0/10/0/10 -I1 1 = foo.nc
gdal_translate -mo AREA_OR_POINT=Point -of GTiff foo.nc foo.tif
gdalinfo --version
GDAL 1.9.2, released 2012/10/08

There is still the reading part but maybe too much work for what it worth, well at least for now.

Yes, it would really be convenient to have grids converted to true pixel registration when writing and read with pixel registration.

#5 Updated by Joaquim over 5 years ago

Strange because this works:
[...]

GDAL had two creation methods: CreateCopy() and Create(). Quite likely is gdal_translate is using the one that allows the AREA_OR_POINT creation option, but we are using the MEM driver (since we are not reading from a file) and the method that it uses (never remember which) does not allow that option.

There is still the reading part but maybe too much work for what it worth, well at least for now.

Yes, it would really be convenient to have grids converted to true pixel registration when writing and read with pixel registration.

Note that it would only matter for the (rare) GeoTiffs with AREA_OR_POINT=Point set. All the other grids have lost their original registration info and had become 'pixel by force' registered.

#6 Updated by Florian over 5 years ago

GDAL had two creation methods: CreateCopy() and Create(). Quite likely is gdal_translate is using the one that allows the AREA_OR_POINT creation option, but we are using the MEM driver (since we are not reading from a file) and the method that it uses (never remember which) does not allow that option.

How sad...

There is still the reading part but maybe too much work for what it worth, well at least for now.

Yes, it would really be convenient to have grids converted to true pixel registration when writing and read with pixel registration.

Note that it would only matter for the (rare) GeoTiffs with AREA_OR_POINT=Point set. All the other grids have lost their original registration info and had become 'pixel by force' registered.

Sure, but since AREA_OR_POINT=Point is not working and 99.9% of all GDAL readable grids are pixel registered it would be nice if GMT also treats all GDAL imported grids as pixel internally per default. I.e., the behavior should be identical to loading a pixel-registered netcdf grid file.

#7 Updated by Joaquim over 5 years ago

How sad...

Well, you can always fill a feature request to see if that is possible. I've seen additions to Create() method (if that's the one in question) when their nonexistence was not based on technical reasons.

Sure, but since AREA_OR_POINT=Point is not working and 99.9% of all GDAL readable grids are pixel registered

No they are not. netCDF grids for instance are all reported to be pixel registered and we both know than much more than 1% netCDF grids trailing around in the world are grid registered.

#8 Updated by Florian over 5 years ago

Sure, but since AREA_OR_POINT=Point is not working and 99.9% of all GDAL readable grids are pixel registered

No they are not. netCDF grids for instance are all reported to be pixel registered and we both know than much more than 1% netCDF grids trailing around in the world are grid registered.

True, but why would you import netcdf via GDAL? This defeats the purpose. We should simply refuse to do so. GDAL import is supposed to read GMT unsupported grids, isn't it? Can you probe the grid format with GDAL prior to reading? Otherwise call GMT_grd_get_format?

#9 Updated by Joaquim over 5 years ago

I don't agree that we should try to impose any kind of censuring on what formats are read via GDAL. That's a pure user choice. And it's not so senseless to read netCDF grids via GDAL. Some have projection reference information stored as attributes and we don't look for that in the GMT functions, but we import it when read via GDAL.
Ah, and BTW the SRTM tiles (before being re-interpolated to please the pixel mania) were also grid registered.

#10 Updated by Joaquim over 5 years ago

Strange because this works:
[...]

Actually you found a GDAL bug. This is what gdalinfo now reports

Corner Coordinates:
Upper Left ( 0.0000000, 10.0000000)
Lower Left ( 0.0000000, -1.0000000)
Upper Right ( 11.0000000, 10.0000000)
Lower Right ( 11.0000000, -1.0000000)
Center ( 5.5000000, 4.5000000)
Band 1 Block=11x11 Type=Float32, ColorInterp=Gray

#11 Updated by Joaquim over 5 years ago

r11436 has an attempt (but #if 0 protected) to identify AREA_OR_POINT=Point cases, but I lack an example file to test it. The gdal_translate example above not only creates a bugged tiff as contains no AREA_OR_POINT metadata (checked with an Hex editor and the code in r11436)

#12 Updated by Florian over 5 years ago

You are right, gdal_translate does not set AREA_OR_POINT at all!

To add more confusion, I attach several files:
gdal.tif created by gdal_translate from ascii
area.tif created by libgeotiff (RasterPixelIsArea)
point.tif created by libgeotiff (RasterPixelIsPoint)

In the following note that:
  1. gdalinfo correctly interprets AREA_OR_POINT
  2. AREA_OR_POINT vanishes with gdal_translate or grdreformat
  3. grdinfo reports all grids written by grdreformat (float?) as gridline registered whereas all other tiffs are reported as pixel registered
  4. grdinfo reports grids written by grdreformat (byte?) as pixel registered

Especially the last two items are very strange.

grdreformat area.tif area_gmt.tif=gd:gtiff
grdreformat point.tif point_gmt.tif=gd:gtiff
grdreformat gdal.tif gmt.tif=gd:gtiff
grdreformat gdal.tif gmtu8.tif=gd:gtiff/u8

gdalinfo gdal.tif
Driver: GTiff/GeoTIFF
Files: gdal.tif
Size is 20, 20
Coordinate System is `'
Origin = (130.000000000000000,32.000000000000000)
Pixel Size = (1.000000000000000,-1.000000000000000)
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( 130.0000000,  32.0000000) 
Lower Left  ( 130.0000000,  12.0000000) 
Upper Right ( 150.0000000,  32.0000000) 
Lower Right ( 150.0000000,  12.0000000) 
Center      ( 140.0000000,  22.0000000) 
Band 1 Block=20x20 Type=Byte, ColorInterp=Gray

gdalinfo area.tif 
Driver: GTiff/GeoTIFF
Files: area.tif
Size is 20, 20
Coordinate System is:
GEOGCS["Everest Ellipsoid Used.",
    DATUM["unknown",
        SPHEROID["Everest 1830 (1967 Definition)",6377298.556,300.8016999999915,
            AUTHORITY["EPSG","7016"]]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433]]
Origin = (130.000000000000000,32.000000000000000)
Pixel Size = (1.000000000000000,-1.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( 130.0000000,  32.0000000) (130d 0' 0.00"E, 32d 0' 0.00"N)
Lower Left  ( 130.0000000,  12.0000000) (130d 0' 0.00"E, 12d 0' 0.00"N)
Upper Right ( 150.0000000,  32.0000000) (150d 0' 0.00"E, 32d 0' 0.00"N)
Lower Right ( 150.0000000,  12.0000000) (150d 0' 0.00"E, 12d 0' 0.00"N)
Center      ( 140.0000000,  22.0000000) (140d 0' 0.00"E, 22d 0' 0.00"N)
Band 1 Block=20x20 Type=Byte, ColorInterp=Gray

gdalinfo point.tif 
Driver: GTiff/GeoTIFF
Files: point.tif
Size is 20, 20
Coordinate System is:
GEOGCS["Everest Ellipsoid Used.",
    DATUM["unknown",
        SPHEROID["Everest 1830 (1967 Definition)",6377298.556,300.8016999999915,
            AUTHORITY["EPSG","7016"]]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433]]
Origin = (129.500000000000000,32.500000000000000)
Pixel Size = (1.000000000000000,-1.000000000000000)
Metadata:
  AREA_OR_POINT=Point
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( 129.5000000,  32.5000000) (129d30' 0.00"E, 32d30' 0.00"N)
Lower Left  ( 129.5000000,  12.5000000) (129d30' 0.00"E, 12d30' 0.00"N)
Upper Right ( 149.5000000,  32.5000000) (149d30' 0.00"E, 32d30' 0.00"N)
Lower Right ( 149.5000000,  12.5000000) (149d30' 0.00"E, 12d30' 0.00"N)
Center      ( 139.5000000,  22.5000000) (139d30' 0.00"E, 22d30' 0.00"N)
Band 1 Block=20x20 Type=Byte, ColorInterp=Gray

grdinfo *.tif
gdal.tif: Title: Grid imported via GDAL
gdal.tif: Command: 
gdal.tif: Remark: 
gdal.tif: Pixel node registration used [Cartesian grid]
gdal.tif: Grid file format: gd = Import/export through GDAL
gdal.tif: x_min: 130 x_max: 150 x_inc: 1 name: x nx: 20
gdal.tif: y_min: 12 y_max: 32 y_inc: 1 name: y ny: 20
gdal.tif: z_min: 0 z_max: 0 name: z
gdal.tif: scale_factor: 1 add_offset: 0

gmt.tif: Title: Grid imported via GDAL
gmt.tif: Command: 
gmt.tif: Remark: 
gmt.tif: Gridline node registration used [Cartesian grid]
gmt.tif: Grid file format: gd = Import/export through GDAL
gmt.tif: x_min: 130.5 x_max: 149.5 x_inc: 1 name: x nx: 20
gmt.tif: y_min: 12.5 y_max: 31.5 y_inc: 1 name: y ny: 20
gmt.tif: z_min: 0 z_max: 0 name: z
gmt.tif: scale_factor: 1 add_offset: 0

gmtu8.tif: Title: Grid imported via GDAL
gmtu8.tif: Command: 
gmtu8.tif: Remark: 
gmtu8.tif: Pixel node registration used [Cartesian grid]
gmtu8.tif: Grid file format: gd = Import/export through GDAL
gmtu8.tif: x_min: 130 x_max: 150 x_inc: 1 name: x nx: 20
gmtu8.tif: y_min: 12 y_max: 32 y_inc: 1 name: y ny: 20
gmtu8.tif: z_min: 0 z_max: 0 name: z
gmtu8.tif: scale_factor: 1 add_offset: 0

area.tif: Title: Grid imported via GDAL
area.tif: Command: 
area.tif: Remark: 
area.tif: Pixel node registration used [Cartesian grid]
area.tif: Grid file format: gd = Import/export through GDAL
area.tif: x_min: 130 x_max: 150 x_inc: 1 name: x nx: 20
area.tif: y_min: 12 y_max: 32 y_inc: 1 name: y ny: 20
area.tif: z_min: 0 z_max: 0 name: z
area.tif: scale_factor: 1 add_offset: 0

area_gmt.tif: Title: Grid imported via GDAL
area_gmt.tif: Command: 
area_gmt.tif: Remark: 
area_gmt.tif: Gridline node registration used [Cartesian grid]
area_gmt.tif: Grid file format: gd = Import/export through GDAL
area_gmt.tif: x_min: 130.5 x_max: 149.5 x_inc: 1 name: x nx: 20
area_gmt.tif: y_min: 12.5 y_max: 31.5 y_inc: 1 name: y ny: 20
area_gmt.tif: z_min: 0 z_max: 0 name: z
area_gmt.tif: scale_factor: 1 add_offset: 0

point.tif: Title: Grid imported via GDAL
point.tif: Command: 
point.tif: Remark: 
point.tif: Pixel node registration used [Cartesian grid]
point.tif: Grid file format: gd = Import/export through GDAL
point.tif: x_min: 129.5 x_max: 149.5 x_inc: 1 name: x nx: 20
point.tif: y_min: 12.5 y_max: 32.5 y_inc: 1 name: y ny: 20
point.tif: z_min: 0 z_max: 0 name: z
point.tif: scale_factor: 1 add_offset: 0

point_gmt.tif: Title: Grid imported via GDAL
point_gmt.tif: Command: 
point_gmt.tif: Remark: 
point_gmt.tif: Gridline node registration used [Cartesian grid]
point_gmt.tif: Grid file format: gd = Import/export through GDAL
point_gmt.tif: x_min: 130 x_max: 149 x_inc: 1 name: x nx: 20
point_gmt.tif: y_min: 13 y_max: 32 y_inc: 1 name: y ny: 20
point_gmt.tif: z_min: 0 z_max: 0 name: z
point_gmt.tif: scale_factor: 1 add_offset: 0

#13 Updated by Joaquim over 5 years ago

Florian, my GMT time is currently booked to the gmtmex work. I'll look into this later.

Also available in: Atom PDF