shift_origin not necessary when converting to pixel as point

[richardhchen]

Context

When comparing the DEM that comes with RTC products (*_dem.tif) against the original NED, I observed an obvious shift between them. After going through get_dem.py, I realized that NED was shifted by half a pixel towards the lower right corner in order to convert to PixelIsPoint.

However, I believe that this manual shift/offset is not necessary because GDAL will handle this offset when you change AREA_OR_POINT. Please see below for the reasons and an example.

AREA_OR_POINT

• Used in GDAL Raster Data Model to indicate whether a pixel value should be assumed to represent a sampling over the region of the pixel or a point sample at the center of the pixel.

• It does not influence the interpretation of georeferencing, which remains AREA oriented. The coefficients in GeoTransform should always be considered as area-based, i.e., (GT(0), GT(3)) position is ALWAYS the top left corner of the top left pixel of the raster.

• In ancient GDAL versions, PixelIsPoint/PixelIsArea of a GeoTIFF file was misinterpreted w.r.t. GeoTransform, causing an offset by half a pixel. This has been fixed in the recent GDAL versions. See RFC 33: GTiff - Fixing PixelIsPoint Interpretation. Users can set GTIFF_POINT_GEO_IGNORE to TRUE to revert back to the behavior of ancient GDAL versions.

• By default (GTIFF_POINT_GEO_IGNORE = FALSE), changing AREA_OR_POINT for a GeoTIFF file will NOT affect the GeoTransform coefficients when read in by GDAL, nor the corner coordinates seen in gdalinfo, as they are both area-based. What is affected is the model tiepoint in ModelTiepointTag of GeoTIFF file, which can seen as Tag 33922 in tiffinfo.

• Therefore, there is no need to manually do an offset of half a pixel when changing AREA_OR_POINT. GDAL will by default adjust the model tiepoint in ModelTiepointTag for a GeoTIFF file.

Example

import shutil
import tifffile
from osgeo import gdal

tif_file_1 = 'n63w144.tiff'
ds1 = gdal.Open(tif_file_1)
print(f"GeoTransform = {ds1.GetGeoTransform()}")
print(f"AREA_OR_POINT = {ds1.GetMetadataItem('AREA_OR_POINT')}")
ds1 = None

GeoTransform = (-144.000556, 9.2592592593e-05, 0.0, 63.000556, 0.0, -9.2592592593e-05)
AREA_OR_POINT = Area

tif1 = tifffile.TiffFile(tif_file_1)
print(f"{tif1.pages[0].tags[33922].name} = {tif1.pages[0].tags[33922].value}")
print(f"GTRasterTypeGeoKey (1: PixelIsArea; 2: PixelIsPoint) = {tif1.pages[0].tags[34735].value[11]}")
tif1 = None

ModelTiepointTag = (0.0, 0.0, 0.0, -144.000556, 63.000556, 0.0)
GTRasterTypeGeoKey (1: PixelIsArea; 2: PixelIsPoint) = 1

We can see that the coordinates are consistent in both GeoTransform and ModelTiePoint, as this raster is PixelIsArea.

Now, we change AREA_OR_POINT to POINT (i.e., PixelIsPoint).

tif_file_2 = 'n63w144_pixel_is_point.tiff'
shutil.copyfile(tif_file_1, tif_file_2)
ds2 = gdal.Open(tif_file_2, gdal.GA_Update)
ds2.SetMetadataItem('AREA_OR_POINT', 'Point')
print(f"GeoTransform = {ds2.GetGeoTransform()}")
print(f"AREA_OR_POINT = {ds2.GetMetadataItem('AREA_OR_POINT')}")
ds2 = None

GeoTransform = (-144.000556, 9.2592592593e-05, 0.0, 63.000556, 0.0, -9.2592592593e-05)
AREA_OR_POINT = Point

tif2 = tifffile.TiffFile(tif_file_2)
print(f"{tif2.pages[0].tags[33922].name} = {tif2.pages[0].tags[33922].value}")
print(f"GTRasterTypeGeoKey (1: PixelIsArea; 2: PixelIsPoint) = {tif2.pages[0].tags[34735].value[11]}")
tif2 = None

ModelTiepointTag = (0.0, 0.0, 0.0, -144.0005097037037, 63.000509703703706, 0.0)
GTRasterTypeGeoKey (1: PixelIsArea; 2: PixelIsPoint) = 2

The coefficients in GeoTransform remain unchanged because they are always area-based in GDAL. However, ModelTiePoint is offset by half a pixel towards the lower right corner, which proves that changing AREA_OR_POINT will make the offset for the tiepoint of GeoTIFF file. There is no need to manually do an extra offset.

Now, we change AREA_OR_POINT to POINT but with GTIFF_POINT_GEO_IGNORE set to TRUE.

tif_file_3 = 'n63w144_pixel_is_point_old_gdal.tiff'
shutil.copyfile(tif_file_1, tif_file_3)
gdal.SetConfigOption('GTIFF_POINT_GEO_IGNORE', 'True')
ds3 = gdal.Open(tif_file_3, gdal.GA_Update)
ds3.SetMetadataItem('AREA_OR_POINT', 'Point')
ds3 = None

Before checking GeoTransform again, set GTIFF_POINT_GEO_IGNORE to FALSE to see the effect under default GDAL configuration.

gdal.SetConfigOption('GTIFF_POINT_GEO_IGNORE', 'False')
ds3 = gdal.Open(tif_file_3)
print(f"GeoTransform = {ds3.GetGeoTransform()}")
print(f"AREA_OR_POINT = {ds3.GetMetadataItem('AREA_OR_POINT')}")
ds3 = None

GeoTransform = (-144.00060229629628, 9.2592592593e-05, 0.0, 63.0006022962963, 0.0, -9.2592592593e-05)
AREA_OR_POINT = Point

tif3 = tifffile.TiffFile(tif_file_3)
print(f"{tif3.pages[0].tags[33922].name} = {tif3.pages[0].tags[33922].value}")
print(f"GTRasterTypeGeoKey (1: PixelIsArea; 2: PixelIsPoint) = {tif3.pages[0].tags[34735].value[11]}")
tif3 = None

ModelTiepointTag = (0.0, 0.0, 0.0, -144.000556, 63.000556, 0.0)
GTRasterTypeGeoKey (1: PixelIsArea; 2: PixelIsPoint) = 2

As seen in GeoTransform, this will cause an offset by half a pixel towards the upper left corner. This behavior is not desired as it changes the spatial extent of the raster.

Summary

In get_dem.py or set_pixel_as_point, the manual offset of half a pixel is not necessary (unless GTIFF_POINT_GEO_IGNORE was set to TRUE in RTC jobs). The extra offset will cause the spatial extent of a raster to be actually shifted by half a pixel, which should not be the case when we just want to express pixels as area or point.

References

1. OGC GeoTIFF Standard

2. GDAL - Raster Data Model

3. GDAL - GeoTIFF File Format