Three bugs in open_geotiff and to_geotiff cause silent data loss during read-write round-trips.
Bug 1: Windowed read ignores PixelIsPoint raster type
open_geotiff with window= always applies a half-pixel coordinate offset (__init__.py:216-217), even for PixelIsPoint files where no offset should exist.
The non-windowed path gets this right. _geo_to_coords() checks geo_info.raster_type and skips the offset for PixelIsPoint. But the windowed override at lines 212-218 unconditionally adds pixel_width * 0.5 and pixel_height * 0.5.
For a 1-arc-second DEM (~30m pixels), this shifts coordinates by ~15m. Enough to misalign stacked layers.
Bug 2: CRS dropped for non-EPSG coordinate systems
to_geotiff converts CRS to an EPSG integer before writing. If the CRS has no EPSG code (custom local grids, compound CRS, some older national systems), _wkt_to_epsg() returns None and write() receives crs_epsg=None. No CRS gets written to the file.
This fails silently, no warning, no error. A round-trip on a file with a custom CRS produces a file with no spatial reference at all.
The root cause: write() and build_geo_tags() only accept EPSG integers. There is no code path for writing raw WKT or PROJ strings to the GeoTIFF's GeoAsciiParamsTag.
Bug 3: NaN not restored to nodata sentinel before writing
open_geotiff replaces nodata sentinel values (e.g. -9999) with NaN on read (lines 289-303). to_geotiff passes the original nodata value to the GDAL_NODATA tag (line 446) but does not convert NaN pixels back to the sentinel value before writing pixel data.
The written file ends up with:
- GDAL_NODATA tag:
-9999.0
- Actual pixel values:
NaN where -9999 should be
Tools that check for the sentinel value rather than NaN will miss the nodata pixels. This also breaks integer-only workflows where the original file had integer nodata.
Expected behavior
- Windowed reads should respect PixelIsPoint raster type (skip half-pixel offset)
- CRS should survive round-trips even without an EPSG code (write WKT to GeoAsciiParamsTag)
- NaN should be converted back to the nodata sentinel before writing pixel data
Three bugs in
open_geotiffandto_geotiffcause silent data loss during read-write round-trips.Bug 1: Windowed read ignores PixelIsPoint raster type
open_geotiffwithwindow=always applies a half-pixel coordinate offset (__init__.py:216-217), even for PixelIsPoint files where no offset should exist.The non-windowed path gets this right.
_geo_to_coords()checksgeo_info.raster_typeand skips the offset for PixelIsPoint. But the windowed override at lines 212-218 unconditionally addspixel_width * 0.5andpixel_height * 0.5.For a 1-arc-second DEM (~30m pixels), this shifts coordinates by ~15m. Enough to misalign stacked layers.
Bug 2: CRS dropped for non-EPSG coordinate systems
to_geotiffconverts CRS to an EPSG integer before writing. If the CRS has no EPSG code (custom local grids, compound CRS, some older national systems),_wkt_to_epsg()returnsNoneandwrite()receivescrs_epsg=None. No CRS gets written to the file.This fails silently, no warning, no error. A round-trip on a file with a custom CRS produces a file with no spatial reference at all.
The root cause:
write()andbuild_geo_tags()only accept EPSG integers. There is no code path for writing raw WKT or PROJ strings to the GeoTIFF's GeoAsciiParamsTag.Bug 3: NaN not restored to nodata sentinel before writing
open_geotiffreplaces nodata sentinel values (e.g. -9999) with NaN on read (lines 289-303).to_geotiffpasses the original nodata value to the GDAL_NODATA tag (line 446) but does not convert NaN pixels back to the sentinel value before writing pixel data.The written file ends up with:
-9999.0NaNwhere -9999 should beTools that check for the sentinel value rather than NaN will miss the nodata pixels. This also breaks integer-only workflows where the original file had integer nodata.
Expected behavior