Supporting TIFF AREA_OR_POINT

[gaopeng]

The two attached images show the problem. The only apparent difference between these images is the AREA_OR_POINT tag. This tag was changed in a custom application to prove a point. They should be 0.5 pixel apart.

I think one way to fix the problem is to apply the 0.5 pixel shift on image coordinates when they are used to compute geotransform.

[gaopeng]

I placed the sample images, GDAL_Ticket_3838.zip, at  ftp://GDAL@ftp.esri.com.

[rsbivand]

Replying to gaopeng:

I placed the sample images, GDAL_Ticket_3838.zip, at  ftp://GDAL@ftp.esri.com.

Please provide them on an open site - this one is password protected. The errant behaviour of ArcMap? in using this support metadata tag for a different raster model is causing havoc. IMHO Arc must rather drop using this tag and simply set up the GeoTransform? correctly.

WRT http://trac.osgeo.org/gdal/ticket/3837, I think that trying to distinguish between ESRI-originated GTiff and regular GDAL ones using the metadata tag as specified isn't possible, so that arbitrarily assuming that all with AREA_OR_POINT=Point are wrong ESRI files isn't tenable.

Can goapeng comment on ESRI usage? Can ESRI emit correctly constructed GeoTransform? values and avoid using this metadata item as a flag to switch between two kinds of georeferencing?

Can we keep all of this on one ticket, please?

Apparently this affects everything subsequent to the named version, so it needs resolution, but probably not by changing GeoTransform? conditional on AREA_OR_POINT unless this is also conditioned on the user requesting the change (that is, this interpretation of AREA_OR_POINT rather than the specified GDAL interpretation). Best would be for ESRI to back out of an apparent misuse of the metadata item to disguise two raster representations, where GDAL seems only to have one.

[gaopeng]

I didn't realize there was a Ticket already about the same issue. I'd like to use this Ticket to get it resolved if it's OK.

ArcGIS simply reads the geotransform or GCPs from GDAL and use it. The issue here is if GDAL will use the AREA_OR_POINT tag in computing its geotransform. I think it should use the tag to adjust the geotransform in cases where it's computed from tie points, and tie points are defined in image column and row numbers. For example, in the following, if TIFFTAG_GEOTIEPOINTS are defined in pixel space, the pixel location (0, 0) will mean different things depending on PixelIsPoint? (center of the pixel) or PixelIsArea? (UL corner of the pixel).

        if( TIFFGetField(hTIFF,TIFFTAG_GEOPIXELSCALE,&nCount,&padfScale )
            && nCount >= 2 
            && padfScale[0] != 0.0 && padfScale[1] != 0.0 )
        {
            adfGeoTransform[1] = padfScale[0];
            adfGeoTransform[5] = - ABS(padfScale[1]);

            if( TIFFGetField(hTIFF,TIFFTAG_GEOTIEPOINTS,&nCount,&padfTiePoints )
                && nCount >= 6 )
            {
                adfGeoTransform[0] =
                    padfTiePoints[3] - padfTiePoints[0] * adfGeoTransform[1];
                adfGeoTransform[3] =
                    padfTiePoints[4] - padfTiePoints[1] * adfGeoTransform[5];

                bGeoTransformValid = TRUE;
            }
        }

[warmerdam]

To my immense consternation section 2.5.2.2 of the geotiff specification does make it clear that the georeferencing semantics are altered by half a pixel in the pixel-as-area case refuting 15 years of stubborn behavior on my part.

I shall attempt to prepare a GDAL RFC to address how this transition will be handled.

Details to follow.

[warmerdam]

I will note there is nothing special about the files in GDAL_Ticket_3838.zip. I don't see much need to make them more public.

[rhijmans]

I think there is a little more to it. If we have:

x - left x

y - lower y

dx - x resolution

dy - y resolution

ncols - number of columns

nrows - number of rows

You could argue that the correction should be:

leftx = x - 0.5 * dx

rightx = x + ncols * dx + 0.5 * dx

lowery = y - 0.5 * dy

uppery = y + nrows * dy + 0.5 * dy

In that case the resolution changes as:

(rightx - leftx) / ncols != dx

The question is whether the resolution in a "AREA_OR_POINT=Point" file was computed from the distance between two points (cells), or (for x) by doing (rightx - leftx) / ncols. The latter is wrong as it should be (rightx - leftx) / (ncols - 1) because you have effectively one column less.

Here is what ArcGIS does, I think.

leftx = x - 0.5 * dx

rightx = x + ncols * dx - 0.5 * dx

lowery = y + 0.5 * dy

uppery = y + nrows * dy + 0.5 * dy

In effect it shifts the entire raster to the left and up. In this case (rightx - leftx) / ncols == dx which looks right, but it is in fact inconsistent because the resolution should change if you compute it this way (a column has been added). And it is strange because you might as well do the below

leftx = x + 0.5 * dx

rightx = x + ncols * dx + 0.5 * dx

lowery = y - 0.5 * dy

uppery = y + nrows * dy - 0.5 * dy

(i.e. shift the whole raster to the right and down) or some other combination.

So this really cannot be right. But can we know for sure what was intended? Is it documented how this should be done?

In conclusion, the change you suggest is incomplete, I think, as the resolution (adfGeoTransform[1] and adfGeoTransform[5]) also need to be changed.

[rsbivand]

Replying to warmerdam:

Frank - do you intend AREA_OR_PIXEL in the RFC to be that, or is it a typo for AREA_OR_POINT (with two possible values, AREA or POINT, not PIXEL?

[rhijmans]

I want to retract the above. The most reasonable assumption is that the resolution reported is correct and that the adjustment should therefore be as gaopeng suggested.

My confusion came from the assumption that the extent (xmin, ymin, xmax, ymax) is known and the resolution is derived; but if geotiff stores xmin/ymin and resolution, we must assume it is all correct. But I do not know for sure what the file stores; all I know is what gdalinfo reports. It would be nice if gdalinfo would also report ymax and xmax _if_ the file stores these, as they could then be available at a higher precision than after computing them if the resolution is a number with a fraction; but that is another matter.

[warmerdam]

Roger, I have updated the RFC - just a bunch of typos.

http://trac.osgeo.org/gdal/wiki/rfc33_gtiff_pixelispoint

[warmerdam]

Fix addressed in:

http://trac.osgeo.org/gdal/wiki/rfc33_gtiff_pixelispoint

Applied in trunk (r21158), 1.7 (r21159), 1.6 (r21160) and 1.6-esri (r21161).

[rouault]

Milestone 1.6.4 deleted