projected raster support, fetch retry

CHANGELOG.md

@@ -5,6 +5,15 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
+## [0.4.4] - 2026-06-25
+
+### Added
+- support for projected input rasters in the raster command
+- add 'save_shift' to transform_raster api/cli
+
+### CHANGED
+- rety failed downloads, such as FES
+
 ## [0.4.3] - 2026-05-04
 
 ### Added

CITATION.cff

@@ -14,6 +14,6 @@ authors:
 
 website: https://transformez.readthedocs.io
 title: "Transformez"
-version: 0.4.3
-date-released: 2026-06-04
+version: 0.5.0
+date-released: 2026-06-06
 url: "https://github.com/continuous-dems/transformez"

README.md

@@ -7,7 +7,7 @@
 <p align="center"><strong>Global vertical datum transformations, simplified.</strong></p>
 
 <p align="center">
-  <a href="https://github.com/continuous-dems/transformez"><img src="https://img.shields.io/badge/version-0.4.3-blue.svg" alt="Version"></a>
+  <a href="https://github.com/continuous-dems/transformez"><img src="https://img.shields.io/badge/version-0.5.0-blue.svg" alt="Version"></a>
   <a href="LICENSE"><img src="https://img.shields.io/badge/license-MIT-green.svg" alt="License"></a>
   <a href="https://www.python.org/"><img src="https://img.shields.io/badge/python-3.12+-yellow.svg" alt="Python"></a>
   <a href="https://badge.fury.io/py/transformez"><img src="https://badge.fury.io/py/transformez.svg" alt="PyPI version"></a>

src/transformez/api.py

@@ -113,6 +113,8 @@ def generate_grid(
     increment: Union[str, float],
     datum_in: str,
     datum_out: str,
+    epoch_in: str = "2010.0",
+    epoch_out: str = "2010.0",
     decay_pixels: int = 100,
     out_fn: Optional[str] = None,
     cache_dir: Optional[str] = None,
@@ -126,6 +128,8 @@ def generate_grid(
         increment: Resolution (e.g., '3s' or 0.0008333).
         datum_in: Source datum (e.g., 'mllw', '5703').
         datum_out: Target datum (e.g., '4979', '6319').
+        epoch_in: Source epoch (e.g., '2010.0')
+        epoch_out: Target epoch (e.g., '2010.0')
         decay_pixels: Set the pixel decay in case extrapolation is required.
         out_fn: If provided, saves the grid to this file (.tif or .gtx).
         cache_dir: Path to store downloaded grids.
@@ -167,6 +171,8 @@ def generate_grid(
         epsg_out=epsg_out,
         geoid_in=geoid_in,
         geoid_out=geoid_out,
+        epoch_in=epoch_in,
+        epoch_out=epoch_out,
         decay_pixels=decay_pixels,
         cache_dir=cache_dir,
         use_stations=use_stations,
@@ -196,6 +202,7 @@ def transform_raster(
     z_unit_in: Optional[str] = "auto",
     z_unit_out: Optional[str] = "auto",
     use_stations: bool = False,
+    save_shift: bool = False,
     verbose: bool = False,
 ) -> Optional[str]:
     """Apply a vertical datum transformation directly to an existing raster file.
@@ -209,23 +216,57 @@ def transform_raster(
         cache_dir: Path to store downloaded grids.
         z_unit_in: Input DEM z units.
         z_unit_out: Output DEM z units.
+        use_stations: Generate the shift grid from available tide stations,
+        safe_shift: Save the generated shift raster to disk.
         verbose: Enable debug logging.
 
     Returns:
         str: The path to the transformed output raster, or None if failed.
     """
 
     import rasterio
+    from rasterio.warp import transform_bounds, reproject, Resampling
+    from rasterio.transform import from_bounds
+    from fetchez.spatial import Region
 
     if not os.path.exists(input_raster):
         logger.error(f"Input raster not found: {input_raster}")
         return None
 
     with rasterio.open(input_raster) as src:
-        bounds = src.bounds
-        region_obj = Region(bounds.left, bounds.right, bounds.bottom, bounds.top)
+        native_crs = src.crs
+        native_bounds = src.bounds
+        native_transform = src.transform
         nx, ny = src.width, src.height
 
+    is_projected = native_crs.is_projected if native_crs else False
+    if is_projected:
+        logger.info(
+            f"Projected CRS detected ({native_crs}). Extracting WGS84 envelope..."
+        )
+        w, s, e, n = transform_bounds(native_crs, "EPSG:4326", *native_bounds)
+
+        buffer = 0.05
+        region_obj = Region(w - buffer, e + buffer, s - buffer, n + buffer)
+        logger.info(f"Using WGS84 region: {region_obj}")
+
+        inc_deg = 3.0 / 3600.0
+        vt_nx = int((region_obj.xmax - region_obj.xmin) / inc_deg)
+        vt_ny = int((region_obj.ymax - region_obj.ymin) / inc_deg)
+    else:
+        region_obj = Region(
+            native_bounds.left,
+            native_bounds.right,
+            native_bounds.bottom,
+            native_bounds.top,
+        )
+        vt_nx, vt_ny = nx, ny
+
+    # with rasterio.open(input_raster) as src:
+    #     bounds = src.bounds
+    #     region_obj = Region(bounds.left, bounds.right, bounds.bottom, bounds.top)
+    #     nx, ny = src.width, src.height
+
     epsg_in, geoid_in = _parse_datum(datum_in)
     epsg_out, geoid_out = _parse_datum(datum_out)
 
@@ -266,6 +307,47 @@ def transform_raster(
         logger.error("Failed to generate shift array for the raster bounds.")
         return None
 
+    if is_projected:
+        logger.info("Warping shift grid back to native raster projection...")
+        wgs_transform = from_bounds(
+            region_obj.xmin,
+            region_obj.ymin,
+            region_obj.xmax,
+            region_obj.ymax,
+            vt_nx,
+            vt_ny,
+        )
+        native_shift_array = np.zeros((ny, nx), dtype=np.float32)
+
+        reproject(
+            source=shift_array,
+            destination=native_shift_array,
+            src_transform=wgs_transform,
+            src_crs="EPSG:4326",
+            dst_transform=native_transform,
+            dst_crs=native_crs,
+            resampling=Resampling.bilinear,
+        )
+
+        shift_array = native_shift_array
+
+    if save_shift:
+        shift_fn = f"{os.path.splitext(output_raster)[0]}_shiftgrid.tif"
+        logger.info(f"Saving aligned shift grid to {shift_fn}...")
+        with rasterio.open(
+            shift_fn,
+            "w",
+            driver="GTiff",
+            height=ny,
+            width=nx,
+            count=1,
+            dtype=shift_array.dtype,
+            crs=native_crs,
+            transform=native_transform,
+            nodata=-9999.0,
+        ) as dst:
+            dst.write(shift_array, 1)
+
     success = GridEngine.apply_vertical_shift(
         input_raster,
         shift_array,

src/transformez/cli.py

@@ -302,6 +302,11 @@ def transform_grid(
     is_flag=True,
     help="Force RBF interpolation using live tide stations instead of global satellite models.",
 )
+@click.option(
+    "--save-shift",
+    is_flag=True,
+    help="Save the aligned vertical shift grid to disk alongside the output DEM.",
+)
 def transform_raster(
     input_file,
     input_datum,
@@ -311,6 +316,7 @@ def transform_raster(
     out,
     decay_pixels,
     use_stations,
+    save_shift,
 ):
     """Apply a vertical datum shift to an existing DEM."""
 
@@ -326,6 +332,7 @@ def transform_raster(
         z_unit_in=in_units,
         z_unit_out=out_units,
         use_stations=use_stations,
+        save_shift=save_shift,
         verbose=True,
     )
 

src/transformez/grid_engine.py

@@ -27,6 +27,12 @@
 logger = logging.getLogger(__name__)
 
 
+class GridCorruptionError(Exception):
+    """Raised when a fetched grid is corrupted and needs re-downloading."""
+
+    pass
+
+
 def plot_grid(grid_array, region, title="Vertical Shift Preview"):
     """Plot the transformation grid using Matplotlib."""
 
@@ -168,6 +174,34 @@ def load_and_interpolate(source_files, target_region, nx, ny, decay_pixels=100):
                     valid_mask = ~np.isnan(temp_buffer)
                     mosaic[valid_mask] = temp_buffer[valid_mask]
 
+            except Exception as e:
+                error_msg = str(e)
+
+                if any(
+                    err in error_msg
+                    for err in [
+                        "-101",
+                        "HDF error",
+                        "not recognized as a supported file format",
+                        "RasterioIOError",
+                    ]
+                ):
+                    logger.error(f" CRITICAL: Corrupted grid chunk detected in {fn}!")
+
+                    real_path = fn.split(":")[1] if fn.startswith("netcdf:") else fn
+                    if os.path.exists(real_path):
+                        logger.warning(
+                            f"Auto-deleting corrupted cache file to force re-fetch: {real_path}"
+                        )
+                        try:
+                            os.remove(real_path)
+                        except OSError:
+                            pass
+
+                    raise GridCorruptionError(f"Corrupted file deleted: {real_path}")
+
+                logger.exception(f"Failed to reproject {fn}: {e}")
+                continue
             # except Exception as e:
             #     error_msg = str(e)
 
@@ -189,9 +223,9 @@ def load_and_interpolate(source_files, target_region, nx, ny, decay_pixels=100):
             #     # For all other normal errors, log and continue as usual
             #     logger.exception(f"Failed to reproject {fn}: {e}")
             #     continue
-            except Exception as e:
-                logger.exception(f"Failed to reproject {fn}: {e}")
-                continue
+            # except Exception as e:
+            #     logger.exception(f"Failed to reproject {fn}: {e}")
+            #     continue
 
         # Fill inland areas (decaying to 0) before we clear the remaining NaNs
         # mosaic = GridEngine.fill_nans(mosaic, decay_pixels=decay_pixels)