shift_origin not necessary when converting to pixel as point

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.

ASFHyP3 / hyp3-lib