[Gdal-dev] Mishandled "pixel is point" in geotiff

[Tim Beckmann]

Hi all (and Frank),

I sent an email to the list a few weeks ago about geotiff files with the 
"Pixel is point" option not being handled properly by gdal.  I only 
received a reply from one person off the list stating that it was a known 
problem. 

I did figure out how to pass the "-mo AREA_OR_POINT=POINT" parameter to 
gdal_translate to tell it to include the pixel is point indication in the 
output geotiff file.  So, my command line looks like:

        gdal_translate -mo AREA_OR_POINT=POINT -co INTERLEAVE=PIXEL \
                -a_srs "+proj=utm +zone=15 +datum=WGS84"  \
                -a_ullr 753600 3462900 801900 3260100 \
                'HDF4_SDS:UNKNOWN:"EO1H0220392005002110P0_HDF.L1T":0' 
testit_b001.tif

This results in the geotiff file properly indicating "pixel is point". 
However, the reported pixel size is incorrect.  It is being calculated 
incorrectly from the provided corners like it was still "pixel is area" 
(so instead of being 30 meters, it is just a little less than 30 meters). 
The hack at the bottom of the email works around it by rescaling the pixel 
size for the "pixel is point" case.  But it doesn't do anything for 
properly fixing the issue throughout gdal... (but it looks like the "pixel 
is point" support itself is a hack...)

Tim


diff -u -r1.157 geotiff.cpp
--- geotiff.cpp 17 Feb 2006 15:46:06 -0000      1.157
+++ geotiff.cpp 21 Mar 2006 20:40:41 -0000
@@ -3786,6 +3786,20 @@
                 adfPixelScale[1] = fabs(adfGeoTransform[5]);
                 adfPixelScale[2] = 0.0;
 
+                /* if the metadata indicates the image is using "Pixel is
+                   point", scale the pixel size since it has been 
calculated
+                   incorrectly.  This is not the proper way to fix it, 
but it
+                   takes more knowledge to find the correct place to fix 
it */
+                if( poSrcDS->GetMetadataItem( GDALMD_AREA_OR_POINT )
+                    && 
EQUAL(poSrcDS->GetMetadataItem(GDALMD_AREA_OR_POINT),
+                                GDALMD_AOP_POINT) )
+                {
+                    int nRasterXSize = poSrcDS->GetRasterXSize();
+                    int nRasterYSize = poSrcDS->GetRasterYSize();
+                    adfPixelScale[0] *= 
(double)nRasterXSize/(nRasterXSize - 1);
+                    adfPixelScale[1] *= 
(double)nRasterYSize/(nRasterYSize - 1);
+                }
+
                 TIFFSetField( hTIFF, TIFFTAG_GEOPIXELSCALE, 3, 
adfPixelScale );
 
                 adfTiePoints[0] = 0.0;
--------------------------------------------------------------------------
Tim Beckmann
Software Project Lead
Science Applications International Corporation (SAIC)
Contractor to the USGS/EROS
Sioux Falls, SD 57198
605-594-2521    Phone
605-594-6940    Fax
tbeckman at usgs.gov