From a988820b448032bfc3355b0fd0701968c3081c51 Mon Sep 17 00:00:00 2001
From: konstntokas <konstantin.ntokas@brockmann-consult.de>
Date: Wed, 8 Jul 2026 10:15:13 +0200
Subject: [PATCH 1/9] make sen1 grd lazy

---
 xarray_eopf/amodes/sentinel1.py | 670 +++++++++++++++++++-------------
 1 file changed, 400 insertions(+), 270 deletions(-)

diff --git a/xarray_eopf/amodes/sentinel1.py b/xarray_eopf/amodes/sentinel1.py
index cb1b29f..7b5bd53 100644
--- a/xarray_eopf/amodes/sentinel1.py
+++ b/xarray_eopf/amodes/sentinel1.py
@@ -9,10 +9,12 @@
 import warnings
 from abc import ABC
 from collections.abc import Iterable, Sequence
-from dataclasses import dataclass
+from dataclasses import dataclass, fields
 from typing import Any, Callable, Literal
+import uuid
 
 import dask.array as da
+import fsspec
 import flox.xarray
 import numpy as np
 import pyproj
@@ -23,7 +25,13 @@
 from xcube_resampling.constants import SpatialAggMethods, SpatialInterpMethods
 from xcube_resampling.gridmapping import GridMapping
 from xcube_resampling.rectify import rectify_dataset
-from xcube_resampling.utils import reproject_bbox, transform_resolution
+from xcube_resampling.utils import (
+    reproject_bbox,
+    resolution_degrees_to_meters,
+    _reorganize_tiled_array,
+    transform_resolution,
+    SourceTileIndexing,
+)
 
 from xarray_eopf.amode import AnalysisMode, AnalysisModeRegistry
 from xarray_eopf.source import get_source_path
@@ -42,60 +50,23 @@
 class GridParams:
     """RTC grid parameters."""
 
-    slr0: float
-    d_slr: float
-    spacing_slr: float
+    gr0: float
+    d_gr: float
+    d_gr_scale: float
     az0: np.datetime64
     d_az: float
+    d_az_scale: float
     spacing_az: float
+    spacing_az_scale: float
 
     def __iter__(self):
-        return iter(("slr0", "d_slr", "spacing_slr", "az0", "d_az", "spacing_az"))
+        return (f.name for f in fields(self))
 
-    def __getitem__(self, key: str) -> Any:
+    def __getitem__(self, key: str):
         return getattr(self, key)
 
-    def __contains__(self, key: object) -> bool:
-        return key in {"slr0", "d_slr", "spacing_slr", "az0", "d_az", "spacing_az"}
-
-
-@dataclass
-class Acquisition:
-    """Simulated acquisition geometry."""
-
-    azimuth_time: xr.DataArray
-    distance: xr.DataArray
-    velocity: xr.DataArray
-    slant_range_time: xr.DataArray
-    gamma_area: xr.DataArray | None = None
-
-    def __iter__(self):
-        keys = ["azimuth_time", "distance", "velocity", "slant_range_time"]
-        if self.gamma_area is not None:
-            keys.append("gamma_area")
-        return iter(keys)
-
-    def __getitem__(self, key: str) -> xr.DataArray:
-        value = getattr(self, key)
-        if value is None:
-            raise KeyError(key)
-        return value
-
-    def __contains__(self, key: object) -> bool:
-        return key in set(iter(self))
-
-    def to_dataset(self) -> xr.Dataset:
-        dataset = xr.Dataset(
-            {
-                "azimuth_time": self.azimuth_time,
-                "distance": self.distance,
-                "velocity": self.velocity,
-                "slant_range_time": self.slant_range_time,
-            }
-        )
-        if self.gamma_area is not None:
-            dataset["gamma_area"] = self.gamma_area
-        return dataset
+    def __contains__(self, key):
+        return key in {f.name for f in fields(self)}
 
 
 class Sen1(AnalysisMode, ABC):
@@ -124,6 +95,8 @@ def process_metadata(self, datatree: xr.DataTree) -> dict:
 
 class Sen1GRD(Sen1):
     product_type = "GRDH"
+    cache_fs: fsspec.AbstractFileSystem | None = None
+    cache_uri: str | None = None
 
     def get_applicable_params(self, **kwargs) -> dict[str, Any]:
         params = {}
@@ -178,6 +151,11 @@ def get_applicable_params(self, **kwargs) -> dict[str, Any]:
             assert_arg_is_instance(apply_rtc, "apply_rtc", bool)
             params.update(apply_rtc=apply_rtc)
 
+        cache_uri = kwargs.get("cache_uri")
+        if cache_uri is not None:
+            assert_arg_is_instance(cache_uri, "cache_uri", str)
+            params.update(cache_uri=cache_uri)
+
         return params
 
     def convert_datatree(
@@ -192,7 +170,17 @@ def convert_datatree(
         footprint_scale_factor: tuple[float, float] = (3.0, 3.0),
         dem: xr.DataArray | None = None,
         apply_rtc: bool = True,
+        cache_uri: str | None = None,
     ) -> xr.Dataset:
+
+        if cache_uri is None:
+            self.cache_fs = fsspec.filesystem("file")
+            self.cache_uri = f"tmp/{uuid.uuid4().hex}"
+        else:
+            cache_uri = cache_uri.rstrip("/")
+            self.cache_fs, _ = fsspec.url_to_fs(cache_uri)
+            self.cache_uri = cache_uri
+
         # get dem data array
         if dem is None:
             if bbox is None:
@@ -243,22 +231,97 @@ def convert_datatree(
             )
 
         orbit = datatree[f"{group}/conditions/orbit"].to_dataset()
-        sat_position = orbit["position"]
+        sat_position = orbit["position"].compute()
 
         gcp = datatree[f"{group}/conditions/gcp"].to_dataset()
         time_slr_gcp = gcp["slant_range_time_gcp"]
 
         grid_params = self._get_grid_parameters(datatree, footprint_scale_factor)
 
-        return terrain_correct(
-            grd,
+        try:
+            return self._terrain_correct(
+                grd,
+                time_slr_gcp,
+                sat_position,
+                dem,
+                grid_params,
+                apply_rtc=apply_rtc,
+                interp_method=interp_methods,
+            )
+
+        except Exception as _:
+            self._cleanup()
+            raise
+
+    def _terrain_correct(
+        self,
+        data: xr.Dataset,
+        time_slr_gcp: xr.DataArray,
+        sat_position: xr.DataArray,
+        dem: xr.DataArray,
+        grid_params: GridParams | None = None,
+        apply_rtc: bool = True,
+        interp_method: Literal["nearest", "bilinear"] = "nearest",
+    ) -> xr.Dataset:
+        """Apply terrain correction to SAR data.
+
+        Args:
+            data: Input SAR dataset.
+            time_slr_gcp: GCP slant-range times.
+            sat_position: Satellite positions over time.
+            dem: DEM for terrain correction.
+            apply_rtc: Whether to apply radiometric terrain correction.
+            grid_params: Grid parameters for RTC.
+            interp_method: Interpolation method.
+
+        Returns:
+            Terrain-corrected dataset.
+
+        Raises:
+            ValueError: If RTC is enabled without grid parameters.
+        """
+        gm_dem = GridMapping.from_dataset(dem.to_dataset(name="dem"))
+        src_loc = get_source_location(
+            dem,
             time_slr_gcp,
             sat_position,
-            dem,
-            grid_params=grid_params,
-            apply_rtc=apply_rtc,
-            interp_method=interp_methods,
+            grid_params,
+            gm_dem,
+            apply_rtc,
         )
+        store = fsspec.get_mapper(f"{self.cache_uri}/src_location.zarr")
+        src_loc.to_zarr(store)
+
+        src_loc = xr.open_zarr(store)
+        geocoded = geocode_data(data, src_loc, grid_params, interp_method)
+
+        if apply_rtc:
+            if interp_method == "bilinear":
+                weights_fn = gamma_weights_bilinear
+            else:  # interp_method == "nearest"
+                weights_fn = gamma_weights_nearest
+            gamma_weights = apply_gamma_weights(src_loc, weights_fn, grid_params)
+            geocoded /= gamma_weights
+            rename_dict = {
+                name: name.replace("beta0", "gamma0") for name in geocoded.data_vars
+            }
+            geocoded = geocoded.rename(rename_dict)
+            for var in geocoded.data_vars:
+                geocoded[var].attrs.update(
+                    long_name="gamma nought backscatter coefficient",
+                    units="1",
+                )
+
+        geocoded = assign_grid_mapping(geocoded)
+        return geocoded
+
+    def _cleanup(self):
+        if not self.cache_uri:
+            return
+
+        fs, path = fsspec.url_to_fs(self.cache_uri)
+        if fs.exists(path):
+            fs.rm(path, recursive=True)
 
     @staticmethod
     def _get_grid_parameters(
@@ -275,21 +338,20 @@ def _get_grid_parameters(
             Grid parameters for terrain correction.
         """
 
-        group_VH = [x for x in dt.children if "VH" in x][0]
-        attrs = dt[f"{group_VH}"].attrs["other_metadata"]["image_annotation"][
+        group_vh = [x for x in dt.children if "VH" in x][0]
+        attrs = dt[f"{group_vh}"].attrs["other_metadata"]["image_annotation"][
             "image_information"
         ]
 
-        slant_range_spacing_m = attrs["range_pixel_spacing"] * footprint_scale_factor[1]
-        slant_range_time_interval_s = slant_range_spacing_m * 2 / _SPEED_OF_LIGHT
-
         return GridParams(
-            slr0=attrs["slant_range_time"],
-            d_slr=slant_range_time_interval_s,
-            spacing_slr=slant_range_spacing_m,
+            gr0=0.0,
+            d_gr=attrs["range_pixel_spacing"],
+            d_gr_scale=attrs["range_pixel_spacing"] * footprint_scale_factor[1],
             az0=np.datetime64(attrs["product_first_line_utc_time"]),
-            d_az=attrs["azimuth_time_interval"] * footprint_scale_factor[0],
-            spacing_az=attrs["azimuth_pixel_spacing"] * footprint_scale_factor[0],
+            d_az=attrs["azimuth_time_interval"],
+            spacing_az=attrs["azimuth_pixel_spacing"],
+            d_az_scale=attrs["azimuth_time_interval"] * footprint_scale_factor[0],
+            spacing_az_scale=attrs["azimuth_pixel_spacing"] * footprint_scale_factor[0],
         )
 
 
@@ -351,20 +413,20 @@ def convert_datatree(
         dataset.update(sub_dt.quality.to_dataset().drop_vars("calibration_constant"))
 
         # correct attributes and encoding
-        def _apply_valid_range(da, *, dtype=None, fill_value=None):
+        def _apply_valid_range(array, *, dtype=None, fill_value=None):
             if dtype is not None:
-                da = da.astype(dtype)
+                array = array.astype(dtype)
 
             if fill_value is not None:
-                da.encoding["_FillValue"] = fill_value
+                array.encoding["_FillValue"] = fill_value
 
-            eopf_attrs = da.attrs["_eopf_attrs"]
-            da.attrs.update(
+            eopf_attrs = array.attrs["_eopf_attrs"]
+            array.attrs.update(
                 valid_min=eopf_attrs["valid_min"],
                 valid_max=eopf_attrs["valid_max"],
             )
 
-            return da
+            return array
 
         dataset["inversion_quality"] = _apply_valid_range(
             dataset.inversion_quality,
@@ -473,8 +535,7 @@ def get_dem(
         bbox_wgs84 = bbox
 
     # get STAC items
-    STAC_URL = "https://stac.dataspace.copernicus.eu/v1"
-    client = pystac_client.Client.open(STAC_URL)
+    client = pystac_client.Client.open("https://stac.dataspace.copernicus.eu/v1")
     search = client.search(
         collections=["cop-dem-glo-30-dged-cog"],
         bbox=list(bbox_wgs84),
@@ -509,7 +570,7 @@ def get_dem(
     return dem
 
 
-def convert_dem_to_ecef(dem: xr.DataArray, gm_dem: GridMapping) -> xr.DataArray:
+def convert_dem_to_ecef(dem: xr.DataArray, gm_dem_params: dict) -> xr.DataArray:
     """Convert a DEM from its native CRS to ECEF coordinates.
 
     Args:
@@ -520,38 +581,21 @@ def convert_dem_to_ecef(dem: xr.DataArray, gm_dem: GridMapping) -> xr.DataArray:
         DEM expressed in ECEF axes.
     """
 
-    x_dim, y_dim = gm_dem.xy_var_names
-    transformer = pyproj.Transformer.from_crs(gm_dem.crs, _CRS_ECEF, always_xy=True)
-
-    def _transform(
-        block_xx: np.ndarray, block_yy: np.ndarray, block_dem: np.ndarray
-    ) -> np.ndarray:
-        x, y, z = transformer.transform(block_xx, block_yy, block_dem)
-        return np.stack([x, y, z], axis=0)
-
-    xx, yy = da.meshgrid(
-        da.from_array(dem[x_dim].values, chunks=dem.data.chunks[1][0]),
-        da.from_array(dem[y_dim].values, chunks=dem.data.chunks[0][0]),
-        indexing="xy",
-    )
+    x_dim, y_dim = gm_dem_params["xy_var_names"]
+    xx, yy = np.meshgrid(dem[x_dim].values, dem[y_dim].values, indexing="xy")
 
-    xyz_transformed = da.map_blocks(
-        _transform,
-        xx,
-        yy,
-        dem.data,
-        dtype=np.float32,
-        chunks=(3, *dem.data.chunks),
+    transformer = pyproj.Transformer.from_crs(
+        gm_dem_params["crs"], _CRS_ECEF, always_xy=True
     )
+    x, y, z = transformer.transform(xx, yy, dem.values)
 
     return xr.DataArray(
-        xyz_transformed,
+        np.stack([x, y, z], axis=0),
         dims=("axis", y_dim, x_dim),
         coords={
             y_dim: dem[y_dim].data,
             x_dim: dem[x_dim].data,
             "axis": ["x", "y", "z"],
-            "spatial_ref": xr.DataArray(0, attrs=gm_dem.crs.to_cf()),
         },
     )
 
@@ -753,15 +797,19 @@ def newton(
 
 
 def backward_geocode(
-    dem_ecef: xr.DataArray,
-    pos_coeff: xr.DataArray,
-    vel_coeff: xr.DataArray,
+    dem: xr.DataArray,
+    pos_coeff: xr.DataArray = None,
+    vel_coeff: xr.DataArray = None,
+    gr_coeff: xr.DataArray = None,
+    grid_params: GridParams = None,
+    apply_rtc: bool = True,
+    gm_dem_params: dict = None,
     method="newton",
     tol=1.0,
     speed=7500.0,
     maxiter=10,
     t_shift=-0.1,
-) -> tuple[xr.DataArray, xr.DataArray, xr.DataArray]:
+) -> xr.Dataset:
     """Compute orbit time and vectors for a DEM using inverse geocoding.
 
     Args:
@@ -775,11 +823,19 @@ def backward_geocode(
         t_shift: Time shift for the secant method.
 
     Returns:
-        Orbit time, distance vector, and velocity vector.
+        Orbit time and distance vector
 
     Raises:
         ValueError: If the method is not supported.
     """
+    assert pos_coeff is not None
+    assert vel_coeff is not None
+    assert gr_coeff is not None
+    assert grid_params is not None
+    assert gm_dem_params is not None
+
+    dem_ecef = convert_dem_to_ecef(dem, gm_dem_params)
+
     f = functools.partial(zero_doppler, dem_ecef, pos_coeff, vel_coeff)
 
     t0 = xr.zeros_like(dem_ecef.sel(axis="x", drop=True), dtype="float64")
@@ -797,48 +853,30 @@ def backward_geocode(
     else:
         raise ValueError("method needs to be either 'secant' or 'newton'")
 
-    dist, vel = payload
-    return time_orbit, dist, vel
-
-
-def simulate_acquisition(
-    dem_ecef: xr.DataArray,
-    gm_dem: GridMapping,
-    sat_position: xr.DataArray,
-    apply_rtc: bool = True,
-) -> Acquisition:
-    """Simulate SAR acquisition geometry for a DEM.
-
-    Args:
-        dem_ecef: DEM in ECEF coordinates.
-        gm_dem: GridMapping of the DEM data array.
-        sat_position: Satellite positions over time.
-        apply_rtc: Whether to compute gamma area.
+    dist, _ = payload
 
-    Returns:
-        Simulated acquisition geometry.
-    """
-    pos_coeff = fit_position(sat_position)
-    vel_coeff = poly_derivative(pos_coeff)
-    time_orbit, dist, vel = backward_geocode(dem_ecef, pos_coeff, vel_coeff)
+    # apply bistatic correction
+    slant_range = np.sqrt((dist**2).sum("axis"))
+    time_orbit += slant_range / _SPEED_OF_LIGHT
 
+    # recalculate slant range
+    sat = xr.polyval(time_orbit, pos_coeff)
+    dist = dem_ecef - sat
     slant_range = np.sqrt((dist**2).sum("axis"))
     time_slr = 2 * slant_range / _SPEED_OF_LIGHT
 
-    out = Acquisition(
-        azimuth_time=orbit_to_az(time_orbit, pos_coeff.attrs["epoch"]),
-        distance=dist,
-        velocity=vel.transpose(*dist.dims),
-        slant_range_time=time_slr,
-    )
-
+    # convert to ground range
+    azimuth_time = orbit_to_az(time_orbit, pos_coeff.attrs["epoch"])
+    ground_range = get_ground_range(gr_coeff, azimuth_time, time_slr)
+    out = xr.Dataset({"azimuth_time": azimuth_time, "ground_range": ground_range})
     if apply_rtc:
-        out.gamma_area = compute_gamma_area(dem_ecef, gm_dem, dist / slant_range)
-
+        out["gamma_area"] = compute_gamma_area(
+            dem_ecef, gm_dem_params, dist / slant_range
+        )
     return out
 
 
-def compute_dem_area(dem_ecef: xr.DataArray, gm_dem: GridMapping) -> xr.DataArray:
+def compute_dem_area(dem_ecef: xr.DataArray, gm_dem_params: dict) -> xr.DataArray:
     """Compute per-pixel surface area on the DEM in ECEF coordinates.
 
     Args:
@@ -848,7 +886,7 @@ def compute_dem_area(dem_ecef: xr.DataArray, gm_dem: GridMapping) -> xr.DataArra
     Returns:
         Area vectors per DEM pixel.
     """
-    x_dim, y_dim = gm_dem.xy_var_names
+    x_dim, y_dim = gm_dem_params["xy_var_names"]
     x = dem_ecef[x_dim]
     y = dem_ecef[y_dim]
 
@@ -870,17 +908,11 @@ def compute_dem_area(dem_ecef: xr.DataArray, gm_dem: GridMapping) -> xr.DataArra
     )
 
     # interpolate DEM to pixel corners
-    chunksizes = {key: val[0] for key, val in dem_ecef.chunksizes.items()}
     xyz_c = dem_ecef.interp(
-        {x_dim: x_corner},
+        {x_dim: x_corner, y_dim: y_corner},
         method="linear",
         kwargs={"fill_value": "extrapolate"},
-    ).chunk({x_dim: chunksizes[x_dim]})
-    xyz_c = xyz_c.interp(
-        {y_dim: y_corner},
-        method="linear",
-        kwargs={"fill_value": "extrapolate"},
-    ).chunk(chunksizes)
+    )
 
     # compute edge vectors
     dx = xyz_c.diff(x_dim)
@@ -893,10 +925,10 @@ def compute_dem_area(dem_ecef: xr.DataArray, gm_dem: GridMapping) -> xr.DataArra
     dy2 = dy.isel({x_dim: slice(None, -1)})
 
     # restore original coords
-    dx1 = dx1.assign_coords(dem_ecef.coords).chunk(chunksizes)
-    dy1 = dy1.assign_coords(dem_ecef.coords).chunk(chunksizes)
-    dx2 = dx2.assign_coords(dem_ecef.coords).chunk(chunksizes)
-    dy2 = dy2.assign_coords(dem_ecef.coords).chunk(chunksizes)
+    dx1 = dx1.assign_coords(dem_ecef.coords)
+    dy1 = dy1.assign_coords(dem_ecef.coords)
+    dx2 = dx2.assign_coords(dem_ecef.coords)
+    dy2 = dy2.assign_coords(dem_ecef.coords)
 
     # compute triangle areas
     cross1 = xr.cross(dx1, dy1, dim="axis") / 2
@@ -910,7 +942,9 @@ def compute_dem_area(dem_ecef: xr.DataArray, gm_dem: GridMapping) -> xr.DataArra
 
 
 def compute_gamma_area(
-    dem_ecef: xr.DataArray, gm_dem: GridMapping, direction: xr.DataArray
+    dem_ecef: xr.DataArray,
+    gm_dem_params: dict,
+    direction: xr.DataArray,
 ) -> xr.DataArray:
     """Compute gamma area by projecting DEM areas onto look direction.
 
@@ -922,7 +956,7 @@ def compute_gamma_area(
     Returns:
         Gamma area for each DEM pixel.
     """
-    area = compute_dem_area(dem_ecef, gm_dem)
+    area = compute_dem_area(dem_ecef, gm_dem_params)
     gamma_area = xr.dot(area, -direction, dim="axis")
     return gamma_area.where(gamma_area > 0, 0)
 
@@ -930,34 +964,34 @@ def compute_gamma_area(
 def sum_weights(
     weights: xr.DataArray,
     az_idx: xr.DataArray,
-    slr_idx: xr.DataArray,
+    gr_idx: xr.DataArray,
 ) -> xr.DataArray:
     """Accumulate weights into the SAR image grid.
 
     Args:
         weights: Weights to accumulate.
         az_idx: Azimuth indices.
-        slr_idx: Slant-range indices.
+        gr_idx: Ground-range indices.
 
     Returns:
         Accumulated weights on the SAR grid.
     """
     reduced = flox.xarray.xarray_reduce(
         weights,
-        slr_idx,
+        gr_idx,
         az_idx,
         func="sum",
         method="map-reduce",
     )
 
     return reduced.interp(
-        slr_idx=slr_idx,
+        gr_idx=gr_idx,
         az_idx=az_idx,
         method="nearest",
-    ).drop_vars(("az_idx", "slr_idx"))
+    ).drop_vars(("az_idx", "gr_idx"))
 
 
-def gamma_weights_bilinear(acq: xr.Dataset) -> xr.DataArray:
+def gamma_weights_bilinear(src_loc: xr.Dataset) -> xr.DataArray:
     """Compute bilinear gamma weights for the acquisition grid.
 
     Args:
@@ -966,29 +1000,29 @@ def gamma_weights_bilinear(acq: xr.Dataset) -> xr.DataArray:
     Returns:
         Gamma weights on the SAR grid.
     """
-    az_idx = acq.az_idx
-    slr_idx = acq.slr_idx
+    az_idx = src_loc.az_idx
+    gr_idx = src_loc.gr_idx
 
     az0 = np.floor(az_idx).astype(np.intp)
     az1 = np.ceil(az_idx).astype(np.intp)
-    slr0 = np.floor(slr_idx).astype(np.intp)
-    slr1 = np.ceil(slr_idx).astype(np.intp)
+    gr0 = np.floor(gr_idx).astype(np.intp)
+    gr1 = np.ceil(gr_idx).astype(np.intp)
 
-    w00 = abs((az1 - az_idx) * (slr1 - slr_idx))
-    w01 = abs((az1 - az_idx) * (slr0 - slr_idx))
-    w10 = abs((az0 - az_idx) * (slr1 - slr_idx))
-    w11 = abs((az0 - az_idx) * (slr0 - slr_idx))
+    w00 = abs((az1 - az_idx) * (gr1 - gr_idx))
+    w01 = abs((az1 - az_idx) * (gr0 - gr_idx))
+    w10 = abs((az0 - az_idx) * (gr1 - gr_idx))
+    w11 = abs((az0 - az_idx) * (gr0 - gr_idx))
 
-    gamma = acq.gamma_area
+    gamma = src_loc.gamma_area
     return (
-        sum_weights(gamma * w00, az0, slr0)
-        + sum_weights(gamma * w01, az0, slr1)
-        + sum_weights(gamma * w10, az1, slr0)
-        + sum_weights(gamma * w11, az1, slr1)
+        sum_weights(gamma * w00, az0, gr0)
+        + sum_weights(gamma * w01, az0, gr1)
+        + sum_weights(gamma * w10, az1, gr0)
+        + sum_weights(gamma * w11, az1, gr1)
     )
 
 
-def gamma_weights_nearest(acq: xr.Dataset) -> xr.DataArray:
+def gamma_weights_nearest(src_loc: xr.Dataset) -> xr.DataArray:
     """Compute nearest-neighbor gamma weights for the acquisition grid.
 
     Args:
@@ -998,13 +1032,13 @@ def gamma_weights_nearest(acq: xr.Dataset) -> xr.DataArray:
         Gamma weights on the SAR grid.
     """
 
-    az_idx = np.round(acq.az_idx).astype(np.intp)
-    slr_idx = np.round(acq.slr_idx).astype(np.intp)
-    return sum_weights(acq.gamma_area, az_idx, slr_idx)
+    az_idx = np.round(src_loc.az_idx).astype(np.intp)
+    gr_idx = np.round(src_loc.gr_idx).astype(np.intp)
+    return sum_weights(src_loc.gamma_area, az_idx, gr_idx)
 
 
 def apply_gamma_weights(
-    acq: Acquisition,
+    src_loc: xr.Dataset,
     func: Callable[..., xr.DataArray],
     params: GridParams,
 ) -> xr.DataArray:
@@ -1018,14 +1052,15 @@ def apply_gamma_weights(
     Returns:
         Gamma-corrected area per pixel.
     """
-    acq_ds = acq.to_dataset()
-    acq_ds["az_idx"] = (acq_ds.azimuth_time - params.az0) / _ONE_SECOND / params.d_az
-    acq_ds["slr_idx"] = (acq_ds.slant_range_time - params.slr0) / params.d_slr
+    src_loc["az_idx"] = (
+        (src_loc.azimuth_time - params.az0) / _ONE_SECOND / params.d_az_scale
+    )
+    src_loc["gr_idx"] = (src_loc.ground_range - params.gr0) / params.d_gr_scale
 
-    template = acq_ds.gamma_area * 0
+    template = src_loc.gamma_area * 0
+    area = xr.map_blocks(func, src_loc, template=template)
 
-    area = xr.map_blocks(func, acq_ds, template=template)
-    return area / (params.spacing_slr * params.spacing_az)
+    return area / (params.d_gr_scale * params.spacing_az_scale)
 
 
 def fit_ground_range(time_slr_gcp: xr.DataArray, deg: int = 8) -> xr.DataArray:
@@ -1045,21 +1080,28 @@ def fit_ground_range(time_slr_gcp: xr.DataArray, deg: int = 8) -> xr.DataArray:
 
     # polynomial fit per azimuth line
     coeff = []
-    for i, time in enumerate(x_gcp.azimuth_time.values):
-        coeff.append(np.polyfit(x_gcp[i, :], x_gcp.ground_range, deg=deg))
+    for i, time in enumerate(x_gcp["azimuth_time"].values):
+        coeff.append(np.polyfit(x_gcp[i, :], x_gcp["ground_range"], deg=deg))
     return xr.DataArray(
         coeff,
-        coords=dict(azimuth_time=x_gcp.azimuth_time, degree=np.arange(deg, -1, -1)),
+        coords=dict(azimuth_time=x_gcp["azimuth_time"], degree=np.arange(deg, -1, -1)),
         dims=("azimuth_time", "degree"),
         attrs=dict(mean=mean, std=std),
     )
 
 
+def get_ground_range(
+    coeff: xr.DataArray, time_az: xr.DataArray, time_slr: xr.DataArray
+) -> xr.DataArray:
+    coeff_interp = coeff.interp(azimuth_time=time_az).drop_vars("azimuth_time")
+    x_tgt = (time_slr - coeff.attrs["mean"]) / coeff.attrs["std"]
+    return (coeff_interp * x_tgt**coeff.degree).sum("degree")
+
+
 def geocode_data(
     data: xr.Dataset,
-    time_az: xr.DataArray,
-    time_slr: xr.DataArray,
-    time_slr_gcp: xr.DataArray,
+    src_loc: xr.Dataset,
+    grid_params: GridParams,
     interp_method: Literal["nearest", "bilinear"],
 ) -> xr.Dataset:
     """Geocode data from SAR grid to map coordinates.
@@ -1067,110 +1109,198 @@ def geocode_data(
     Args:
         data: Input dataset on the SAR grid.
         time_az: Target azimuth times.
-        time_slr: Target slant-range times.
-        time_slr_gcp: GCP slant-range times.
+        ground_range: Target ground range.
         interp_method: Interpolation method.
 
     Returns:
         Geocoded dataset.
     """
-
-    coeff = fit_ground_range(time_slr_gcp)
-    method = "linear" if interp_method == "bilinear" else "nearest"
-
-    def _interp_block(block):
-        coeff_interp = coeff.interp(azimuth_time=block.time_az)
-        x_tgt = (block.time_slr - coeff.attrs["mean"]) / coeff.attrs["std"]
-        ground_range = (coeff_interp * x_tgt**coeff.degree).sum("degree")
-        return data.interp(
-            azimuth_time=block.time_az,
-            ground_range=ground_range,
-            method=method,
+    az_idx = (src_loc.azimuth_time - grid_params.az0) / _ONE_SECOND / grid_params.d_az
+    gr_idx = (src_loc.ground_range - grid_params.gr0) / grid_params.d_gr
+    scr_indexing = _compute_indexing(data, az_idx, gr_idx)
+    temp_ij_bboxes = scr_indexing.ij_bboxes.copy()
+    temp_ij_bboxes[[1, 3]] -= scr_indexing.pad_width[0][0]
+    temp_ij_bboxes[[0, 2]] -= scr_indexing.pad_width[1][0]
+    tile_size = tuple(chunk[0] for chunk in gr_idx.chunks)
+    for j in range(temp_ij_bboxes.shape[1]):
+        for i in range(temp_ij_bboxes.shape[2]):
+            i_min = tile_size[1] * i
+            i_max = tile_size[1] * (i + 1)
+            j_min = tile_size[0] * j
+            j_max = tile_size[0] * (j + 1)
+            gr_idx[j_min:j_max, i_min:i_max] -= temp_ij_bboxes[0, j, i]
+            az_idx[j_min:j_max, i_min:i_max] -= temp_ij_bboxes[1, j, i]
+
+    target_ds = xr.Dataset(coords=az_idx.coords)
+    for var_name, data_array in data.items():
+        tiled = _reorganize_tiled_array(data_array.data, scr_indexing, np.nan)
+        resampled = da.map_blocks(
+            _sample_array_at_indices,
+            tiled,
+            gr_idx.data,
+            az_idx.data,
+            interp_method=interp_method,
+            dtype=data.dtype,
+            chunks=gr_idx.data.chunks,
         )
+        target_ds[var_name] = (az_idx.dims, resampled)
 
-    # Build template with new coordinates
-    chunksizes = {}
-    for val in [time_az, time_slr]:
-        for dim in val.dims:
-            chunksizes[dim] = val.chunksizes[dim]
-    coeff_interp = coeff.interp(azimuth_time=time_az)
-    x_tgt = (time_slr - coeff.attrs["mean"]) / coeff.attrs["std"]
-    ground_range = (coeff_interp * x_tgt**coeff.degree).sum("degree")
-    template = data.interp(
-        azimuth_time=time_az,
-        ground_range=ground_range,
-    ).chunk(chunksizes)
-
-    target_coords = xr.Dataset({"time_az": time_az, "time_slr": time_slr})
-    return xr.map_blocks(_interp_block, target_coords, template=template)
+    return target_ds
 
 
-def terrain_correct(
-    data: xr.Dataset,
+def get_source_location(
+    dem: xr.DataArray,
     time_slr_gcp: xr.DataArray,
     sat_position: xr.DataArray,
-    dem: xr.DataArray,
-    apply_rtc: bool = True,
-    grid_params: GridParams | None = None,
-    interp_method: Literal["nearest", "bilinear"] = "nearest",
+    grid_params: GridParams,
+    gm_dem: GridMapping,
+    apply_rtc: bool,
 ) -> xr.Dataset:
-    """Apply terrain correction to SAR data.
 
-    Args:
-        data: Input SAR dataset.
-        time_slr_gcp: GCP slant-range times.
-        sat_position: Satellite positions over time.
-        dem: DEM for terrain correction.
-        apply_rtc: Whether to apply radiometric terrain correction.
-        grid_params: Grid parameters for RTC.
-        interp_method: Interpolation method.
-
-    Returns:
-        Terrain-corrected dataset.
-
-    Raises:
-        ValueError: If RTC is enabled without grid parameters.
-    """
-    gm_dem = GridMapping.from_dataset(dem.to_dataset(name="dem"))
-    dem_ecef = convert_dem_to_ecef(dem, gm_dem)
+    # get polynomial coefficients to convert from slant range to ground range
+    gr_coeff = fit_ground_range(time_slr_gcp)
 
-    acquisition = simulate_acquisition(
-        dem_ecef, gm_dem, sat_position, apply_rtc=apply_rtc
-    )
+    # get polynomial coefficient to convert from azimuth time to position and velocity
+    pos_coeff = fit_position(sat_position)
+    vel_coeff = poly_derivative(pos_coeff)
 
-    geocoded = geocode_data(
-        data,
-        acquisition.azimuth_time,
-        acquisition.slant_range_time,
-        time_slr_gcp,
-        interp_method,
+    data_array = xr.zeros_like(dem, dtype="float32").drop_vars("spatial_ref")
+    data_array_az = xr.zeros_like(dem, dtype="datetime64[ns]").drop_vars("spatial_ref")
+    template = xr.Dataset(
+        {"azimuth_time": data_array_az, "ground_range": data_array},
     )
-
     if apply_rtc:
-        if grid_params is None:
-            raise ValueError("grid parameters required for RTC")
-
-        if interp_method == "bilinear":
-            weights_fn = gamma_weights_bilinear
-        else:  # interp_method == "nearest"
-            weights_fn = gamma_weights_nearest
-        beta_sim = apply_gamma_weights(acquisition, weights_fn, grid_params)
-        geocoded = geocoded / beta_sim
-        rename_dict = {
-            name: name.replace("beta0", "gamma0") for name in geocoded.data_vars
-        }
-        geocoded = geocoded.rename(rename_dict)
-        for var in geocoded.data_vars:
-            geocoded[var].attrs.update(
-                long_name="gamma nought backscatter coefficient",
-                units="1",
-            )
-
-    geocoded = assign_grid_mapping(geocoded)
-    return geocoded
+        template["gamma_area"] = data_array
+    gm_dem_params = {
+        "crs": gm_dem.crs.to_wkt(),
+        "xy_var_names": gm_dem.xy_var_names,
+    }
+    out = xr.map_blocks(
+        backward_geocode,
+        dem,
+        kwargs={
+            "pos_coeff": pos_coeff,
+            "vel_coeff": vel_coeff,
+            "gr_coeff": gr_coeff,
+            "grid_params": grid_params,
+            "apply_rtc": apply_rtc,
+            "gm_dem_params": gm_dem_params,
+        },
+        template=template,
+    )
+    out["spatial_ref"] = xr.DataArray(0, attrs=gm_dem.crs.to_cf())
+    return out
 
 
 def assign_grid_mapping(dataset: xr.Dataset) -> xr.Dataset:
     for var_name, data_var in dataset.data_vars.items():
         dataset[var_name].attrs["grid_mapping"] = "spatial_ref"
     return dataset
+
+
+# INTERPOLATION -> xcube-resampling?
+def _xy_bbox_block(x_coords: np.ndarray, y_coords: np.ndarray):
+    x_edges = np.concatenate([x_coords[:, 0], x_coords[:, -1]])
+    y_edges = np.concatenate([y_coords[0, :], y_coords[-1, :]])
+    bbox = np.array(
+        [
+            np.floor(x_edges.min()),
+            np.floor(y_edges.min()),
+            np.ceil(x_edges.max()),
+            np.ceil(y_edges.max()),
+        ],
+        dtype=np.int32,
+    )
+    return bbox[:, None, None]
+
+
+def _compute_indexing(
+    data: xr.Dataset,
+    az_ix: xr.DataArray,
+    gr_idx: xr.DataArray,
+) -> SourceTileIndexing:
+
+    src_ij_bboxes = da.map_blocks(
+        _xy_bbox_block,
+        gr_idx.data,
+        az_ix.data,
+        dtype=gr_idx.dtype,
+        chunks=(4, 1, 1),
+    )
+    src_ij_bboxes = src_ij_bboxes.compute()
+
+    # Extend bounding box indices to match the largest bounding box.
+    # This ensures uniform chunk sizes, which are required for da.map_blocks.
+    i_diff = src_ij_bboxes[2] - src_ij_bboxes[0]
+    j_diff = src_ij_bboxes[3] - src_ij_bboxes[1]
+    i_diff_max = np.nanmax(i_diff) + 1
+    j_diff_max = np.nanmax(j_diff) + 1
+    i_half = (i_diff_max - i_diff) // 2
+    j_half = (j_diff_max - j_diff) // 2
+    src_ij_bboxes[0] -= i_half
+    src_ij_bboxes[2] = src_ij_bboxes[0] + i_diff_max
+    src_ij_bboxes[1] -= j_half
+    src_ij_bboxes[3] = src_ij_bboxes[1] + j_diff_max
+
+    # assign padding if needed
+    i_min = np.nanmin(src_ij_bboxes[0])
+    i_max = np.nanmax(src_ij_bboxes[2])
+    j_min = np.nanmin(src_ij_bboxes[[1, 3]])
+    j_max = np.nanmax(src_ij_bboxes[[1, 3]])
+    pad_width = (
+        (-min(0, int(j_min)), max(0, int(j_max - data.sizes["azimuth_time"]))),
+        (-min(0, int(i_min)), max(0, int(i_max - data.sizes["ground_range"]))),
+    )
+    src_ij_bboxes[[1, 3]] += pad_width[0][0]
+    src_ij_bboxes[[0, 2]] += pad_width[1][0]
+
+    tile_size = (int(j_diff_max), int(i_diff_max))
+    size = (
+        int(j_diff_max * src_ij_bboxes.shape[1]),
+        int(i_diff_max * src_ij_bboxes.shape[2]),
+    )
+
+    return SourceTileIndexing(
+        ij_bboxes=src_ij_bboxes,
+        pad_width=pad_width,
+        output_size=size,
+        tile_size=tile_size,
+    )
+
+
+def _sample_array_at_indices(
+    data: np.ndarray,
+    ix: np.ndarray,
+    iy: np.ndarray,
+    interp_method: Literal["nearest", "bilinear"] | None = None,
+) -> np.ndarray:
+    """
+    Sample a 3d array at fractional indices (iy, ix).
+    """
+    if interp_method == "nearest":
+        ix_i = np.ceil(ix - 0.5).astype(np.intp)
+        iy_i = np.ceil(iy - 0.5).astype(np.intp)
+        return data[iy_i, ix_i]
+
+    ix_floor = np.floor(ix).astype(np.intp)
+    iy_floor = np.floor(iy).astype(np.intp)
+    ix_ceil = np.ceil(ix).astype(np.intp)
+    iy_ceil = np.ceil(iy).astype(np.intp)
+
+    dx = ix - ix_floor
+    dy = iy - iy_floor
+
+    v00 = data[iy_floor, ix_floor]
+    v01 = data[iy_floor, ix_ceil]
+    v10 = data[iy_ceil, ix_floor]
+    v11 = data[iy_ceil, ix_ceil]
+
+    if interp_method == "bilinear":
+        u0 = v00 + dx * (v01 - v00)
+        u1 = v10 + dx * (v11 - v10)
+        return u0 + dy * (u1 - u0)
+
+    raise NotImplementedError(
+        f"interp_methods must be one of 'nearest', 'bilinear', "
+        f"was '{interp_method}'."
+    )

From 97627a818c7c2c5e4ba59d916feea0e1b38c91f9 Mon Sep 17 00:00:00 2001
From: konstntokas <konstantin.ntokas@brockmann-consult.de>
Date: Wed, 8 Jul 2026 10:17:00 +0200
Subject: [PATCH 2/9] backuo

---
 xarray_eopf/amodes/sentinel1.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/xarray_eopf/amodes/sentinel1.py b/xarray_eopf/amodes/sentinel1.py
index 7b5bd53..a2d8350 100644
--- a/xarray_eopf/amodes/sentinel1.py
+++ b/xarray_eopf/amodes/sentinel1.py
@@ -873,6 +873,8 @@ def backward_geocode(
         out["gamma_area"] = compute_gamma_area(
             dem_ecef, gm_dem_params, dist / slant_range
         )
+    print(out["gamma_area"].min().values)
+    print(out["gamma_area"].max().values)
     return out
 
 

From 6bfbeddc57e5618884f4641ae46da9297d2a3b9f Mon Sep 17 00:00:00 2001
From: konstntokas <konstantin.ntokas@brockmann-consult.de>
Date: Wed, 8 Jul 2026 10:20:46 +0200
Subject: [PATCH 3/9] backup

---
 xarray_eopf/amodes/sentinel1.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/xarray_eopf/amodes/sentinel1.py b/xarray_eopf/amodes/sentinel1.py
index a2d8350..9969087 100644
--- a/xarray_eopf/amodes/sentinel1.py
+++ b/xarray_eopf/amodes/sentinel1.py
@@ -301,6 +301,8 @@ def _terrain_correct(
             else:  # interp_method == "nearest"
                 weights_fn = gamma_weights_nearest
             gamma_weights = apply_gamma_weights(src_loc, weights_fn, grid_params)
+            print(gamma_weights.min().values)
+            print(gamma_weights.max().values)
             geocoded /= gamma_weights
             rename_dict = {
                 name: name.replace("beta0", "gamma0") for name in geocoded.data_vars
@@ -873,8 +875,6 @@ def backward_geocode(
         out["gamma_area"] = compute_gamma_area(
             dem_ecef, gm_dem_params, dist / slant_range
         )
-    print(out["gamma_area"].min().values)
-    print(out["gamma_area"].max().values)
     return out
 
 

From 94c70e12a3b70959b5e9552c157a170cc7e8d022 Mon Sep 17 00:00:00 2001
From: konstntokas <konstantin.ntokas@brockmann-consult.de>
Date: Wed, 8 Jul 2026 15:37:10 +0200
Subject: [PATCH 4/9] ready for review

---
 CHANGES.md                        |   7 +
 docs/guide.md                     |   6 +-
 integration/test_sen1_analysis.py |   8 +-
 integration/test_sen1_native.py   |  12 +-
 tests/amodes/test_sentinel1.py    | 526 ++++++++++--------------------
 xarray_eopf/amodes/sentinel1.py   | 148 +++++----
 xarray_eopf/backend.py            |   6 +
 7 files changed, 293 insertions(+), 420 deletions(-)

diff --git a/CHANGES.md b/CHANGES.md
index 7dd06f5..1b7cae4 100644
--- a/CHANGES.md
+++ b/CHANGES.md
@@ -1,9 +1,16 @@
+## Changes in 0.2.10 (in development)
+
+- Sentinel-1 GRD analysis mode is now fully lazy, enabling seamless execution on 
+  local and distributed Dask clusters.
+
+
 ## Changes in 0.2.9 (from 2026-06-03)
 
 - Added support for Sentinel-1 Level-2 OCN analysis mode.
 - Fixed an issue in Sentinel-1 GRD analysis mode that could produce NaN values along 
   the edges of the bounding box.
 
+
 ## Changes in 0.2.8 (from 2026-05-08)
 
 - Fix package discovery in `pyproject.toml` to ensure only `xarray_eopf` 
diff --git a/docs/guide.md b/docs/guide.md
index 35e78a7..7cdf1a9 100644
--- a/docs/guide.md
+++ b/docs/guide.md
@@ -110,7 +110,11 @@ Sentinel-1 Level-1 GRD data is provided in radar geometry, defined by the coordi
   `nearest`, `bilinear`.
 - `footprint_scale_factor`: Defines how radar pixels contribute to the output grid.
   Default: `(3.0, 3.0)`, accounting for resolution differences (e.g., ~10 m GRD
-  vs. ~30 m DEM).
+  vs. ~30 m DEM). 
+- `cache_uri`: Temporary path used to store intermediate results from the
+  backward geocoding step in the Sentinel-1 processing workflow. The cache is 
+  automatically removed when the Python process exits. If None, a temporary 
+  directory with a unique UUID-based name is created.
 
 Examples:  
 
diff --git a/integration/test_sen1_analysis.py b/integration/test_sen1_analysis.py
index 8649046..d22796b 100644
--- a/integration/test_sen1_analysis.py
+++ b/integration/test_sen1_analysis.py
@@ -15,7 +15,7 @@
 show_chunking = False
 
 
-class Sentinel2AnalysisTest(TestCase):
+class Sentinel1AnalysisTest(TestCase):
     def test_open_dataset_sen1_grd(self):
         dem_path = Path(__file__).resolve().parent / "test_data" / "dem_small.zarr.zip"
         store = zarr.ZipStore(str(dem_path), mode="r")
@@ -28,9 +28,9 @@ def test_open_dataset_sen1_grd(self):
         dem = dem.dem
 
         url = (
-            "https://objects.eodc.eu/e05ab01a9d56408d82ac32d69a5aae2a:202603-"
-            "s01siwgrh-global/19/products/cpm_v262/S1A_IW_GRDH_1SDV_20260319"
-            "T102725_20260319T102758_063695_0801D3_2EC6.zarr"
+            "https://objects.eodc.eu/e05ab01a9d56408d82ac32d69a5aae2a:202606"
+            "-s01siwgrh-global/23/products/cpm_v270/S1D_IW_GRDH_1SDV_20260623"
+            "T225558_20260623T225623_003369_005EC5_B4C2.zarr"
         )
         with timeit("open " + url) as result:
             # noinspection PyTypeChecker
diff --git a/integration/test_sen1_native.py b/integration/test_sen1_native.py
index f11a574..abdf482 100644
--- a/integration/test_sen1_native.py
+++ b/integration/test_sen1_native.py
@@ -10,19 +10,19 @@
 class Sentinel1NativeTest(TestCase):
     def test_open_datatree_sen1_grd(self):
         path = (
-            "https://objects.eodc.eu/e05ab01a9d56408d82ac32d69a5aae2a:202603-"
-            "s01siwgrh-global/19/products/cpm_v262/S1A_IW_GRDH_1SDV_20260319"
-            "T102725_20260319T102758_063695_0801D3_2EC6.zarr"
+            "https://objects.eodc.eu/e05ab01a9d56408d82ac32d69a5aae2a:202606"
+            "-s01siwgrh-global/23/products/cpm_v270/S1D_IW_GRDH_1SDV_20260623"
+            "T225558_20260623T225623_003369_005EC5_B4C2.zarr"
         )
         # noinspection PyTypeChecker
         dt = xr.open_datatree(path, engine="eopf-zarr", op_mode="native")
         self.assertEqual(25, len(dt.groups))
         self.assertIn(
-            "/S01SIWGRD_20260319T102725_0033_A364_2EC6_0801D3_VH/measurements",
+            "/S01SIWGRD_20260623T225558_0025_D019_B4C2_005EC5_VH/measurements",
             dt.groups,
         )
-        ds = dt.S01SIWGRD_20260319T102725_0033_A364_2EC6_0801D3_VH.measurements
-        self.assertEqual({"azimuth_time": 22290, "ground_range": 25223}, ds.sizes)
+        ds = dt.S01SIWGRD_20260623T225558_0025_D019_B4C2_005EC5_VH.measurements
+        self.assertEqual({"azimuth_time": 16802, "ground_range": 25319}, ds.sizes)
 
     def test_open_datatree_sen1_slc(self):
         path = (
diff --git a/tests/amodes/test_sentinel1.py b/tests/amodes/test_sentinel1.py
index 7e23ebe..3c88005 100644
--- a/tests/amodes/test_sentinel1.py
+++ b/tests/amodes/test_sentinel1.py
@@ -17,7 +17,7 @@
 from tests.helpers import make_s1_grd_datatree, make_s1_ocn_datatree
 from xarray_eopf.amode import AnalysisModeRegistry
 from xarray_eopf.amodes import sentinel1 as sen1
-from xarray_eopf.amodes.sentinel1 import Sen1GRD, Sen1OCN, register
+from xarray_eopf.amodes.sentinel1 import Sen1GRD, Sen1OCN, register, GridParams
 
 
 class Sentinel1AnalysisModeTest(TestCase):
@@ -68,12 +68,16 @@ def test_is_not_valid_source(self):
     def test_get_grid_parameters(self):
         params = self.mode._get_grid_parameters(self.dt, (2.0, 3.0))
 
-        self.assertEqual(1.0e-4, params["slr0"])
-        self.assertEqual(30.0, params["spacing_slr"])
-        self.assertAlmostEqual(30.0 * 2.0 / sen1._SPEED_OF_LIGHT, params["d_slr"])
+        self.assertEqual(0.0, params["gr0"])
+        self.assertEqual(20.0, params["spacing_az"])
+        self.assertEqual(10.0, params["d_gr"])
+        self.assertEqual(30.0, params["d_gr_scale"])
         self.assertEqual(np.datetime64("2024-01-01T00:00:00"), params["az0"])
-        self.assertEqual(1.0, params["d_az"])
-        self.assertEqual(40.0, params["spacing_az"])
+        self.assertEqual(0.5, params["d_az"])
+        self.assertEqual(40.0, params["spacing_az_scale"])
+        self.assertEqual(
+            np.datetime64("2024-01-01T00:00:00.250000000"), params["az0_scale"]
+        )
 
     def test_get_applicable_params(self: TestCase):
         dem = xr.DataArray(np.ones((2, 2)), dims=("lat", "lon"))
@@ -108,7 +112,9 @@ def test_convert_datatree(self):
             {"vv": xr.DataArray(np.ones((2, 2)), dims=("lat", "lon"))}
         )
 
-        with patch.object(sen1, "terrain_correct", return_value=expected) as mocked:
+        with patch.object(
+            self.mode, "_terrain_correct", return_value=expected
+        ) as mocked:
             out = self.mode.convert_datatree(self.dt, includes=["vv"], dem=self.dem)
 
         self.assertIs(out, expected)
@@ -117,7 +123,7 @@ def test_convert_datatree(self):
         self.assertIs(args[3], self.dem)
         self.assertEqual("bilinear", kwargs["interp_method"])
         self.assertTrue(kwargs["apply_rtc"])
-        self.assertIn("slr0", kwargs["grid_params"])
+        self.assertIn("gr0", args[4])
 
     def test_convert_datatree_uses_get_dem(self):
         expected = xr.Dataset(
@@ -125,7 +131,7 @@ def test_convert_datatree_uses_get_dem(self):
         )
 
         with patch.object(sen1, "get_dem", return_value=self.dem) as get_dem_mock:
-            with patch.object(sen1, "terrain_correct", return_value=expected):
+            with patch.object(self.mode, "_terrain_correct", return_value=expected):
                 _ = self.mode.convert_datatree(self.dt, includes=["vv"])
 
         get_dem_mock.assert_called_once()
@@ -227,41 +233,69 @@ def test_convert_datatree_fail(self):
 
 
 class Sentinel1FunctionsTest(TestCase):
-    def test_gridparams_iter(self):
-        params = sen1.GridParams(
-            slr0=0.0,
-            d_slr=1.0,
-            spacing_slr=2.0,
+
+    def setUp(self):
+        self.gm_dem_params = {
+            "crs": "EPSG:4326",
+            "xy_var_names": ("lat", "lon"),
+        }
+        self.grid_params = sen1.GridParams(
+            gr0=0.0,
+            gr0_scale=0.0,
+            d_gr=1.0,
+            d_gr_scale=1.0,
             az0=np.datetime64("2024-01-01T00:00:00"),
+            az0_scale=np.datetime64("2024-01-01T00:00:00"),
             d_az=1.0,
-            spacing_az=3.0,
+            d_az_scale=1.0,
+            spacing_az=2.0,
+            spacing_az_scale=2.0,
         )
-        self.assertEqual(
-            ["slr0", "d_slr", "spacing_slr", "az0", "d_az", "spacing_az"],
-            list(iter(params)),
+        self.dem = xr.DataArray(
+            np.ones((2, 2), dtype="float64"),
+            dims=("lat", "lon"),
+            coords={"lat": [0.0, 0.1], "lon": [0.0, 0.1]},
         )
-
-    def test_acquisition_getitem_raises_for_missing_gamma_area(self):
-        acquisition = sen1.Acquisition(
-            azimuth_time=xr.DataArray(np.array([0]), dims=("lat",)),
-            distance=xr.DataArray(np.ones((3, 1, 1)), dims=("axis", "lat", "lon")),
-            velocity=xr.DataArray(np.ones((3, 1, 1)), dims=("axis", "lat", "lon")),
-            slant_range_time=xr.DataArray(np.array([0.0]), dims=("lon",)),
-            gamma_area=None,
+        self.dem_ecef = xr.DataArray(
+            np.ones((3, 2, 2), dtype="float64"),
+            dims=("axis", "lat", "lon"),
+            coords={"axis": ["x", "y", "z"], "lat": [0.0, 0.1], "lon": [0.0, 0.1]},
+        )
+        self.posvel_coeff = xr.DataArray(
+            np.zeros((2, 3)),
+            dims=("degree", "axis"),
+            coords={"degree": [1, 0], "axis": ["x", "y", "z"]},
+            attrs={"epoch": np.datetime64("2024-01-01T00:00:00")},
+        )
+        self.gr_coeff = xr.DataArray(
+            np.zeros((2, 9)),
+            dims=("azimuth_time", "degree"),
+            coords={
+                "degree": np.arange(8, -1, -1),
+                "azimuth_time": np.array(
+                    ["2024-01-01T00:00:00", "2024-01-01T00:00:02"],
+                    dtype="datetime64[ns]",
+                ),
+            },
+            attrs=dict(mean=1, std=1),
         )
-        with pytest.raises(KeyError, match="gamma_area"):
-            _ = acquisition["gamma_area"]
 
-    def test_acquisition_getitem_returns_existing_field(self):
-        distance = xr.DataArray(np.ones((3, 1, 1)), dims=("axis", "lat", "lon"))
-        acquisition = sen1.Acquisition(
-            azimuth_time=xr.DataArray(np.array([0]), dims=("lat",)),
-            distance=distance,
-            velocity=xr.DataArray(np.ones((3, 1, 1)), dims=("axis", "lat", "lon")),
-            slant_range_time=xr.DataArray(np.array([0.0]), dims=("lon",)),
-            gamma_area=None,
+    def test_gridparams_iter(self):
+        self.assertEqual(
+            [
+                "gr0",
+                "gr0_scale",
+                "d_gr",
+                "d_gr_scale",
+                "az0",
+                "az0_scale",
+                "d_az",
+                "d_az_scale",
+                "spacing_az",
+                "spacing_az_scale",
+            ],
+            list(iter(self.grid_params)),
         )
-        self.assertIs(acquisition["distance"], distance)
 
     def test_get_dem_requires_credentials(self):
         with patch.dict(os.environ, {}, clear=True):
@@ -392,7 +426,10 @@ def test_convert_dem_to_ecef(self):
             coords={"lat": [0, 1, 2, 3], "lon": [0, 1, 2, 3]},
         )
         gm_dem = GridMapping.from_dataset(dem.to_dataset(name="dem"))
-        out = sen1.convert_dem_to_ecef(dem, gm_dem)
+        out = sen1.convert_dem_to_ecef(
+            dem,
+            {"crs": gm_dem.crs.to_wkt(), "xy_var_names": gm_dem.xy_var_names},
+        )
         self.assertEqual(("axis", "lat", "lon"), out.dims)
         self.assertEqual(3, out.sizes["axis"])
 
@@ -482,40 +519,28 @@ def func_p(t, _payload):
         self.assertEqual(0, k)
 
     def test_backward_geocode_invalid_method(self):
-        dem_ecef = xr.DataArray(
-            np.ones((3, 1, 1), dtype="float64"),
-            dims=("axis", "lat", "lon"),
-            coords={"axis": ["x", "y", "z"], "lat": [0.0], "lon": [0.0]},
-        )
-        coeff = xr.DataArray(
-            np.zeros((2, 3)),
-            dims=("degree", "axis"),
-            coords={"degree": [1, 0], "axis": ["x", "y", "z"]},
-        )
         with pytest.raises(ValueError, match="method needs to be either"):
-            sen1.backward_geocode(dem_ecef, coeff, coeff, method="x")
+            sen1.backward_geocode(
+                self.dem,
+                pos_coeff=self.posvel_coeff,
+                vel_coeff=self.posvel_coeff,
+                gr_coeff=self.gr_coeff,
+                gm_dem_params=self.gm_dem_params,
+                grid_params=self.grid_params,
+                method="x",
+            )
 
     def test_backward_geocode_secant_and_newton_paths(self):
-        dem_ecef = xr.DataArray(
-            np.ones((3, 1, 1), dtype="float64"),
-            dims=("axis", "lat", "lon"),
-            coords={"axis": ["x", "y", "z"], "lat": [0.0], "lon": [0.0]},
-        )
-        coeff = xr.DataArray(
-            np.zeros((2, 3)),
-            dims=("degree", "axis"),
-            coords={"degree": [1, 0], "axis": ["x", "y", "z"]},
-        )
         payload = (
-            xr.DataArray(np.ones((3, 1, 1)), dims=("axis", "lat", "lon")),
-            xr.DataArray(np.ones((3, 1, 1)), dims=("axis", "lat", "lon")),
+            xr.DataArray(np.ones((3, 2, 2)), dims=("axis", "lat", "lon")),
+            xr.DataArray(np.ones((3, 2, 2)), dims=("axis", "lat", "lon")),
         )
         with (
             patch.object(
                 sen1,
                 "secant",
                 return_value=(
-                    xr.DataArray([[0.0]], dims=("lat", "lon")),
+                    xr.DataArray([[0.0, 0.0], [0.0, 0.0]], dims=("lat", "lon")),
                     None,
                     None,
                     0,
@@ -526,88 +551,33 @@ def test_backward_geocode_secant_and_newton_paths(self):
                 sen1,
                 "newton",
                 return_value=(
-                    xr.DataArray([[0.0]], dims=("lat", "lon")),
+                    xr.DataArray([[0.0, 0.0], [0.0, 0.0]], dims=("lat", "lon")),
                     None,
                     0,
                     payload,
                 ),
             ),
         ):
-            out_secant = sen1.backward_geocode(dem_ecef, coeff, coeff, method="secant")
-            out_newton = sen1.backward_geocode(dem_ecef, coeff, coeff, method="newton")
-        self.assertEqual(3, len(out_secant))
-        self.assertEqual(3, len(out_newton))
-
-    def test_simulate_acquisition_without_rtc(self):
-        dem_ecef = xr.DataArray(
-            np.ones((3, 2, 2), dtype="float64"),
-            dims=("axis", "lat", "lon"),
-            coords={"axis": ["x", "y", "z"], "lat": [0.0, 1.0], "lon": [0.0, 1.0]},
-        )
-        gm_dem = GridMapping.from_dataset(dem_ecef.to_dataset(name="dem"))
-        sat_position = xr.DataArray(
-            np.zeros((2, 3)),
-            dims=("azimuth_time", "axis"),
-            coords={
-                "azimuth_time": np.array(
-                    [
-                        np.datetime64("2024-01-01T00:00:00"),
-                        np.datetime64("2024-01-01T00:00:01"),
-                    ]
-                ),
-                "axis": ["x", "y", "z"],
-            },
-        )
-        dist = xr.DataArray(
-            np.ones((3, 2, 2), dtype="float64"),
-            dims=("axis", "lat", "lon"),
-            coords={"axis": ["x", "y", "z"], "lat": [0.0, 1.0], "lon": [0.0, 1.0]},
-        )
-        vel = dist.copy()
-        t = xr.DataArray(np.zeros((2, 2), dtype="float64"), dims=("lat", "lon"))
-        with patch.object(sen1, "backward_geocode", return_value=(t, dist, vel)):
-            acq = sen1.simulate_acquisition(
-                dem_ecef, gm_dem, sat_position, apply_rtc=False
+            out_secant = sen1.backward_geocode(
+                self.dem,
+                pos_coeff=self.posvel_coeff,
+                vel_coeff=self.posvel_coeff,
+                gr_coeff=self.gr_coeff,
+                grid_params=self.grid_params,
+                gm_dem_params=self.gm_dem_params,
+                method="secant",
             )
-        self.assertIn("slant_range_time", acq)
-        self.assertNotIn("gamma_area", acq)
-
-    def test_simulate_acquisition_with_rtc(self):
-        dem_ecef = xr.DataArray(
-            np.ones((3, 2, 2), dtype="float64"),
-            dims=("axis", "lat", "lon"),
-            coords={"axis": ["x", "y", "z"], "lat": [0.0, 1.0], "lon": [0.0, 1.0]},
-        )
-        gm_dem = GridMapping.from_dataset(dem_ecef.to_dataset(name="dem"))
-        sat_position = xr.DataArray(
-            np.zeros((2, 3)),
-            dims=("azimuth_time", "axis"),
-            coords={
-                "azimuth_time": np.array(
-                    [
-                        np.datetime64("2024-01-01T00:00:00"),
-                        np.datetime64("2024-01-01T00:00:01"),
-                    ]
-                ),
-                "axis": ["x", "y", "z"],
-            },
-        )
-        dist = xr.DataArray(
-            np.ones((3, 2, 2), dtype="float64"),
-            dims=("axis", "lat", "lon"),
-            coords={"axis": ["x", "y", "z"], "lat": [0.0, 1.0], "lon": [0.0, 1.0]},
-        )
-        vel = dist.copy()
-        t = xr.DataArray(np.zeros((2, 2), dtype="float64"), dims=("lat", "lon"))
-        gamma = xr.DataArray(np.ones((2, 2)), dims=("lat", "lon"))
-        with (
-            patch.object(sen1, "backward_geocode", return_value=(t, dist, vel)),
-            patch.object(sen1, "compute_gamma_area", return_value=gamma),
-        ):
-            acq = sen1.simulate_acquisition(
-                dem_ecef, gm_dem, sat_position, apply_rtc=True
+            out_newton = sen1.backward_geocode(
+                self.dem,
+                pos_coeff=self.posvel_coeff,
+                vel_coeff=self.posvel_coeff,
+                gr_coeff=self.gr_coeff,
+                grid_params=self.grid_params,
+                gm_dem_params=self.gm_dem_params,
+                method="newton",
             )
-        self.assertIn("gamma_area", acq)
+        self.assertEqual(3, len(out_secant))
+        self.assertEqual(3, len(out_newton))
 
     def test_compute_gamma_area_clips_negative(self):
         dem_ecef = xr.DataArray(
@@ -615,7 +585,12 @@ def test_compute_gamma_area_clips_negative(self):
             dims=("axis", "lat", "lon"),
             coords={"axis": ["x", "y", "z"], "lat": [0.0, 1.0], "lon": [0.0, 1.0]},
         )
-        gm_dem = GridMapping.from_dataset(dem_ecef.to_dataset(name="dem"))
+        gm_dem = {
+            "crs": GridMapping.from_dataset(
+                dem_ecef.to_dataset(name="dem")
+            ).crs.to_wkt(),
+            "xy_var_names": ("lon", "lat"),
+        }
         area = xr.DataArray(
             np.array(
                 [
@@ -651,14 +626,11 @@ def test_compute_dem_area(self):
         y = 20.0 + lat2d
         z = 30.0 + lon2d + lat2d
         dem_ecef = xr.DataArray(
-            da.from_array(
-                np.stack([x, y, z], axis=0).astype("float32"), chunks=(3, 2, 2)
-            ),
+            np.stack([x, y, z], axis=0).astype("float32"),
             dims=("axis", "lat", "lon"),
             coords={"axis": ["x", "y", "z"], "lat": lat, "lon": lon},
         )
-        gm_dem = GridMapping.from_dataset(dem_ecef.to_dataset(name="dem"))
-        area = sen1.compute_dem_area(dem_ecef, gm_dem)
+        area = sen1.compute_dem_area(dem_ecef, self.gm_dem_params)
         self.assertEqual(("axis", "lat", "lon"), area.dims)
 
     def test_sum_weights_and_gamma_weight_helpers(self):
@@ -666,16 +638,16 @@ def test_sum_weights_and_gamma_weight_helpers(self):
             {
                 "gamma_area": xr.DataArray([[1.0, 2.0]], dims=("lat", "lon")),
                 "az_idx": xr.DataArray([[0.2, 0.8]], dims=("lat", "lon")),
-                "slr_idx": xr.DataArray([[1.2, 1.8]], dims=("lat", "lon")),
+                "gr_idx": xr.DataArray([[1.2, 1.8]], dims=("lat", "lon")),
             }
         )
         reduced = xr.DataArray(
             [[3.0]],
-            dims=("slr_idx", "az_idx"),
-            coords={"slr_idx": [1], "az_idx": [1]},
+            dims=("gr_idx", "az_idx"),
+            coords={"gr_idx": [1], "az_idx": [1]},
         )
         with patch.object(sen1.flox.xarray, "xarray_reduce", return_value=reduced):
-            summed = sen1.sum_weights(acq.gamma_area, acq.az_idx, acq.slr_idx)
+            summed = sen1.sum_weights(acq.gamma_area, acq.az_idx, acq.gr_idx)
         self.assertEqual(("lat", "lon"), summed.dims)
         self.assertEqual((1, 2), summed.data.shape)
 
@@ -691,39 +663,42 @@ def test_sum_weights_and_gamma_weight_helpers(self):
             sw.assert_called_once()
 
     def test_apply_gamma_weights(self):
-        azimuth_time = xr.DataArray(
-            np.array(
-                ["2024-01-01T00:00:00", "2024-01-01T00:00:01"], dtype="datetime64[ns]"
-            ),
-            dims=("lat",),
-        )
-        slant_range_time = xr.DataArray(np.array([0.0, 2.0]), dims=("lon",))
-        distance = xr.DataArray(
-            np.ones((3, 2, 2), dtype="float64"),
-            dims=("axis", "lat", "lon"),
-            coords={"axis": ["x", "y", "z"]},
+        params = sen1.GridParams(
+            gr0=0.0,
+            gr0_scale=0.0,
+            d_gr=1.0,
+            d_gr_scale=1.0,
+            az0=np.datetime64("2024-01-01T00:00:00"),
+            az0_scale=np.datetime64("2024-01-01T00:00:00"),
+            d_az=1.0,
+            d_az_scale=1.0,
+            spacing_az=2.0,
+            spacing_az_scale=2.0,
         )
-        acq = sen1.Acquisition(
-            azimuth_time=azimuth_time,
-            distance=distance,
-            velocity=distance,
-            slant_range_time=slant_range_time,
-            gamma_area=xr.DataArray(np.ones((2, 2)), dims=("lat", "lon")),
+        src_loc = xr.Dataset(
+            {
+                "azimuth_time": xr.DataArray(
+                    np.array(
+                        [
+                            ["2024-01-01T00:00:00", "2024-01-01T00:00:01"],
+                            ["2024-01-01T00:00:00", "2024-01-01T00:00:01"],
+                        ],
+                        dtype="datetime64[ns]",
+                    ),
+                    dims=("lat", "lon"),
+                ),
+                "ground_range": xr.DataArray(
+                    np.array([[0.0, 1.0], [0.0, 1.0]]), dims=("lat", "lon")
+                ),
+                "gamma_area": xr.DataArray(np.ones((2, 2)), dims=("lat", "lon")),
+            }
         )
 
         def passthrough(ds):
             return ds.gamma_area * 2
 
-        params = sen1.GridParams(
-            az0=azimuth_time.values[0],
-            d_az=1.0,
-            slr0=0.0,
-            d_slr=1.0,
-            spacing_slr=2.0,
-            spacing_az=2.0,
-        )
-        out = sen1.apply_gamma_weights(acq, passthrough, params)
-        self.assertTrue(np.allclose(out.values, 0.5))
+        out = sen1.apply_gamma_weights(src_loc, passthrough, params)
+        self.assertTrue(np.allclose(out.values, 1.0))
 
     def test_fit_ground_range(self):
         time_slr_gcp = xr.DataArray(
@@ -745,188 +720,49 @@ def test_geocode_data(self):
         data = xr.Dataset(
             {
                 "vv": xr.DataArray(
-                    np.ones((2, 3)), dims=("azimuth_time", "ground_range")
+                    da.ones((2, 3), chunks=(2, 3)),
+                    dims=("azimuth_time", "ground_range"),
                 )
             },
-            coords={"azimuth_time": [0, 1], "ground_range": [0, 1, 2]},
+            coords={
+                "azimuth_time": np.array(
+                    ["2023-12-31T23:59:50", "2024-01-01T00:00:20"],
+                    dtype="datetime64[ns]",
+                ),
+                "ground_range": [0, 3, 6],
+            },
         )
         time_az = xr.DataArray(
-            da.from_array(np.array([[0, 0], [1, 1]]), chunks=(2, 2)),
+            da.from_array(
+                np.array(
+                    [
+                        ["2024-01-01T00:00:00", "2024-01-01T00:00:02"],
+                        ["2024-01-01T00:00:04", "2024-01-01T00:00:06"],
+                    ],
+                    dtype="datetime64[ns]",
+                ),
+                chunks=(2, 2),
+            ),
             dims=("lat", "lon"),
         )
-        time_slr = xr.DataArray(
+        ground_range = xr.DataArray(
             da.from_array(np.array([[2, 3], [2, 3]]), chunks=(2, 2)),
             dims=("lat", "lon"),
         )
-        time_slr_gcp = xr.DataArray(
-            np.array([[1.0, 2.0, 3.0], [1.0, 2.0, 3.0]]),
-            dims=("azimuth_time", "ground_range"),
-            coords={"azimuth_time": [0, 1], "ground_range": [0, 1, 2]},
+        src_loc = xr.Dataset({"azimuth_time": time_az, "ground_range": ground_range})
+        grid_params = sen1.GridParams(
+            gr0=0.0,
+            gr0_scale=0.0,
+            d_gr=3.0,
+            d_gr_scale=3.0,
+            az0=np.datetime64("2023-12-31T23:59:50"),
+            az0_scale=np.datetime64("2023-12-31T23:59:50"),
+            d_az=20.0,
+            d_az_scale=20.0,
+            spacing_az=3.0,
+            spacing_az_scale=3.0,
         )
 
-        out = sen1.geocode_data(data, time_az, time_slr, time_slr_gcp, "nearest")
+        out = sen1.geocode_data(data, src_loc, grid_params, "nearest")
         self.assertIn("vv", out.data_vars)
         np.testing.assert_allclose(out.vv.values, np.ones((2, 2), dtype=float))
-
-    def test_terrain_correct_paths(self):
-        data = xr.Dataset(
-            {"vv": xr.DataArray(np.ones((2, 2)), dims=("azimuth_time", "ground_range"))}
-        )
-        time_slr_gcp = xr.DataArray(
-            np.ones((2, 2)),
-            dims=("azimuth_time", "ground_range"),
-            coords={"azimuth_time": [0, 1], "ground_range": [0, 1]},
-        )
-        sat_position = xr.DataArray(
-            np.ones((2, 3)),
-            dims=("azimuth_time", "axis"),
-            coords={"azimuth_time": [0, 1], "axis": ["x", "y", "z"]},
-        )
-        dem = xr.DataArray(
-            np.ones((2, 2)),
-            dims=("lat", "lon"),
-            coords={"lat": [0, 1], "lon": [0, 1]},
-        )
-        distance = xr.DataArray(
-            np.ones((3, 2, 2), dtype="float64"),
-            dims=("axis", "lat", "lon"),
-            coords={"axis": ["x", "y", "z"]},
-        )
-        acquisition = sen1.Acquisition(
-            azimuth_time=xr.DataArray(np.array([0, 1]), dims=("lat",)),
-            slant_range_time=xr.DataArray(np.array([0.0, 1.0]), dims=("lon",)),
-            distance=distance,
-            velocity=distance,
-            gamma_area=xr.DataArray(np.ones((2, 2)), dims=("lat", "lon")),
-        )
-        geocoded = xr.Dataset(
-            {"vv": xr.DataArray(np.full((2, 2), 2.0), dims=("lat", "lon"))}
-        )
-        with (
-            patch.object(sen1, "convert_dem_to_ecef", return_value=dem),
-            patch.object(
-                sen1,
-                "fit_position",
-                return_value=xr.DataArray(
-                    np.ones((2, 3)),
-                    dims=("degree", "axis"),
-                    coords={"degree": [1, 0], "axis": ["x", "y", "z"]},
-                ),
-            ),
-            patch.object(
-                sen1,
-                "poly_derivative",
-                return_value=xr.DataArray(
-                    np.ones((1, 3)),
-                    dims=("degree", "axis"),
-                    coords={"degree": [0], "axis": ["x", "y", "z"]},
-                ),
-            ),
-            patch.object(sen1, "simulate_acquisition", return_value=acquisition),
-            patch.object(sen1, "geocode_data", return_value=geocoded),
-        ):
-            out = sen1.terrain_correct(
-                data, time_slr_gcp, sat_position, dem, apply_rtc=False
-            )
-            self.assertIn("vv", out)
-            self.assertEqual(2, out.vv.shape[0])
-            self.assertEqual(2, out.vv.values[0, 0])
-            with pytest.raises(ValueError, match="grid parameters required for RTC"):
-                sen1.terrain_correct(
-                    data, time_slr_gcp, sat_position, dem, apply_rtc=True
-                )
-
-        with (
-            patch.object(sen1, "convert_dem_to_ecef", return_value=dem),
-            patch.object(
-                sen1,
-                "fit_position",
-                return_value=xr.DataArray(
-                    np.ones((2, 3)),
-                    dims=("degree", "axis"),
-                    coords={"degree": [1, 0], "axis": ["x", "y", "z"]},
-                ),
-            ),
-            patch.object(
-                sen1,
-                "poly_derivative",
-                return_value=xr.DataArray(
-                    np.ones((1, 3)),
-                    dims=("degree", "axis"),
-                    coords={"degree": [0], "axis": ["x", "y", "z"]},
-                ),
-            ),
-            patch.object(sen1, "simulate_acquisition", return_value=acquisition),
-            patch.object(sen1, "geocode_data", return_value=geocoded),
-            patch.object(
-                sen1,
-                "apply_gamma_weights",
-                return_value=xr.DataArray(np.ones((2, 2)), dims=("lat", "lon")),
-            ) as agw,
-        ):
-            out = sen1.terrain_correct(
-                data,
-                time_slr_gcp,
-                sat_position,
-                dem,
-                apply_rtc=True,
-                grid_params=sen1.GridParams(
-                    slr0=0.0,
-                    d_slr=1.0,
-                    spacing_slr=1.0,
-                    az0=np.datetime64("2024-01-01T00:00:00"),
-                    d_az=1.0,
-                    spacing_az=1.0,
-                ),
-                interp_method="bilinear",
-            )
-            self.assertIn("vv", out)
-            args, _ = agw.call_args
-            self.assertIs(args[1], sen1.gamma_weights_bilinear)
-
-        with (
-            patch.object(sen1, "convert_dem_to_ecef", return_value=dem),
-            patch.object(
-                sen1,
-                "fit_position",
-                return_value=xr.DataArray(
-                    np.ones((2, 3)),
-                    dims=("degree", "axis"),
-                    coords={"degree": [1, 0], "axis": ["x", "y", "z"]},
-                ),
-            ),
-            patch.object(
-                sen1,
-                "poly_derivative",
-                return_value=xr.DataArray(
-                    np.ones((1, 3)),
-                    dims=("degree", "axis"),
-                    coords={"degree": [0], "axis": ["x", "y", "z"]},
-                ),
-            ),
-            patch.object(sen1, "simulate_acquisition", return_value=acquisition),
-            patch.object(sen1, "geocode_data", return_value=geocoded),
-            patch.object(
-                sen1,
-                "apply_gamma_weights",
-                return_value=xr.DataArray(np.ones((2, 2)), dims=("lat", "lon")),
-            ) as agw,
-        ):
-            _ = sen1.terrain_correct(
-                data,
-                time_slr_gcp,
-                sat_position,
-                dem,
-                apply_rtc=True,
-                grid_params=sen1.GridParams(
-                    slr0=0.0,
-                    d_slr=1.0,
-                    spacing_slr=1.0,
-                    az0=np.datetime64("2024-01-01T00:00:00"),
-                    d_az=1.0,
-                    spacing_az=1.0,
-                ),
-                interp_method="nearest",
-            )
-            args, _ = agw.call_args
-            self.assertIs(args[1], sen1.gamma_weights_nearest)
diff --git a/xarray_eopf/amodes/sentinel1.py b/xarray_eopf/amodes/sentinel1.py
index 9969087..5b9cc5f 100644
--- a/xarray_eopf/amodes/sentinel1.py
+++ b/xarray_eopf/amodes/sentinel1.py
@@ -2,20 +2,20 @@
 #  Permissions are hereby granted under the terms of the Apache 2.0 License:
 #  https://opensource.org/license/apache-2-0.
 
-
+import atexit
 import functools
 import os
 import re
+import uuid
 import warnings
 from abc import ABC
 from collections.abc import Iterable, Sequence
 from dataclasses import dataclass, fields
 from typing import Any, Callable, Literal
-import uuid
 
 import dask.array as da
-import fsspec
 import flox.xarray
+import fsspec
 import numpy as np
 import pyproj
 import pystac_client
@@ -25,12 +25,13 @@
 from xcube_resampling.constants import SpatialAggMethods, SpatialInterpMethods
 from xcube_resampling.gridmapping import GridMapping
 from xcube_resampling.rectify import rectify_dataset
+
+# noinspection PyProtectedMember
 from xcube_resampling.utils import (
-    reproject_bbox,
-    resolution_degrees_to_meters,
+    SourceTileIndexing,
     _reorganize_tiled_array,
+    reproject_bbox,
     transform_resolution,
-    SourceTileIndexing,
 )
 
 from xarray_eopf.amode import AnalysisMode, AnalysisModeRegistry
@@ -51,9 +52,11 @@ class GridParams:
     """RTC grid parameters."""
 
     gr0: float
+    gr0_scale: float
     d_gr: float
     d_gr_scale: float
     az0: np.datetime64
+    az0_scale: np.datetime64
     d_az: float
     d_az_scale: float
     spacing_az: float
@@ -97,6 +100,7 @@ class Sen1GRD(Sen1):
     product_type = "GRDH"
     cache_fs: fsspec.AbstractFileSystem | None = None
     cache_uri: str | None = None
+    _cleanup_registered: bool = False
 
     def get_applicable_params(self, **kwargs) -> dict[str, Any]:
         params = {}
@@ -175,11 +179,14 @@ def convert_datatree(
 
         if cache_uri is None:
             self.cache_fs = fsspec.filesystem("file")
-            self.cache_uri = f"tmp/{uuid.uuid4().hex}"
+            self.cache_uri = f"tmp_{uuid.uuid4().hex}"
         else:
             cache_uri = cache_uri.rstrip("/")
             self.cache_fs, _ = fsspec.url_to_fs(cache_uri)
             self.cache_uri = cache_uri
+        if not getattr(self, "_cleanup_registered", False):
+            atexit.register(self._cleanup)
+            self._cleanup_registered = True
 
         # get dem data array
         if dem is None:
@@ -301,11 +308,10 @@ def _terrain_correct(
             else:  # interp_method == "nearest"
                 weights_fn = gamma_weights_nearest
             gamma_weights = apply_gamma_weights(src_loc, weights_fn, grid_params)
-            print(gamma_weights.min().values)
-            print(gamma_weights.max().values)
             geocoded /= gamma_weights
             rename_dict = {
-                name: name.replace("beta0", "gamma0") for name in geocoded.data_vars
+                name: str(name).replace("beta0", "gamma0")
+                for name in geocoded.data_vars
             }
             geocoded = geocoded.rename(rename_dict)
             for var in geocoded.data_vars:
@@ -327,33 +333,42 @@ def _cleanup(self):
 
     @staticmethod
     def _get_grid_parameters(
-        dt: xr.DataTree,
+        datatree: xr.DataTree,
         footprint_scale_factor: tuple[float, float],
     ) -> GridParams:
         """Build grid parameters for RTC from Sentinel-1 metadata.
 
         Args:
-            dt: Source data tree.
+            datatree: Source data tree.
             footprint_scale_factor: Scaling for SAR footprint spacing.
 
         Returns:
             Grid parameters for terrain correction.
         """
 
-        group_vh = [x for x in dt.children if "VH" in x][0]
-        attrs = dt[f"{group_vh}"].attrs["other_metadata"]["image_annotation"][
+        group_vh = [x for x in datatree.children if "VH" in x][0]
+        attrs = datatree[f"{group_vh}"].attrs["other_metadata"]["image_annotation"][
             "image_information"
         ]
 
+        az_scale, gr_scale = footprint_scale_factor
+        gr0 = 0.0
+        d_gr = attrs["range_pixel_spacing"]
+        az0 = np.datetime64(attrs["product_first_line_utc_time"])
+        d_az = attrs["azimuth_time_interval"]
+        spacing_az = attrs["azimuth_pixel_spacing"]
+
         return GridParams(
-            gr0=0.0,
-            d_gr=attrs["range_pixel_spacing"],
-            d_gr_scale=attrs["range_pixel_spacing"] * footprint_scale_factor[1],
-            az0=np.datetime64(attrs["product_first_line_utc_time"]),
-            d_az=attrs["azimuth_time_interval"],
-            spacing_az=attrs["azimuth_pixel_spacing"],
-            d_az_scale=attrs["azimuth_time_interval"] * footprint_scale_factor[0],
-            spacing_az_scale=attrs["azimuth_pixel_spacing"] * footprint_scale_factor[0],
+            gr0=gr0,
+            d_gr=d_gr,
+            d_gr_scale=d_gr * gr_scale,
+            gr0_scale=(gr0 - (0.5 * d_gr) + (0.5 * d_gr * gr_scale)),
+            az0=az0,
+            d_az=d_az,
+            spacing_az=spacing_az,
+            d_az_scale=d_az * az_scale,
+            az0_scale=(az0 + (-(0.5 * d_az) + (0.5 * d_az * az_scale)) * _ONE_SECOND),
+            spacing_az_scale=spacing_az * az_scale,
         )
 
 
@@ -577,7 +592,7 @@ def convert_dem_to_ecef(dem: xr.DataArray, gm_dem_params: dict) -> xr.DataArray:
 
     Args:
         dem: DEM data array.
-        gm_dem: GridMapping of the DEM data array.
+        gm_dem_params: GridMapping metadata of the DEM data array.
 
     Returns:
         DEM expressed in ECEF axes.
@@ -815,9 +830,13 @@ def backward_geocode(
     """Compute orbit time and vectors for a DEM using inverse geocoding.
 
     Args:
-        dem_ecef: DEM in ECEF coordinates.
+        dem: Digital elevation model.
         pos_coeff: Position polynomial coefficients.
         vel_coeff: Velocity polynomial coefficients.
+        gr_coeff: Ground-range polynomial coefficients.
+        grid_params: Grid parameters for RTC.
+        apply_rtc: Whether to compute RTC gamma area.
+        gm_dem_params: DEM grid metadata for ECEF conversion.
         method: Root-finding method.
         tol: Function tolerance.
         speed: Nominal platform speed for tolerance scaling.
@@ -825,7 +844,8 @@ def backward_geocode(
         t_shift: Time shift for the secant method.
 
     Returns:
-        Orbit time and distance vector
+        A dataset containging the optimized ground_range and azimuth time
+        for each target pixel and optinal the gamma area needed for RTC.
 
     Raises:
         ValueError: If the method is not supported.
@@ -883,7 +903,7 @@ def compute_dem_area(dem_ecef: xr.DataArray, gm_dem_params: dict) -> xr.DataArra
 
     Args:
         dem_ecef: DEM in ECEF coordinates.
-        gm_dem: GridMapping of the DEM data array.
+        gm_dem_params: GridMapping metadata of the DEM data array.
 
     Returns:
         Area vectors per DEM pixel.
@@ -952,7 +972,7 @@ def compute_gamma_area(
 
     Args:
         dem_ecef: DEM in ECEF coordinates.
-        gm_dem: GridMapping of the DEM data array.
+        gm_dem_params: GridMapping metadata of the DEM data array.
         direction: Look direction vectors.
 
     Returns:
@@ -997,7 +1017,7 @@ def gamma_weights_bilinear(src_loc: xr.Dataset) -> xr.DataArray:
     """Compute bilinear gamma weights for the acquisition grid.
 
     Args:
-        acq: Acquisition dataset with indices and gamma area.
+        src_loc: Source location dataset with indices and gamma area.
 
     Returns:
         Gamma weights on the SAR grid.
@@ -1028,7 +1048,7 @@ def gamma_weights_nearest(src_loc: xr.Dataset) -> xr.DataArray:
     """Compute nearest-neighbor gamma weights for the acquisition grid.
 
     Args:
-        acq: Acquisition dataset with indices and gamma area.
+        src_loc: Source location dataset with indices and gamma area.
 
     Returns:
         Gamma weights on the SAR grid.
@@ -1047,7 +1067,7 @@ def apply_gamma_weights(
     """Apply gamma weighting block-wise.
 
     Args:
-        acq: Acquisition dataset with geometry.
+        src_loc: Source location dataset with geometry.
         func: Weighting function.
         params: Grid parameters for index conversion.
 
@@ -1055,9 +1075,9 @@ def apply_gamma_weights(
         Gamma-corrected area per pixel.
     """
     src_loc["az_idx"] = (
-        (src_loc.azimuth_time - params.az0) / _ONE_SECOND / params.d_az_scale
+        (src_loc.azimuth_time - params.az0_scale) / _ONE_SECOND / params.d_az_scale
     )
-    src_loc["gr_idx"] = (src_loc.ground_range - params.gr0) / params.d_gr_scale
+    src_loc["gr_idx"] = (src_loc.ground_range - params.gr0_scale) / params.d_gr_scale
 
     template = src_loc.gamma_area * 0
     area = xr.map_blocks(func, src_loc, template=template)
@@ -1082,7 +1102,7 @@ def fit_ground_range(time_slr_gcp: xr.DataArray, deg: int = 8) -> xr.DataArray:
 
     # polynomial fit per azimuth line
     coeff = []
-    for i, time in enumerate(x_gcp["azimuth_time"].values):
+    for i, time in enumerate(x_gcp["azimuth_time"].data):
         coeff.append(np.polyfit(x_gcp[i, :], x_gcp["ground_range"], deg=deg))
     return xr.DataArray(
         coeff,
@@ -1110,8 +1130,8 @@ def geocode_data(
 
     Args:
         data: Input dataset on the SAR grid.
-        time_az: Target azimuth times.
-        ground_range: Target ground range.
+        src_loc: Source location dataset with target coordinates.
+        grid_params: Grid parameters for index conversion.
         interp_method: Interpolation method.
 
     Returns:
@@ -1142,7 +1162,7 @@ def geocode_data(
             gr_idx.data,
             az_idx.data,
             interp_method=interp_method,
-            dtype=data.dtype,
+            dtype=data_array.dtype,
             chunks=gr_idx.data.chunks,
         )
         target_ds[var_name] = (az_idx.dims, resampled)
@@ -1155,7 +1175,7 @@ def get_source_location(
     time_slr_gcp: xr.DataArray,
     sat_position: xr.DataArray,
     grid_params: GridParams,
-    gm_dem: GridMapping,
+    gm_dem_grid: GridMapping,
     apply_rtc: bool,
 ) -> xr.Dataset:
 
@@ -1166,16 +1186,18 @@ def get_source_location(
     pos_coeff = fit_position(sat_position)
     vel_coeff = poly_derivative(pos_coeff)
 
-    data_array = xr.zeros_like(dem, dtype="float32").drop_vars("spatial_ref")
-    data_array_az = xr.zeros_like(dem, dtype="datetime64[ns]").drop_vars("spatial_ref")
+    data_array = xr.zeros_like(dem, dtype="float32")
+    azimuth_data = xr.zeros_like(dem, dtype="datetime64[ns]")
     template = xr.Dataset(
-        {"azimuth_time": data_array_az, "ground_range": data_array},
+        {"azimuth_time": azimuth_data, "ground_range": data_array},
     )
     if apply_rtc:
         template["gamma_area"] = data_array
+    if "spatial_ref" in template:
+        template = template.drop_vars("spatial_ref")
     gm_dem_params = {
-        "crs": gm_dem.crs.to_wkt(),
-        "xy_var_names": gm_dem.xy_var_names,
+        "crs": gm_dem_grid.crs.to_wkt(),
+        "xy_var_names": gm_dem_grid.xy_var_names,
     }
     out = xr.map_blocks(
         backward_geocode,
@@ -1190,7 +1212,7 @@ def get_source_location(
         },
         template=template,
     )
-    out["spatial_ref"] = xr.DataArray(0, attrs=gm_dem.crs.to_cf())
+    out.coords["spatial_ref"] = xr.DataArray(0, attrs=gm_dem_grid.crs.to_cf())
     return out
 
 
@@ -1272,30 +1294,28 @@ def _compute_indexing(
 
 def _sample_array_at_indices(
     data: np.ndarray,
-    ix: np.ndarray,
-    iy: np.ndarray,
+    x_idx: np.ndarray,
+    y_idx: np.ndarray,
     interp_method: Literal["nearest", "bilinear"] | None = None,
 ) -> np.ndarray:
-    """
-    Sample a 3d array at fractional indices (iy, ix).
-    """
+    """Sample a 3D array at fractional indices (y_idx, x_idx)."""
     if interp_method == "nearest":
-        ix_i = np.ceil(ix - 0.5).astype(np.intp)
-        iy_i = np.ceil(iy - 0.5).astype(np.intp)
-        return data[iy_i, ix_i]
-
-    ix_floor = np.floor(ix).astype(np.intp)
-    iy_floor = np.floor(iy).astype(np.intp)
-    ix_ceil = np.ceil(ix).astype(np.intp)
-    iy_ceil = np.ceil(iy).astype(np.intp)
-
-    dx = ix - ix_floor
-    dy = iy - iy_floor
-
-    v00 = data[iy_floor, ix_floor]
-    v01 = data[iy_floor, ix_ceil]
-    v10 = data[iy_ceil, ix_floor]
-    v11 = data[iy_ceil, ix_ceil]
+        x_i = np.ceil(x_idx - 0.5).astype(np.intp)
+        y_i = np.ceil(y_idx - 0.5).astype(np.intp)
+        return data[y_i, x_i]
+
+    x_floor = np.floor(x_idx).astype(np.intp)
+    y_floor = np.floor(y_idx).astype(np.intp)
+    x_ceil = np.ceil(x_idx).astype(np.intp)
+    y_ceil = np.ceil(y_idx).astype(np.intp)
+
+    dx = x_idx - x_floor
+    dy = y_idx - y_floor
+
+    v00 = data[y_floor, x_floor]
+    v01 = data[y_floor, x_ceil]
+    v10 = data[y_ceil, x_floor]
+    v11 = data[y_ceil, x_ceil]
 
     if interp_method == "bilinear":
         u0 = v00 + dx * (v01 - v00)
diff --git a/xarray_eopf/backend.py b/xarray_eopf/backend.py
index 16763b8..afb565f 100644
--- a/xarray_eopf/backend.py
+++ b/xarray_eopf/backend.py
@@ -127,9 +127,11 @@ def open_dataset(
         crs: pyproj.CRS | str | None = None,
         interp_methods: SpatialInterpMethods | Sen1InterpMethods | None = None,
         agg_methods: SpatialAggMethods | None = None,
+        # params for Sentinel-1 specifically
         dem: xr.DataArray | None = None,
         footprint_scale_factor: tuple[float | int, float | int] | None = None,
         apply_rtc: bool = True,
+        cache_uri: str | None = None,
         # params required by xarray backend interface
         drop_variables: str | Iterable[str] | None = None,
         # params for other reasons
@@ -207,6 +209,10 @@ def open_dataset(
                 when not provided.
             apply_rtc: Whether to apply radiometric terrain correction (RTC) for
                 Sentinel-1 analysis mode. Defaults to `True`.
+            cache_uri: Temporary path used to store intermediate results from the
+                backward geocoding step in the Sentinel-1 processing workflow. The cache
+                is automatically removed when the Python process exits. If None, a
+                temporary directory with a unique UUID-based name is created.
             variables: Variables to include in the dataset. Can be a name or
                 regex pattern or iterable of the latter.
             drop_variables: Variable name or iterable of variable names

From 11bb3b8854d4c6bc822b0525ef8da6b418991a50 Mon Sep 17 00:00:00 2001
From: konstntokas <konstantin.ntokas@brockmann-consult.de>
Date: Wed, 8 Jul 2026 16:56:30 +0200
Subject: [PATCH 5/9] test coverage 100

---
 tests/amodes/test_sentinel1.py | 327 ++++++++++++++++++++++++++++-----
 1 file changed, 280 insertions(+), 47 deletions(-)

diff --git a/tests/amodes/test_sentinel1.py b/tests/amodes/test_sentinel1.py
index 3c88005..9677875 100644
--- a/tests/amodes/test_sentinel1.py
+++ b/tests/amodes/test_sentinel1.py
@@ -3,6 +3,7 @@
 #  https://opensource.org/license/apache-2-0.
 
 import os
+import uuid
 from types import SimpleNamespace
 from unittest import TestCase
 from unittest.mock import patch
@@ -17,7 +18,7 @@
 from tests.helpers import make_s1_grd_datatree, make_s1_ocn_datatree
 from xarray_eopf.amode import AnalysisModeRegistry
 from xarray_eopf.amodes import sentinel1 as sen1
-from xarray_eopf.amodes.sentinel1 import Sen1GRD, Sen1OCN, register, GridParams
+from xarray_eopf.amodes.sentinel1 import Sen1GRD, Sen1OCN, register
 
 
 class Sentinel1AnalysisModeTest(TestCase):
@@ -53,9 +54,37 @@ def test_transform_dataset(self: TestCase):
 
 
 class Sen1GRDTest(Sen1TestMixin, TestCase):
-    mode = Sen1GRD()
-    dem = xr.DataArray(np.ones((2, 2), dtype="float32"), dims=("lat", "lon"))
-    dt = make_s1_grd_datatree()
+
+    def setUp(self):
+        self.mode = Sen1GRD()
+        self.dem = xr.DataArray(
+            np.ones((2, 2), dtype="float32"),
+            dims=("lat", "lon"),
+            coords={"lat": [0.0, 1.0], "lon": [0.0, 1.0]},
+        )
+        self.dt = make_s1_grd_datatree()
+        self.expected_vv = xr.Dataset(
+            {"vv": xr.DataArray(np.ones((2, 2)), dims=("lat", "lon"))}
+        )
+        self.expected_beta0_vv = xr.Dataset(
+            {
+                "beta0_vv": xr.DataArray(
+                    np.ones((2, 2)),
+                    dims=("lat", "lon"),
+                    coords={"lat": [0.0, 1.0], "lon": [0.0, 1.0]},
+                )
+            },
+        )
+        self.src_loc = xr.Dataset(
+            {
+                "azimuth_time": (
+                    ("lat", "lon"),
+                    np.zeros((2, 2), dtype="datetime64[ns]"),
+                ),
+                "ground_range": (("lat", "lon"), np.zeros((2, 2))),
+                "gamma_area": (("lat", "lon"), np.ones((2, 2))),
+            }
+        )
 
     def test_is_valid_source_ok(self):
         self.assertTrue(self.mode.is_valid_source("data/S1A_IW_GRDH_20240201.zarr"))
@@ -80,32 +109,38 @@ def test_get_grid_parameters(self):
         )
 
     def test_get_applicable_params(self: TestCase):
-        dem = xr.DataArray(np.ones((2, 2)), dims=("lat", "lon"))
+
         self.assertEqual({}, self.mode.get_applicable_params())
         self.assertEqual(
             {
                 "resolution": 10,
                 "bbox": [1, 3, 4, 5],
                 "crs": pyproj.CRS.from_string("EPSG:4326"),
-                "dem": dem,
+                "dem": self.dem,
                 "interp_methods": "nearest",
                 "footprint_scale_factor": (2.0, 3.0),
                 "apply_rtc": False,
+                "cache_uri": "file:///tmp/cache",
             },
             self.mode.get_applicable_params(
                 resolution=10,
                 bbox=[1, 3, 4, 5],
                 crs="EPSG:4326",
-                dem=dem,
+                dem=self.dem,
                 interp_methods="nearest",
                 footprint_scale_factor=(2.0, 3.0),
                 apply_rtc=False,
+                cache_uri="file:///tmp/cache",
             ),
         )
         with pytest.raises(TypeError, match="interp_methods"):
             self.mode.get_applicable_params(interp_methods="cubic")
         with pytest.raises(TypeError, match="footprint_scale_factor"):
             self.mode.get_applicable_params(footprint_scale_factor=(1.0, "x"))
+        with pytest.raises(TypeError, match="apply_rtc"):
+            self.mode.get_applicable_params(apply_rtc="yes")
+        with pytest.raises(TypeError, match="cache_uri"):
+            self.mode.get_applicable_params(cache_uri=123)
 
     def test_convert_datatree(self):
         expected = xr.Dataset(
@@ -125,13 +160,29 @@ def test_convert_datatree(self):
         self.assertTrue(kwargs["apply_rtc"])
         self.assertIn("gr0", args[4])
 
-    def test_convert_datatree_uses_get_dem(self):
+    def test_convert_datatree_with_cache_uri_uses_fs(self):
         expected = xr.Dataset(
             {"vv": xr.DataArray(np.ones((2, 2)), dims=("lat", "lon"))}
         )
+        fs = SimpleNamespace()
+        with (
+            patch.object(
+                sen1.fsspec, "url_to_fs", return_value=(fs, "/cache")
+            ) as url_to_fs,
+            patch.object(self.mode, "_terrain_correct", return_value=expected),
+        ):
+            _ = self.mode.convert_datatree(
+                self.dt, includes=["vv"], dem=self.dem, cache_uri="file:///cache/"
+            )
+        url_to_fs.assert_called_once_with("file:///cache")
+        self.assertEqual("file:///cache", self.mode.cache_uri)
+
+    def test_convert_datatree_uses_get_dem(self):
 
         with patch.object(sen1, "get_dem", return_value=self.dem) as get_dem_mock:
-            with patch.object(self.mode, "_terrain_correct", return_value=expected):
+            with patch.object(
+                self.mode, "_terrain_correct", return_value=self.expected_beta0_vv
+            ):
                 _ = self.mode.convert_datatree(self.dt, includes=["vv"])
 
         get_dem_mock.assert_called_once()
@@ -142,6 +193,79 @@ def test_convert_datatree_fail(self):
         with pytest.raises(ValueError, match="No valid variable names"):
             self.mode.convert_datatree(self.dt, includes="bibo", dem=self.dem)
 
+    def test_convert_datatree_cleans_up_on_failure(self):
+        with (
+            patch.object(
+                self.mode, "_terrain_correct", side_effect=RuntimeError("boom")
+            ),
+            patch.object(self.mode, "_cleanup") as cleanup,
+        ):
+            with pytest.raises(RuntimeError, match="boom"):
+                self.mode.convert_datatree(self.dt, includes=["vv"], dem=self.dem)
+        cleanup.assert_called_once()
+
+    def test_terrain_correct_with_rtc_nearest(self):
+        with (
+            patch.object(sen1, "get_source_location", return_value=self.src_loc),
+            patch.object(sen1, "geocode_data", return_value=self.expected_beta0_vv),
+            patch.object(
+                sen1,
+                "apply_gamma_weights",
+                return_value=xr.ones_like(self.src_loc.gamma_area),
+            ) as gamma_mock,
+            patch.object(sen1, "assign_grid_mapping", side_effect=lambda ds: ds),
+        ):
+            self.mode.cache_uri = f"tmp_{uuid.uuid4().hex}"
+            out = self.mode._terrain_correct(
+                self.expected_beta0_vv,
+                xr.DataArray(np.zeros((2, 2)), dims=("lat", "lon")),
+                xr.DataArray(np.zeros((2, 2, 3)), dims=("lat", "lon", "axis")),
+                self.dem,
+                self.mode._get_grid_parameters(self.dt, (3.0, 3.0)),
+                apply_rtc=True,
+                interp_method="nearest",
+            )
+        gamma_mock.assert_called_once()
+        self.assertIn("gamma0_vv", out.data_vars)
+
+    def test_terrain_correct_with_rtc_bilinear(self):
+        with (
+            patch.object(sen1, "get_source_location", return_value=self.src_loc),
+            patch.object(sen1, "geocode_data", return_value=self.expected_beta0_vv),
+            patch.object(
+                sen1,
+                "apply_gamma_weights",
+                return_value=xr.ones_like(self.src_loc.gamma_area),
+            ) as gamma_mock,
+            patch.object(sen1, "assign_grid_mapping", side_effect=lambda ds: ds),
+        ):
+            self.mode.cache_uri = f"tmp_{uuid.uuid4().hex}"
+            out = self.mode._terrain_correct(
+                self.expected_beta0_vv,
+                xr.DataArray(np.zeros((2, 2)), dims=("lat", "lon")),
+                xr.DataArray(np.zeros((2, 2, 3)), dims=("lat", "lon", "axis")),
+                self.dem,
+                self.mode._get_grid_parameters(self.dt, (3.0, 3.0)),
+                apply_rtc=True,
+                interp_method="bilinear",
+            )
+        gamma_mock.assert_called_once()
+        self.assertIn("gamma0_vv", out.data_vars)
+
+    def test_cleanup_removes_cache(self):
+        with patch.object(sen1.fsspec, "url_to_fs") as url_to_fs:
+            self.mode.cache_uri = "file:///tmp/fake-cache"
+            fs = SimpleNamespace(
+                exists=lambda path: True, rm=lambda path, recursive: None
+            )
+            url_to_fs.return_value = (fs, "/tmp/fake-cache")
+            self.mode._cleanup()
+            url_to_fs.assert_called_once_with("file:///tmp/fake-cache")
+
+    def test_cleanup_without_cache_uri_is_noop(self):
+        self.mode.cache_uri = None
+        self.mode._cleanup()
+
 
 class Sen1OCNTest(Sen1TestMixin, TestCase):
     mode = Sen1OCN()
@@ -254,7 +378,13 @@ def setUp(self):
         self.dem = xr.DataArray(
             np.ones((2, 2), dtype="float64"),
             dims=("lat", "lon"),
-            coords={"lat": [0.0, 0.1], "lon": [0.0, 0.1]},
+            coords={
+                "lat": [0.0, 0.1],
+                "lon": [0.0, 0.1],
+                "spatial_ref": xr.DataArray(
+                    0, attrs=pyproj.CRS.from_epsg(4326).to_cf()
+                ),
+            },
         )
         self.dem_ecef = xr.DataArray(
             np.ones((3, 2, 2), dtype="float64"),
@@ -279,6 +409,16 @@ def setUp(self):
             },
             attrs=dict(mean=1, std=1),
         )
+        self.time_slr = xr.DataArray(
+            np.ones((2, 2), dtype="float64"),
+            dims=("azimuth_time", "ground_range"),
+            coords={"azimuth_time": [0, 1], "ground_range": [0, 1]},
+        )
+        self.sat_position = xr.DataArray(
+            np.ones((2, 3), dtype="float64"),
+            dims=("azimuth_time", "axis"),
+            coords={"azimuth_time": [0, 1], "axis": ["x", "y", "z"]},
+        )
 
     def test_gridparams_iter(self):
         self.assertEqual(
@@ -456,27 +596,127 @@ def test_fit_position_and_poly_derivative(self):
         deriv = sen1.poly_derivative(coeff)
         self.assertEqual(coeff.sizes["degree"] - 1, deriv.sizes["degree"])
 
-    def test_zero_doppler_and_prime(self):
-        dem_ecef = xr.DataArray(
-            np.ones((3, 2, 2), dtype="float64"),
-            dims=("axis", "lat", "lon"),
-            coords={"axis": ["x", "y", "z"], "lat": [0.0, 1.0], "lon": [0.0, 1.0]},
+    def test_get_source_location_and_assign_grid_mapping(self):
+        dem = xr.DataArray(
+            np.ones((2, 2), dtype="float64"),
+            dims=("lat", "lon"),
+            coords={"lat": [0.0, 1.0], "lon": [0.0, 1.0]},
         )
-        time_orbit = xr.DataArray(np.zeros((2, 2)), dims=("lat", "lon"))
-        pos_coeff = xr.DataArray(
-            np.zeros((2, 3)),
-            dims=("degree", "axis"),
-            coords={"degree": [1, 0], "axis": ["x", "y", "z"]},
+        time_slr = xr.DataArray(
+            np.ones((2, 2), dtype="float64"),
+            dims=("azimuth_time", "ground_range"),
+            coords={"azimuth_time": [0, 1], "ground_range": [0, 1]},
         )
-        vel_coeff = xr.DataArray(
-            np.zeros((2, 3)),
-            dims=("degree", "axis"),
-            coords={"degree": [1, 0], "axis": ["x", "y", "z"]},
+        sat_position = xr.DataArray(
+            np.ones((2, 3), dtype="float64"),
+            dims=("azimuth_time", "axis"),
+            coords={"azimuth_time": [0, 1], "axis": ["x", "y", "z"]},
+        )
+        with (
+            patch.object(sen1, "fit_ground_range", return_value=self.gr_coeff),
+            patch.object(sen1, "fit_position", return_value=self.posvel_coeff),
+            patch.object(sen1, "backward_geocode") as bg,
+        ):
+            bg.return_value = xr.Dataset(
+                {
+                    "azimuth_time": xr.DataArray(
+                        np.zeros((2, 2), dtype="datetime64[ns]"), dims=("lat", "lon")
+                    ),
+                    "ground_range": xr.DataArray(np.zeros((2, 2)), dims=("lat", "lon")),
+                    "gamma_area": xr.DataArray(np.ones((2, 2)), dims=("lat", "lon")),
+                }
+            )
+            gm_dem = GridMapping.from_dataset(dem.to_dataset(name="dem"))
+            out = sen1.get_source_location(
+                dem,
+                time_slr,
+                sat_position,
+                self.grid_params,
+                gm_dem,
+                True,
+            )
+        self.assertIn("spatial_ref", out.coords)
+        out = sen1.assign_grid_mapping(
+            xr.Dataset({"a": xr.DataArray([1], dims=("x",))})
+        )
+        self.assertEqual("spatial_ref", out["a"].attrs["grid_mapping"])
+
+    def test_get_source_location_without_rtc(self):
+
+        gm_dem = GridMapping.from_dataset(self.dem.to_dataset(name="dem"))
+        with (
+            patch.object(
+                sen1,
+                "backward_geocode",
+                return_value=xr.Dataset(
+                    {
+                        "azimuth_time": xr.DataArray(
+                            np.zeros((2, 2), dtype="datetime64[ns]"),
+                            dims=("lat", "lon"),
+                        ),
+                        "ground_range": xr.DataArray(
+                            np.zeros((2, 2)), dims=("lat", "lon")
+                        ),
+                    }
+                ),
+            ) as bg,
+            patch.object(sen1, "fit_ground_range", return_value=self.gr_coeff),
+            patch.object(sen1, "fit_position", return_value=self.posvel_coeff),
+        ):
+            out = sen1.get_source_location(
+                self.dem,
+                self.time_slr,
+                self.sat_position,
+                self.grid_params,
+                gm_dem,
+                False,
+            )
+        bg.assert_called_once()
+        self.assertNotIn("gamma_area", out.data_vars)
+
+    def test_compute_indexing_and_sample_array_errors(self):
+        data = xr.Dataset(
+            {
+                "a": xr.DataArray(
+                    np.arange(9).reshape(3, 3), dims=("azimuth_time", "ground_range")
+                )
+            }
+        )
+        az_idx = xr.DataArray(
+            da.from_array(np.array([[0.2, 1.2], [0.2, 1.2]]), chunks=(2, 2)),
+            dims=("lat", "lon"),
+        )
+        gr_idx = xr.DataArray(
+            da.from_array(np.array([[0.2, 1.2], [0.2, 1.2]]), chunks=(2, 2)),
+            dims=("lat", "lon"),
+        )
+        indexing = sen1._compute_indexing(data, az_idx, gr_idx)
+        np.testing.assert_array_equal(
+            indexing.ij_bboxes,
+            np.array([[[0]], [[0]], [[3]], [[3]]], dtype=np.int32),
+        )
+
+        arr = np.array([[1.0, 2.0], [3.0, 4.0]])
+        nearest = sen1._sample_array_at_indices(
+            arr, np.array([[0.6]]), np.array([[1.4]]), "nearest"
+        )
+        np.testing.assert_array_equal(nearest, np.array([[4.0]]))
+        bilinear = sen1._sample_array_at_indices(
+            arr, np.array([[0.5]]), np.array([[0.5]]), "bilinear"
+        )
+        np.testing.assert_allclose(bilinear, np.array([[2.5]]))
+        with pytest.raises(NotImplementedError, match="interp_methods"):
+            sen1._sample_array_at_indices(
+                np.zeros((2, 2)), np.zeros((2, 2)), np.zeros((2, 2)), "cubic"
+            )
+
+    def test_zero_doppler_and_prime(self):
+        time_orbit = xr.DataArray(np.zeros((2, 2)), dims=("lat", "lon"))
+        f, payload = sen1.zero_doppler(
+            self.dem_ecef, self.posvel_coeff, self.posvel_coeff, time_orbit
         )
-        vel_coeff.loc[{"degree": 0, "axis": "x"}] = 1.0
-        f, payload = sen1.zero_doppler(dem_ecef, pos_coeff, vel_coeff, time_orbit)
         self.assertEqual(("lat", "lon"), f.dims)
-        fp = sen1.zero_doppler_prime(vel_coeff, time_orbit, payload)
+        fp = sen1.zero_doppler_prime(self.posvel_coeff, time_orbit, payload)
         self.assertEqual(("lat", "lon"), fp.dims)
 
     def test_secant_and_newton(self):
@@ -580,17 +820,6 @@ def test_backward_geocode_secant_and_newton_paths(self):
         self.assertEqual(3, len(out_newton))
 
     def test_compute_gamma_area_clips_negative(self):
-        dem_ecef = xr.DataArray(
-            np.ones((3, 2, 2), dtype="float64"),
-            dims=("axis", "lat", "lon"),
-            coords={"axis": ["x", "y", "z"], "lat": [0.0, 1.0], "lon": [0.0, 1.0]},
-        )
-        gm_dem = {
-            "crs": GridMapping.from_dataset(
-                dem_ecef.to_dataset(name="dem")
-            ).crs.to_wkt(),
-            "xy_var_names": ("lon", "lat"),
-        }
         area = xr.DataArray(
             np.array(
                 [
@@ -600,7 +829,7 @@ def test_compute_gamma_area_clips_negative(self):
                 ]
             ),
             dims=("axis", "lat", "lon"),
-            coords=dem_ecef.coords,
+            coords=self.dem_ecef.coords,
         )
         direction = xr.DataArray(
             np.array(
@@ -611,10 +840,12 @@ def test_compute_gamma_area_clips_negative(self):
                 ]
             ),
             dims=("axis", "lat", "lon"),
-            coords=dem_ecef.coords,
+            coords=self.dem_ecef.coords,
         )
         with patch.object(sen1, "compute_dem_area", return_value=area):
-            gamma = sen1.compute_gamma_area(dem_ecef, gm_dem, direction)
+            gamma = sen1.compute_gamma_area(
+                self.dem_ecef, self.gm_dem_params, direction
+            )
         self.assertTrue(np.all(gamma.values >= 0))
         self.assertEqual(0.0, float(gamma.values[0, 1]))
 
@@ -634,7 +865,7 @@ def test_compute_dem_area(self):
         self.assertEqual(("axis", "lat", "lon"), area.dims)
 
     def test_sum_weights_and_gamma_weight_helpers(self):
-        acq = xr.Dataset(
+        scr_indices = xr.Dataset(
             {
                 "gamma_area": xr.DataArray([[1.0, 2.0]], dims=("lat", "lon")),
                 "az_idx": xr.DataArray([[0.2, 0.8]], dims=("lat", "lon")),
@@ -647,19 +878,21 @@ def test_sum_weights_and_gamma_weight_helpers(self):
             coords={"gr_idx": [1], "az_idx": [1]},
         )
         with patch.object(sen1.flox.xarray, "xarray_reduce", return_value=reduced):
-            summed = sen1.sum_weights(acq.gamma_area, acq.az_idx, acq.gr_idx)
+            summed = sen1.sum_weights(
+                scr_indices.gamma_area, scr_indices.az_idx, scr_indices.gr_idx
+            )
         self.assertEqual(("lat", "lon"), summed.dims)
         self.assertEqual((1, 2), summed.data.shape)
 
         with patch.object(
-            sen1, "sum_weights", return_value=xr.zeros_like(acq.gamma_area)
+            sen1, "sum_weights", return_value=xr.zeros_like(scr_indices.gamma_area)
         ) as sw:
-            _ = sen1.gamma_weights_bilinear(acq)
+            _ = sen1.gamma_weights_bilinear(scr_indices)
             self.assertEqual(4, sw.call_count)
         with patch.object(
-            sen1, "sum_weights", return_value=xr.zeros_like(acq.gamma_area)
+            sen1, "sum_weights", return_value=xr.zeros_like(scr_indices.gamma_area)
         ) as sw:
-            _ = sen1.gamma_weights_nearest(acq)
+            _ = sen1.gamma_weights_nearest(scr_indices)
             sw.assert_called_once()
 
     def test_apply_gamma_weights(self):

From 5ad5c7bf727d46bc5eaae1fa6a6055ec1a4ddeee Mon Sep 17 00:00:00 2001
From: konstntokas <konstantin.ntokas@brockmann-consult.de>
Date: Wed, 8 Jul 2026 17:05:16 +0200
Subject: [PATCH 6/9] added cleanup after test

---
 tests/amodes/test_sentinel1.py | 20 ++++++++++----------
 1 file changed, 10 insertions(+), 10 deletions(-)

diff --git a/tests/amodes/test_sentinel1.py b/tests/amodes/test_sentinel1.py
index 9677875..d71de41 100644
--- a/tests/amodes/test_sentinel1.py
+++ b/tests/amodes/test_sentinel1.py
@@ -143,39 +143,34 @@ def test_get_applicable_params(self: TestCase):
             self.mode.get_applicable_params(cache_uri=123)
 
     def test_convert_datatree(self):
-        expected = xr.Dataset(
-            {"vv": xr.DataArray(np.ones((2, 2)), dims=("lat", "lon"))}
-        )
-
         with patch.object(
-            self.mode, "_terrain_correct", return_value=expected
+            self.mode, "_terrain_correct", return_value=self.expected_vv
         ) as mocked:
             out = self.mode.convert_datatree(self.dt, includes=["vv"], dem=self.dem)
 
-        self.assertIs(out, expected)
+        self.assertIs(out, self.expected_vv)
         args, kwargs = mocked.call_args
         self.assertEqual(["beta0_vv"], list(args[0].data_vars))
         self.assertIs(args[3], self.dem)
         self.assertEqual("bilinear", kwargs["interp_method"])
         self.assertTrue(kwargs["apply_rtc"])
         self.assertIn("gr0", args[4])
+        self.mode._cleanup()
 
     def test_convert_datatree_with_cache_uri_uses_fs(self):
-        expected = xr.Dataset(
-            {"vv": xr.DataArray(np.ones((2, 2)), dims=("lat", "lon"))}
-        )
         fs = SimpleNamespace()
         with (
             patch.object(
                 sen1.fsspec, "url_to_fs", return_value=(fs, "/cache")
             ) as url_to_fs,
-            patch.object(self.mode, "_terrain_correct", return_value=expected),
+            patch.object(self.mode, "_terrain_correct", return_value=self.expected_vv),
         ):
             _ = self.mode.convert_datatree(
                 self.dt, includes=["vv"], dem=self.dem, cache_uri="file:///cache/"
             )
         url_to_fs.assert_called_once_with("file:///cache")
         self.assertEqual("file:///cache", self.mode.cache_uri)
+        self.mode._cleanup()
 
     def test_convert_datatree_uses_get_dem(self):
 
@@ -188,10 +183,12 @@ def test_convert_datatree_uses_get_dem(self):
         get_dem_mock.assert_called_once()
         args, _ = get_dem_mock.call_args
         self.assertEqual(self.dt.attrs["stac_discovery"]["bbox"], args[0])
+        self.mode._cleanup()
 
     def test_convert_datatree_fail(self):
         with pytest.raises(ValueError, match="No valid variable names"):
             self.mode.convert_datatree(self.dt, includes="bibo", dem=self.dem)
+        self.mode._cleanup()
 
     def test_convert_datatree_cleans_up_on_failure(self):
         with (
@@ -203,6 +200,7 @@ def test_convert_datatree_cleans_up_on_failure(self):
             with pytest.raises(RuntimeError, match="boom"):
                 self.mode.convert_datatree(self.dt, includes=["vv"], dem=self.dem)
         cleanup.assert_called_once()
+        self.mode._cleanup()
 
     def test_terrain_correct_with_rtc_nearest(self):
         with (
@@ -227,6 +225,7 @@ def test_terrain_correct_with_rtc_nearest(self):
             )
         gamma_mock.assert_called_once()
         self.assertIn("gamma0_vv", out.data_vars)
+        self.mode._cleanup()
 
     def test_terrain_correct_with_rtc_bilinear(self):
         with (
@@ -251,6 +250,7 @@ def test_terrain_correct_with_rtc_bilinear(self):
             )
         gamma_mock.assert_called_once()
         self.assertIn("gamma0_vv", out.data_vars)
+        self.mode._cleanup()
 
     def test_cleanup_removes_cache(self):
         with patch.object(sen1.fsspec, "url_to_fs") as url_to_fs:

From 198752aca46c232949d6f0e889ee220537cb3f5a Mon Sep 17 00:00:00 2001
From: Konstantin Ntokas <38956538+konstntokas@users.noreply.github.com>
Date: Fri, 10 Jul 2026 08:12:24 +0200
Subject: [PATCH 7/9] Update xarray_eopf/amodes/sentinel1.py

Co-authored-by: Pontus Lurcock <pont@talvi.net>
---
 xarray_eopf/amodes/sentinel1.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/xarray_eopf/amodes/sentinel1.py b/xarray_eopf/amodes/sentinel1.py
index 5b9cc5f..1bc1423 100644
--- a/xarray_eopf/amodes/sentinel1.py
+++ b/xarray_eopf/amodes/sentinel1.py
@@ -844,7 +844,7 @@ def backward_geocode(
         t_shift: Time shift for the secant method.
 
     Returns:
-        A dataset containging the optimized ground_range and azimuth time
+        A dataset containing the optimized ground_range and azimuth time
         for each target pixel and optinal the gamma area needed for RTC.
 
     Raises:

From 2078d72a9db3942cb3c785407e236dc2bac2d3de Mon Sep 17 00:00:00 2001
From: Konstantin Ntokas <38956538+konstntokas@users.noreply.github.com>
Date: Fri, 10 Jul 2026 08:12:40 +0200
Subject: [PATCH 8/9] Update xarray_eopf/amodes/sentinel1.py

Co-authored-by: Pontus Lurcock <pont@talvi.net>
---
 xarray_eopf/amodes/sentinel1.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/xarray_eopf/amodes/sentinel1.py b/xarray_eopf/amodes/sentinel1.py
index 1bc1423..b4e0f3f 100644
--- a/xarray_eopf/amodes/sentinel1.py
+++ b/xarray_eopf/amodes/sentinel1.py
@@ -845,7 +845,7 @@ def backward_geocode(
 
     Returns:
         A dataset containing the optimized ground_range and azimuth time
-        for each target pixel and optinal the gamma area needed for RTC.
+        for each target pixel and optionally the gamma area needed for RTC.
 
     Raises:
         ValueError: If the method is not supported.

From 8033e96cd56438c2487bfbd26fbe61cbdf3e9f79 Mon Sep 17 00:00:00 2001
From: konstntokas <konstantin.ntokas@brockmann-consult.de>
Date: Fri, 10 Jul 2026 08:26:40 +0200
Subject: [PATCH 9/9] address comments

---
 xarray_eopf/backend.py | 9 +++++----
 1 file changed, 5 insertions(+), 4 deletions(-)

diff --git a/xarray_eopf/backend.py b/xarray_eopf/backend.py
index afb565f..1526465 100644
--- a/xarray_eopf/backend.py
+++ b/xarray_eopf/backend.py
@@ -209,10 +209,11 @@ def open_dataset(
                 when not provided.
             apply_rtc: Whether to apply radiometric terrain correction (RTC) for
                 Sentinel-1 analysis mode. Defaults to `True`.
-            cache_uri: Temporary path used to store intermediate results from the
-                backward geocoding step in the Sentinel-1 processing workflow. The cache
-                is automatically removed when the Python process exits. If None, a
-                temporary directory with a unique UUID-based name is created.
+            cache_uri: Temporary path, interpreted as a fsspec `urlpath`, where
+                intermediate results from the backward geocoding step of the Sentinel-1
+                processing workflow are stored. The cache is automatically deleted
+                when the Python process exits. If `None`, a temporary directory with
+                a unique UUID-based name is created automatically.
             variables: Variables to include in the dataset. Can be a name or
                 regex pattern or iterable of the latter.
             drop_variables: Variable name or iterable of variable names