open_geotiff overview_level coord offset wrong for PixelIsPoint COGs

[brendancol]

Summary

open_geotiff(path, overview_level=N) returns coords that are silently
off by half a level-0 pixel (multiplied by the overview reduction
factor) when the source COG declares raster_type = PixelIsPoint
(GeoKey 1025 = 2).

PR #1641 inherits the level-0 georef on overview reads, scaling the
pixel size by width_full / width_overview. It keeps the level-0
origin_x / origin_y unchanged. That is correct for the default
PixelIsArea convention, where the origin is the upper-left corner of
pixel (0, 0) and the overview's upper-left corner sits in the same
place as level 0's. It is wrong for PixelIsPoint, where the origin
is the center of pixel (0, 0): an overview pixel that spans the
first scale_x columns of level 0 has its center at
origin_x + (scale_x - 1) * 0.5 * pixel_width_lvl0, not at
origin_x.

Reproduction

import numpy as np
import xarray as xr
from xrspatial.geotiff import to_geotiff, open_geotiff

H = W = 1024
data = np.arange(H * W, dtype=np.float32).reshape(H, W)
# PixelIsPoint: coords are pixel CENTERS. Pixel (0, 0) center is at (0, 0).
x = np.arange(W) * 10.0
y = -(np.arange(H) * 10.0)
da = xr.DataArray(data, coords={'y': y, 'x': x}, dims=('y', 'x'))
da.attrs['raster_type'] = 'point'
da.attrs['crs'] = 'EPSG:32610'

to_geotiff(da, 'pp.tif', cog=True, overview_levels=[1, 2])

da0 = open_geotiff('pp.tif', overview_level=0)
da1 = open_geotiff('pp.tif', overview_level=1)
print('level0 x[:3]:', da0.x.values[:3])  # [ 0. 10. 20.]
print('level1 x[:3]:', da1.x.values[:3])  # [ 0. 20. 40.]  WRONG
# Level-1 pixel 0 covers original pixels 0 and 1 (centers 0 and 10),
# so its center sits at x = 5, not x = 0.

Expected level1 x[:3] for PixelIsPoint: [5. 25. 45.].

Impact

Cat 5 (backend / convention inconsistency) plus Cat 3 (off-by-one /
boundary handling).

• da.sel(x=v, method='nearest') on the overview snaps to the wrong
  pixel for any v near a boundary because the coord vector is
  shifted by half an overview pixel.
• da.interp(...) produces silently-wrong interpolated values.
• Reproject / hillshade / slope-and-aspect chains that consume an
  overview-read DataArray pick up the wrong pixel positions and
  produce silently-wrong results downstream.

The bug only affects raster_type = PixelIsPoint COGs. That convention
is mainstream for DEMs (USGS, OpenTopography, Copernicus DEM all emit
it). PixelIsArea is unaffected.

Proposed fix

In extract_geo_info_with_overview_inheritance (xrspatial/geotiff/_geotags.py),
after computing scale_x and scale_y, branch on the inherited
raster_type:

• PixelIsArea: keep origin_x / origin_y (existing behaviour).
• PixelIsPoint: shift origin by (scale - 1) * 0.5 * pixel_size_lvl0
  along each axis so the overview pixel-0 center matches the centroid
  of the level-0 pixels it covers.

Add regression coverage in
xrspatial/geotiff/tests/test_overview_geo_inheritance_1640.py
(or a follow-on file) for raster_type='point' across all four
backends and at least two reduction levels.

Categories: 3, 5 (accuracy sweep).