[Git][debian-gis-team/xarray-eopf][upstream] New upstream version 0.3.0
Antonio Valentino (@antonio.valentino)
gitlab at salsa.debian.org
Wed Jul 29 22:12:51 BST 2026
Antonio Valentino pushed to branch upstream at Debian GIS Project / xarray-eopf
Commits:
df45a2b4 by Antonio Valentino at 2026-07-29T20:55:20+00:00
New upstream version 0.3.0
- - - - -
22 changed files:
- CHANGES.md
- docs/examples/introduction.ipynb
- docs/examples/sentinel_1_analysis.ipynb
- docs/examples/sentinel_1_native.ipynb
- docs/guide.md
- docs/index.md
- examples/introduction.ipynb
- examples/sentinel_1_analysis.ipynb
- examples/sentinel_1_native.ipynb
- integration/test_data/dem_small.zarr.zip → integration/test_data/dem_grd.zarr.zip
- + integration/test_data/dem_slc.zarr.zip
- integration/test_sen1_analysis.py
- integration/test_sen1_native.py
- mkdocs.yml
- tests/amodes/test_sentinel1.py
- tests/amodes/test_sentinel3.py
- tests/helpers/__init__.py
- tests/helpers/sentinel1.py
- xarray_eopf/amodes/sentinel1.py
- xarray_eopf/amodes/sentinel3.py
- xarray_eopf/backend.py
- xarray_eopf/version.py
Changes:
=====================================
CHANGES.md
=====================================
@@ -1,9 +1,17 @@
+## Changes in 0.3.0 (from 2026-07-22)
+
+- Sentinel-1 GRD analysis mode is now fully lazy, enabling seamless execution on
+ local and distributed Dask clusters.
+- Add support for Sentinel-1 Level-1 SLC analysis mode.
+
+
## 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`
=====================================
docs/examples/introduction.ipynb
=====================================
@@ -4,9 +4,7 @@
"cell_type": "markdown",
"id": "c04818c2-3b78-4522-8f7d-aedda2a57c5d",
"metadata": {},
- "source": [
- "# Intoduction to the xarray EOPF backend"
- ]
+ "source": "# Introduction to the xarray EOPF backend"
},
{
"cell_type": "markdown",
@@ -37,7 +35,7 @@
"source": [
"## Install the xarray-eopf Backend\n",
"\n",
- "The backend is implemented as an **xarray plugin** and can be installed using either `pip` or `conda/mamba` from the conda-forge channel.\n",
+ "The backend is implemented as a **xarray plugin** and can be installed using either `pip` or `conda/mamba` from the conda-forge channel.\n",
"\n",
"- 📦 **PyPI:** [xarray-eopf on PyPI](https://pypi.org/project/xarray-eopf/) `pip install xarray-eopf`\n",
"- 🐍 **Conda (conda-forge):** [xarray-eopf on Anaconda](https://anaconda.org/conda-forge/xarray-eopf) `conda install -c conda-forge xarray-eopf`\n",
@@ -89,7 +87,7 @@
"source": [
"## Main Features of the xarray-eopf Backend\n",
"\n",
- "The xarray-eopf backend for EOPF data products can be selecterd by setting `engine=\"eopf-zarr\"` in `xarray.open_dataset(..)` and `xarray.open_datatree(..)` method. All data access is *lazy*, meaning that data is only loaded when required—for example, during plotting or when writing to storage. It supports two modes of operation:\n",
+ "The xarray-eopf backend for EOPF data products can be selected by setting `engine=\"eopf-zarr\"` in `xarray.open_dataset(..)` and `xarray.open_datatree(..)` method. All data access is *lazy*, meaning that data is only loaded when required—for example, during plotting or when writing to storage. It supports two modes of operation:\n",
"\n",
"- **Analysis Mode** *(default)*\n",
"- **Native Mode**\n",
=====================================
docs/examples/sentinel_1_analysis.ipynb
=====================================
The diff for this file was not included because it is too large.
=====================================
docs/examples/sentinel_1_native.ipynb
=====================================
The diff for this file was not included because it is too large.
=====================================
docs/guide.md
=====================================
@@ -3,8 +3,8 @@ namely _analysis mode_ (the default) and _native mode_, which are described in
the following.
An introductory example notebook is available at :
-- [Docs - Intoduction to the xarray EOPF backend](https://eopf-sample-service.github.io/xarray-eopf/examples/introduction/)
-- [Notebook Gallery - Intoduction to the xarray EOPF backend](https://eopf-sample-service.github.io/eopf-sample-notebooks/introduction/)
+- [Docs - Introduction to the xarray EOPF backend](https://eopf-sample-service.github.io/xarray-eopf/examples/introduction/)
+- [Notebook Gallery - Introduction to the xarray EOPF backend](https://eopf-sample-service.github.io/eopf-sample-notebooks/introduction/)
---
@@ -58,12 +58,12 @@ Additional parameters specific to each Sentinel mission are described below.
Processing workflows differ significantly across Sentinel-1 product types. Therefore,
each product family is documented in its own dedicated section.
-> **Note:** Support for SLC products is planned for a future release.
-
##### Sentinel-1 Level-1 GRD
> Note: Support for Sentinel-1 GRD products in analysis mode is
-> currently experimental and undergoing validation.
+> currently experimental and undergoing validation. Some conversion parameters
+> are missing in the new EOPF product, which are currently estimated. Newer EOPF
+> product version will include these parameters.
Sentinel-1 Level-1 GRD data is provided in radar geometry, defined by the coordinates
(`azimuth_time`, `ground_range`). To transform this data into an
@@ -110,12 +110,80 @@ 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:
- [Docs – Sentinel-1 Analysis Mode](https://eopf-sample-service.github.io/xarray-eopf/examples/sentinel_1_analysis/)
+##### Sentinel-1 Level-1 SLC
+
+> Note: Support for Sentinel-1 SLC products in analysis mode is
+> currently experimental and undergoing validation.
+
+Sentinel-1 Level-1 SLC data is provided in radar geometry, defined by the coordinates
+(`azimuth_time`, `slant_range_time`) and organized in bursts and swaths. To transform
+this data into an **analysis-ready dataset**, the following processing steps are applied:
+
+1. **Radiometric Calibration:** For each burst complex SLC measurements are
+ converted into `beta0` backscatter values using the `beta_nought` calibration
+ lookup table (LUT).
+2. **Burst and Swath Merging:** Valid burst regions are extracted using burst metadata,
+ merged along azimuth time, aligned across swaths, and then merged along slant range
+ to produce one continuous acquisition grid per selected polarization.
+3. **Geometric Terrain Correction (GTC):** Using a Digital Elevation Model (DEM),
+ the processor performs **inverse geocoding** by solving the zero-Doppler equation
+ based on satellite orbit information and terrain elevation. This step maps the
+ data from radar geometry to a georeferenced grid.
+4. **Radiometric Terrain Correction (RTC):** (Optional) RTC compensates for
+ terrain-induced radiometric distortions such as foreshortening, layover,
+ and slope-dependent brightness variations.
+
+📖 [D. Small, *Flattening Gamma: Radiometric Terrain Correction for SAR Imagery*](https://ieeexplore.ieee.org/document/5752845)
+
+**Supported Products:**
+
+- Sentinel-1 Level-1 SLC
+
+**Supported Variables**
+
+- **Polarization bands**:
+ `vv`, `vh`, `hh`, `hv` *(each SLC product contains only a subset of these bands)*
+
+**Specific Sentinel-1 Level-1 SLC parameters `**kwargs`:**
+
+- `crs`: Coordinate reference system of the output dataset. Can be provided as a
+ `str` or a `pyproj.CRS` object. If a string is given, it will be parsed using
+ [`pyproj.crs.CRS.from_string`](https://pyproj4.github.io/pyproj/dev/api/crs/crs.html#pyproj.crs.CRS.from_string).
+ If not specified, [EPSG:4326](https://epsg.io/4326) is used.
+- `resolution`: Target resolution for all spatial variables expressed in the units
+ of the specified `crs`. If not specified, the resolution is derived from the DEM
+ (see `dem` below).
+- `dem`: Digital Elevation Model (DEM) as a CF-compliant `xarray.DataArray` used for
+ terrain correction. If provided, the parameters `crs`, `bbox`, and `resolution` are
+ ignored, and the target grid is derived from the DEM. If not provided, the
+ [CopDEM COG (30 m)](https://browser.stac.dataspace.copernicus.eu/collections/cop-dem-glo-30-dged-cog) in geographic coordinates
+ is retrieved automatically via the CDSE STAC API. This requires
+ [CDSE S3 credentials](https://documentation.dataspace.copernicus.eu/APIs/S3.html#generate-secrets).
+- `apply_rtc`: Enable or disable radiometric terrain correction (RTC). Default is `True`.
+- `interp_methods`: Interpolation method used during GTC and RTC. Supported methods:
+ `nearest`, `bilinear`.
+- `footprint_scale_factor`: Defines how radar pixels contribute to the output grid.
+ Default: `(3.0, 15.0)`, reflecting the different scaling used for azimuth and
+ slant-range processing of SLC data.
+- `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:
+
+- [Docs – Sentinel-1 Analysis Mode](https://eopf-sample-service.github.io/xarray-eopf/examples/sentinel_1_analysis/)
+
##### Sentinel-1 Level-2 OCN
Sentinel-1 Level-2 OCN products are geolocated datasets provided on their
@@ -182,7 +250,7 @@ Sentinel-2 provides multi-spectral imagery at different native resolutions:
The analysis mode enables resampling between these different resolutions, bringing
bands from multiple resolutions onto the same grid using [affine transformation via xcube-resampling](https://xcube-dev.github.io/xcube-resampling/guide/#1-affine-transformation).
-**Suported Products:**
+**Supported Products:**
- [Sentinel-2 Level-1C](https://stac.browser.user.eopf.eodc.eu/collections/sentinel-2-l1c)
- [Sentinel-2 Level-2A](https://stac.browser.user.eopf.eodc.eu/collections/sentinel-2-l2a)
@@ -254,7 +322,7 @@ For OLCI products, no additional terrain correction is required, as it is alread
incorporated in the Level-1 data. See the [OLCI Level-1 product description](https://sentiwiki.copernicus.eu/web/olci-products#OLCIProducts-L1BProducts-ObservationModeS3-OLCI-Products-L1B-OM)
for details.
-**Suported Products:**
+**Supported Products:**
- [Sentinel-3 OLCI Level-1 EFR](https://stac.browser.user.eopf.eodc.eu/collections/sentinel-3-olci-l1-efr)
- [Sentinel-3 OLCI Level-1 ERR](https://stac.browser.user.eopf.eodc.eu/collections/sentinel-3-olci-l1-err)
=====================================
docs/index.md
=====================================
@@ -55,7 +55,7 @@ datatree = xr.open_datatree(url_or_path, engine="eopf-zarr", op_mode="native")
More information on the two modes are given in the [User Guide](guide.md)
Data variables will always be represented as chunked Dask arrays for
-efficient out-of core computations and visualisations.
+efficient out-of core computations and visualizations.
The package has minimal core dependencies: `xarray`, `zarr`, and `dask`.
Packages for accessing remote filesystems are optional, e.g., you will need `s3fs`
@@ -64,5 +64,5 @@ if you need to access EOPF data products in S3-compatible remote object storages
## License
The package is open source and released under the
-[Apache 2.0 license](https://www.apache.org/licenses/LICENSE-2.0.html) license. :heart:
+[Apache 2.0](https://www.apache.org/licenses/LICENSE-2.0.html) license. :heart:
=====================================
examples/introduction.ipynb
=====================================
The diff for this file was not included because it is too large.
=====================================
examples/sentinel_1_analysis.ipynb
=====================================
The diff for this file was not included because it is too large.
=====================================
examples/sentinel_1_native.ipynb
=====================================
The diff for this file was not included because it is too large.
=====================================
integration/test_data/dem_small.zarr.zip → integration/test_data/dem_grd.zarr.zip
=====================================
=====================================
integration/test_data/dem_slc.zarr.zip
=====================================
Binary files /dev/null and b/integration/test_data/dem_slc.zarr.zip differ
=====================================
integration/test_sen1_analysis.py
=====================================
@@ -15,22 +15,17 @@ allowed_open_time = 1000 # seconds
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"
+ dem_path = Path(__file__).resolve().parent / "test_data" / "dem_grd.zarr.zip"
store = zarr.ZipStore(str(dem_path), mode="r")
- dem = xr.open_zarr(
- store,
- group="dem_small.zarr",
- consolidated=False, # also important
- chunks={},
- )
+ dem = xr.open_zarr(store, group="dem_small.zarr", consolidated=False, chunks={})
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
@@ -50,6 +45,35 @@ class Sentinel2AnalysisTest(TestCase):
for var_name in ds.data_vars:
self.assertEqual((541, 1081), ds[var_name].shape, msg=var_name)
+ def test_open_dataset_sen1_slc(self):
+ dem_path = Path(__file__).resolve().parent / "test_data" / "dem_slc.zarr.zip"
+ store = zarr.ZipStore(str(dem_path), mode="r")
+ dem = xr.open_zarr(store, group="dem_small.zarr", consolidated=False, chunks={})
+ dem = dem.dem
+
+ url = (
+ "https://data.eodc.eu/collections/EOPF_ZARR/products/cpm_v270/"
+ "S01SIWSLC/2026/05/31/S1D_IW_SLC__1SDV_20260531T171503_"
+ "20260531T171530_003031_00539E_C573.zarr"
+ )
+ with timeit("open " + url) as result:
+ # noinspection PyTypeChecker
+ ds = xr.open_dataset(
+ url,
+ engine="eopf-zarr",
+ op_mode="analysis",
+ dem=dem,
+ chunks={},
+ )
+ self.assertTrue(result.time_delta < allowed_open_time)
+
+ self.assertIn("gamma0_vv", ds)
+ self.assertIn("gamma0_vh", ds)
+
+ assert_dataset_is_chunked(self, ds, verbose=show_chunking)
+ for var_name in ds.data_vars:
+ self.assertEqual((361, 361), ds[var_name].shape, msg=var_name)
+
def test_open_datatree_sen1_onc(self):
url = (
"https://objects.eodc.eu/e05ab01a9d56408d82ac32d69a5aae2a:202507-s01siwocn"
=====================================
integration/test_sen1_native.py
=====================================
@@ -10,19 +10,19 @@ import xarray as xr
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 = (
=====================================
mkdocs.yml
=====================================
@@ -2,7 +2,7 @@ site_name: EOPF xarray backend
repo_url: https://github.com/EOPF-Sample-Service/xarray-eopf
repo_name: EOPF-Sample-Service/xarray-eopf
-copyright: Copyright © 2025 EOPF Sample Service Team
+copyright: Copyright © 2026 EOPF Sample Service Team
nav:
- Overview: index.md
=====================================
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
@@ -14,18 +15,23 @@ import pytest
import xarray as xr
from xcube_resampling.gridmapping import GridMapping
-from tests.helpers import make_s1_grd_datatree, make_s1_ocn_datatree
+from tests.helpers import (
+ make_s1_grd_datatree,
+ make_s1_ocn_datatree,
+ make_s1_slc_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, Sen1SLC, register
class Sentinel1AnalysisModeTest(TestCase):
def test_register(self):
registry = AnalysisModeRegistry()
register(registry)
- self.assertEqual(2, len(list(registry.keys())))
+ self.assertEqual(3, len(list(registry.keys())))
self.assertIn(Sen1GRD.product_type, registry.keys())
+ self.assertIn(Sen1SLC.product_type, registry.keys())
self.assertIn(Sen1OCN.product_type, registry.keys())
@@ -53,9 +59,41 @@ class Sen1TestMixin:
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):
+ sen1._REGISTERED_CACHE_URIS.clear()
+ 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 tearDown(self):
+ sen1._REGISTERED_CACHE_URIS.clear()
def test_is_valid_source_ok(self):
self.assertTrue(self.mode.is_valid_source("data/S1A_IW_GRDH_20240201.zarr"))
@@ -66,66 +104,104 @@ class Sen1GRDTest(Sen1TestMixin, TestCase):
self.assertFalse(self.mode.is_valid_source(dict()))
def test_get_grid_parameters(self):
- params = self.mode._get_grid_parameters(self.dt, (2.0, 3.0))
+ params = sen1._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["range0"])
+ self.assertEqual(10.0, params["spacing_range"])
+ self.assertEqual(20.0, params["spacing_az"])
+ self.assertEqual(10.0, params["d_range"])
+ self.assertEqual(30.0, params["d_range_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(30.0, params["spacing_range_scale"])
+ 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"))
+
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="resolution argument must contain exactly two"
+ ):
+ self.mode.get_applicable_params(resolution=(1, "x"))
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(
- {"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=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.assertIs(args[0], self.dt)
+ self.assertEqual(["beta0_vv"], list(args[1].data_vars))
+ self.assertIs(args[2], self.dem)
self.assertEqual("bilinear", kwargs["interp_method"])
self.assertTrue(kwargs["apply_rtc"])
- self.assertIn("slr0", kwargs["grid_params"])
+
+ def test_convert_datatree_updates_footprint_scale_factor(self):
+ with patch.object(self.mode, "_terrain_correct", return_value=self.expected_vv):
+ self.mode.convert_datatree(
+ self.dt,
+ includes=["vv"],
+ dem=self.dem,
+ footprint_scale_factor=(2.0, 4.0),
+ )
+ self.assertEqual((2.0, 4.0), self.mode.footprint_scale_factor)
+
+ def test_convert_datatree_with_cache_uri_uses_fs(self):
+ 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=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.assertEqual(["file:///cache"], sen1._REGISTERED_CACHE_URIS)
def test_convert_datatree_uses_get_dem(self):
- expected = xr.Dataset(
- {"vv": xr.DataArray(np.ones((2, 2)), dims=("lat", "lon"))}
- )
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=self.expected_beta0_vv
+ ):
_ = self.mode.convert_datatree(self.dt, includes=["vv"])
get_dem_mock.assert_called_once()
@@ -136,6 +212,440 @@ class Sen1GRDTest(Sen1TestMixin, TestCase):
with pytest.raises(ValueError, match="No valid variable names"):
self.mode.convert_datatree(self.dt, includes="bibo", dem=self.dem)
+ def test_open_data_fails_for_duplicate_measurement_groups(self):
+ azimuth_time = np.array(
+ ["2024-01-01T00:00:00", "2024-01-01T00:00:01"], dtype="datetime64[ns]"
+ )
+ ground_range = np.array([0.0, 10.0])
+ invalid_dt = xr.DataTree.from_dict(
+ {
+ "S1A_IW_GRDH_TEST_VV/measurements": xr.Dataset(
+ {
+ "grd": xr.DataArray(
+ np.ones((2, 2), dtype="float32"),
+ dims=("azimuth_time", "ground_range"),
+ coords={
+ "azimuth_time": azimuth_time,
+ "ground_range": ground_range,
+ },
+ )
+ }
+ ),
+ "S1A_IW_GRDH_TEST_VV_DUP/measurements": xr.Dataset(
+ {
+ "grd": xr.DataArray(
+ np.ones((2, 2), dtype="float32"),
+ dims=("azimuth_time", "ground_range"),
+ coords={
+ "azimuth_time": azimuth_time,
+ "ground_range": ground_range,
+ },
+ )
+ }
+ ),
+ }
+ )
+
+ with pytest.raises(
+ ValueError,
+ match="Expected exactly one measurement group for polarization 'VV'",
+ ):
+ self.mode._open_data(invalid_dt, includes=["vv"])
+
+ def test_convert_datatree_warns_when_processing_full_product(self):
+ with (
+ self.assertWarns(UserWarning) as cm,
+ patch.object(sen1, "get_dem", return_value=self.dem),
+ patch.object(self.mode, "_terrain_correct", return_value=self.expected_vv),
+ ):
+ self.mode.convert_datatree(self.dt, includes=["vv"])
+ self.assertIn("No bounding box specified", str(cm.warning))
+
+ 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),
+ patch.object(xr.Dataset, "to_zarr", return_value=None),
+ patch.object(sen1.xr, "open_zarr", return_value=self.src_loc),
+ ):
+ self.mode.cache_uri = f"tmp_{uuid.uuid4().hex}"
+ out = self.mode._terrain_correct(
+ self.dt,
+ self.expected_beta0_vv,
+ self.dem,
+ 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),
+ patch.object(xr.Dataset, "to_zarr", return_value=None),
+ patch.object(sen1.xr, "open_zarr", return_value=self.src_loc),
+ ):
+ self.mode.cache_uri = f"tmp_{uuid.uuid4().hex}"
+ out = self.mode._terrain_correct(
+ self.dt,
+ self.expected_beta0_vv,
+ self.dem,
+ apply_rtc=True,
+ interp_method="bilinear",
+ )
+ gamma_mock.assert_called_once()
+ self.assertIn("gamma0_vv", out.data_vars)
+
+ def test_terrain_correct_without_rtc(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") as gamma_mock,
+ patch.object(sen1, "assign_grid_mapping", side_effect=lambda ds: ds),
+ patch.object(xr.Dataset, "to_zarr", return_value=None),
+ patch.object(sen1.xr, "open_zarr", return_value=self.src_loc),
+ ):
+ self.mode.cache_uri = f"tmp_{uuid.uuid4().hex}"
+ out = self.mode._terrain_correct(
+ self.dt,
+ self.expected_beta0_vv,
+ self.dem,
+ apply_rtc=False,
+ )
+ gamma_mock.assert_not_called()
+ self.assertIn("beta0_vv", out.data_vars)
+
+ def test_cleanup_registered_cache_uris_removes_all_caches(self):
+ sen1._register_cache_uri("file:///tmp/cache-a")
+ sen1._register_cache_uri("file:///tmp/cache-b")
+
+ fs = SimpleNamespace(exists=lambda path: True, rm=lambda path, recursive: None)
+ with patch.object(
+ sen1.fsspec,
+ "url_to_fs",
+ return_value=(fs, "/tmp/fake-cache"),
+ ) as url_to_fs:
+ sen1._cleanup_registered_cache_uris()
+
+ self.assertEqual(2, url_to_fs.call_count)
+ url_to_fs.assert_any_call("file:///tmp/cache-a")
+ url_to_fs.assert_any_call("file:///tmp/cache-b")
+ self.assertEqual(
+ ["file:///tmp/cache-a", "file:///tmp/cache-b"],
+ sen1._REGISTERED_CACHE_URIS,
+ )
+
+
+class Sen1SLCTest(Sen1TestMixin, TestCase):
+
+ def setUp(self):
+ sen1._REGISTERED_CACHE_URIS.clear()
+ self.mode = Sen1SLC()
+ self.dt = make_s1_slc_datatree()
+ 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.expected_beta0 = 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]"),
+ ),
+ "slant_range_time": (("lat", "lon"), np.zeros((2, 2))),
+ "gamma_area": (("lat", "lon"), np.ones((2, 2))),
+ }
+ )
+
+ def tearDown(self):
+ sen1._REGISTERED_CACHE_URIS.clear()
+
+ @staticmethod
+ def _make_slc_dataset(
+ azimuth_time: np.ndarray, slant_range_time: np.ndarray
+ ) -> xr.Dataset:
+ return xr.Dataset(
+ {
+ "beta0_vv": xr.DataArray(
+ np.ones((len(azimuth_time), len(slant_range_time))),
+ dims=("azimuth_time", "slant_range_time"),
+ coords={
+ "azimuth_time": azimuth_time,
+ "slant_range_time": slant_range_time,
+ },
+ )
+ }
+ )
+
+ @staticmethod
+ def _as_object_array(*datasets: xr.Dataset) -> np.ndarray:
+ out = np.empty(len(datasets), dtype=object)
+ for idx, dataset in enumerate(datasets):
+ out[idx] = dataset
+ return out
+
+ def test_is_valid_source_ok(self):
+ self.assertTrue(self.mode.is_valid_source("data/S1A_IW_SLC_20240201.zarr"))
+ self.assertTrue(self.mode.is_valid_source("S1D_SM_SLC_TEST"))
+
+ def test_is_not_valid_source(self):
+ self.assertFalse(self.mode.is_valid_source("data/S1A_IW_GRDH_20240201.zarr"))
+ self.assertFalse(self.mode.is_valid_source(dict()))
+
+ def test_get_groups(self):
+ groups = self.mode._get_groups(self.dt)
+ self.assertEqual(2, len(groups))
+ self.assertEqual(2, len(groups[0]))
+ self.assertEqual(3, len(groups[0][0]))
+ self.assertEqual(1, len(groups[0][0][0]))
+ self.assertEqual("S1A_IW_SLC_TEST_VV_IW1_0", groups[0][0][0][0])
+ self.assertEqual(["VH", "VV"], [str(mode) for mode in groups[1]])
+
+ def test_get_grid_parameters(self):
+ params = sen1._get_grid_parameters(
+ self.dt, (2.0, 3.0), range_coord="slant_range_time"
+ )
+ self.assertEqual(0.0, params["range0"])
+ self.assertEqual(1.0, params["d_range"])
+ self.assertEqual(10.0, params["spacing_range"])
+ self.assertEqual(3.0, params["d_range_scale"])
+ self.assertEqual(30.0, params["spacing_range_scale"])
+ self.assertEqual(1.0, params["range0_scale"])
+
+ def test_calibrate_burst_and_extract_valid_region(self):
+ burst = self.dt["S1A_IW_SLC_TEST_VV_IW1_0"]
+ beta0 = self.mode._calibrate_burst(burst)
+ trimmed = self.mode._extract_valid_region(beta0, burst)
+ self.assertEqual(("azimuth_time", "slant_range_time"), trimmed.slc.dims)
+ self.assertEqual(2, trimmed.sizes["azimuth_time"])
+ self.assertEqual(2, trimmed.sizes["slant_range_time"])
+
+ def test_open_data(self):
+ out = self.mode._open_data(self.dt, includes=["vv", "vh"])
+ self.assertCountEqual(["beta0_vv", "beta0_vh"], out.data_vars)
+ self.assertNotIn("line", out.coords)
+ self.assertNotIn("pixel", out.coords)
+ self.assertEqual({"azimuth_time": 2, "slant_range_time": 4}, out.sizes)
+
+ def test_open_data_accepts_string_include(self):
+ out = self.mode._open_data(self.dt, includes="vv")
+ self.assertEqual(["beta0_vv"], list(out.data_vars))
+
+ def test_open_data_fails_when_no_variables_match(self):
+ with pytest.raises(
+ ValueError, match="No valid variable names found in dataset"
+ ):
+ self.mode._open_data(self.dt, includes="does_not_exist")
+
+ def test_convert_datatree(self):
+ with patch.object(
+ self.mode, "_terrain_correct", return_value=self.expected_beta0
+ ) as mocked:
+ out = self.mode.convert_datatree(self.dt, includes=["vv"], dem=self.dem)
+ self.assertIs(out, self.expected_beta0)
+ args, kwargs = mocked.call_args
+ self.assertIs(args[0], self.dt)
+ self.assertEqual(["beta0_vv"], list(args[1].data_vars))
+ self.assertIs(args[2], self.dem)
+ self.assertEqual("bilinear", kwargs["interp_method"])
+
+ def test_terrain_correct_uses_slant_range_path(self):
+ with (
+ patch.object(
+ sen1, "get_source_location", return_value=self.src_loc
+ ) as src_mock,
+ patch.object(sen1, "geocode_data", return_value=self.expected_beta0),
+ 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),
+ patch.object(xr.Dataset, "to_zarr", return_value=None),
+ patch.object(sen1.xr, "open_zarr", return_value=self.src_loc),
+ ):
+ self.mode.cache_uri = f"tmp_{uuid.uuid4().hex}"
+ out = self.mode._terrain_correct(
+ self.dt,
+ self.expected_beta0,
+ self.dem,
+ apply_rtc=True,
+ interp_method="nearest",
+ )
+ self.assertIn("gamma0_vv", out.data_vars)
+ self.assertIsNone(src_mock.call_args.kwargs["time_slr_gcp"])
+ self.assertEqual("slant_range_time", src_mock.call_args.kwargs["range_coord"])
+ self.assertEqual("slant_range_time", gamma_mock.call_args.kwargs["range_coord"])
+
+ def test_merge_bursts_warns_on_irregular_spacing(self):
+ ds0 = self._make_slc_dataset(
+ np.array(
+ [
+ "2024-01-01T00:00:00",
+ "2024-01-01T00:00:01",
+ "2024-01-01T00:00:02",
+ ],
+ dtype="datetime64[ns]",
+ ),
+ np.array([0.0, 1.0]),
+ )
+ ds1 = self._make_slc_dataset(
+ np.array(
+ [
+ "2024-01-01T00:00:01.800000000",
+ "2024-01-01T00:00:03.800000000",
+ "2024-01-01T00:00:05.800000000",
+ ],
+ dtype="datetime64[ns]",
+ ),
+ np.array([0.0, 1.0]),
+ )
+ with pytest.warns(UserWarning, match="Azimuth time spacing is not regular"):
+ out = self.mode._merge_bursts(self._as_object_array(ds0, ds1))
+ self.assertEqual(5, out.sizes["azimuth_time"])
+
+ def test_merge_bursts_raises_without_overlap(self):
+ ds0 = self._make_slc_dataset(
+ np.array(
+ ["2024-01-01T00:00:00", "2024-01-01T00:00:01"],
+ dtype="datetime64[ns]",
+ ),
+ np.array([0.0, 1.0]),
+ )
+ ds1 = self._make_slc_dataset(
+ np.array(
+ ["2024-01-01T00:00:10", "2024-01-01T00:00:11"],
+ dtype="datetime64[ns]",
+ ),
+ np.array([0.0, 1.0]),
+ )
+ with pytest.raises(ValueError, match="No overlap found"):
+ self.mode._merge_bursts(self._as_object_array(ds0, ds1))
+
+ def test_align_azimuth_warns_on_irregular_spacing(self):
+ ds0 = self._make_slc_dataset(
+ np.array(
+ [
+ "2024-01-01T00:00:00",
+ "2024-01-01T00:00:01",
+ "2024-01-01T00:00:02",
+ "2024-01-01T00:00:03",
+ ],
+ dtype="datetime64[ns]",
+ ),
+ np.array([0.0, 1.0]),
+ )
+ ds1 = self._make_slc_dataset(
+ np.array(
+ [
+ "2024-01-01T00:00:00",
+ "2024-01-01T00:00:01.100000000",
+ "2024-01-01T00:00:02.100000000",
+ "2024-01-01T00:00:03.100000000",
+ ],
+ dtype="datetime64[ns]",
+ ),
+ np.array([0.0, 1.0]),
+ )
+ with pytest.warns(UserWarning, match="Azimuth spacing is irregular"):
+ out = self.mode._align_azimuth(self._as_object_array(ds0, ds1))
+ self.assertEqual(4, out[0].sizes["azimuth_time"])
+ self.assertTrue(
+ np.array_equal(out[0].azimuth_time.values, out[1].azimuth_time.values)
+ )
+
+ def test_align_azimuth_warns_on_size_difference(self):
+ ds0 = self._make_slc_dataset(
+ np.array(
+ [
+ "2024-01-01T00:00:00",
+ "2024-01-01T00:00:01",
+ "2024-01-01T00:00:02",
+ "2024-01-01T00:00:03",
+ "2024-01-01T00:00:04",
+ ],
+ dtype="datetime64[ns]",
+ ),
+ np.array([0.0, 1.0]),
+ )
+ ds1 = self._make_slc_dataset(
+ np.array(
+ [
+ "2024-01-01T00:00:01",
+ "2024-01-01T00:00:02",
+ "2024-01-01T00:00:03.300000000",
+ "2024-01-01T00:00:04.300000000",
+ ],
+ dtype="datetime64[ns]",
+ ),
+ np.array([0.0, 1.0]),
+ )
+ with pytest.warns(
+ UserWarning, match="Aligned swaths have different azimuth sizes"
+ ):
+ out = self.mode._align_azimuth(self._as_object_array(ds0, ds1))
+ self.assertEqual(3, out[0].sizes["azimuth_time"])
+
+ def test_merge_swaths_warns_on_irregular_spacing(self):
+ ds0 = self._make_slc_dataset(
+ np.array(
+ ["2024-01-01T00:00:00", "2024-01-01T00:00:01"],
+ dtype="datetime64[ns]",
+ ),
+ np.array([0.0, 1.0, 2.0]),
+ )
+ ds1 = self._make_slc_dataset(
+ np.array(
+ ["2024-01-01T00:00:00", "2024-01-01T00:00:01"],
+ dtype="datetime64[ns]",
+ ),
+ np.array([1.8, 3.8, 5.8]),
+ )
+ with pytest.warns(UserWarning, match="Slant range spacing is irregular"):
+ out = self.mode._merge_swaths([ds0, ds1])
+ self.assertEqual(5, out.sizes["slant_range_time"])
+
+ def test_merge_swaths_raises_without_overlap(self):
+ ds0 = self._make_slc_dataset(
+ np.array(
+ ["2024-01-01T00:00:00", "2024-01-01T00:00:01"],
+ dtype="datetime64[ns]",
+ ),
+ np.array([0.0, 1.0]),
+ )
+ ds1 = self._make_slc_dataset(
+ np.array(
+ ["2024-01-01T00:00:00", "2024-01-01T00:00:01"],
+ dtype="datetime64[ns]",
+ ),
+ np.array([10.0, 11.0]),
+ )
+ with pytest.raises(ValueError, match="No overlap found"):
+ self.mode._merge_swaths([ds0, ds1])
+
class Sen1OCNTest(Sen1TestMixin, TestCase):
mode = Sen1OCN()
@@ -167,6 +677,10 @@ class Sen1OCNTest(Sen1TestMixin, TestCase):
agg_methods="nearest",
),
)
+ with pytest.raises(
+ TypeError, match="resolution argument must contain exactly two"
+ ):
+ self.mode.get_applicable_params(resolution=(1, "x"))
with pytest.raises(TypeError):
self.mode.get_applicable_params(interp_methods="cubic")
@@ -227,48 +741,100 @@ class Sen1OCNTest(Sen1TestMixin, TestCase):
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(
+ range0=0.0,
+ range0_scale=0.0,
+ d_range=1.0,
+ d_range_scale=1.0,
+ spacing_range=1.0,
+ spacing_range_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],
+ "spatial_ref": xr.DataArray(
+ 0, attrs=pyproj.CRS.from_epsg(4326).to_cf()
+ ),
+ },
)
-
- 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),
+ )
+ 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"]},
)
- 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(
+ [
+ "range0",
+ "range0_scale",
+ "d_range",
+ "d_range_scale",
+ "spacing_range",
+ "spacing_range_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_gridparams_contains(self):
+ self.assertIn("range0", self.grid_params)
+ self.assertNotIn("gr0", self.grid_params)
def test_get_dem_requires_credentials(self):
with patch.dict(os.environ, {}, clear=True):
with pytest.raises(ValueError, match="Missing AWS credentials"):
sen1.get_dem([0, 50, 1, 51])
- def test_get_dem_resolution_required_if_crs_given(self):
+ def test_get_dem_with_projected_crs_uses_inferred_resolution(self):
with patch.dict(
os.environ,
{"AWS_ACCESS_KEY_ID": "k", "AWS_SECRET_ACCESS_KEY": "s"},
@@ -277,6 +843,18 @@ class Sentinel1FunctionsTest(TestCase):
with (
patch.object(sen1.pystac_client.Client, "open") as client_open,
patch.object(sen1.rioxarray, "open_rasterio") as open_rasterio,
+ patch.object(sen1, "transform_resolution", return_value=30.0) as tr,
+ patch.object(
+ sen1,
+ "resample_in_space",
+ return_value=SimpleNamespace(
+ dem=xr.DataArray(
+ np.ones((2, 2), dtype="float32"),
+ dims=("lat", "lon"),
+ coords={"lat": [1.0, 0.0], "lon": [0.0, 1.0]},
+ )
+ ),
+ ) as resample,
):
fake_item = SimpleNamespace(assets={"data": SimpleNamespace(href="x")})
search = SimpleNamespace(items=lambda: [fake_item])
@@ -287,10 +865,12 @@ class Sentinel1FunctionsTest(TestCase):
coords={"band": [1], "y": [3, 2, 1, 0], "x": [0, 1, 2, 3]},
)
- with pytest.raises(
- ValueError, match="Resolution must be provided if CRS is not None"
- ):
- sen1.get_dem([0, 0, 1, 1], crs=pyproj.CRS.from_epsg(32632))
+ out = sen1.get_dem([0, 0, 900, 900], crs=pyproj.CRS.from_epsg(32632))
+
+ tr.assert_called_once()
+ resample.assert_called_once()
+ self.assertIsInstance(out, xr.DataArray)
+ self.assertEqual((2, 2), out.shape)
def test_get_dem_reprojects_bbox_and_resamples(self):
with patch.dict(
@@ -392,7 +972,10 @@ class Sentinel1FunctionsTest(TestCase):
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"])
@@ -419,27 +1002,127 @@ class Sentinel1FunctionsTest(TestCase):
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):
@@ -482,40 +1165,28 @@ class Sentinel1FunctionsTest(TestCase):
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,96 +1197,35 @@ class Sentinel1FunctionsTest(TestCase):
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(
- 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"))
area = xr.DataArray(
np.array(
[
@@ -625,7 +1235,7 @@ class Sentinel1FunctionsTest(TestCase):
]
),
dims=("axis", "lat", "lon"),
- coords=dem_ecef.coords,
+ coords=self.dem_ecef.coords,
)
direction = xr.DataArray(
np.array(
@@ -636,10 +1246,12 @@ class Sentinel1FunctionsTest(TestCase):
]
),
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]))
@@ -651,79 +1263,83 @@ class Sentinel1FunctionsTest(TestCase):
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):
- 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")),
- "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(
+ 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):
- 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(
+ range0=0.0,
+ range0_scale=0.0,
+ d_range=1.0,
+ d_range_scale=1.0,
+ spacing_range=1.0,
+ spacing_range_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 +1361,51 @@ class Sentinel1FunctionsTest(TestCase):
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(
+ range0=0.0,
+ range0_scale=0.0,
+ d_range=3.0,
+ d_range_scale=3.0,
+ spacing_range=3.0,
+ spacing_range_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)
=====================================
tests/amodes/test_sentinel3.py
=====================================
@@ -49,6 +49,12 @@ class Sen3TestMixin:
agg_methods={"scl": "mode"},
),
)
+ with pytest.raises(
+ TypeError, match="resolution argument must contain exactly two"
+ ):
+ self.mode.get_applicable_params(resolution=(1, "x"))
+ with pytest.raises(TypeError):
+ self.mode.get_applicable_params(interp_methods=["nearest"])
def test_process_metadata(self: TestCase):
self.assertEqual({}, self.mode.process_metadata(xr.DataTree()))
=====================================
tests/helpers/__init__.py
=====================================
@@ -2,13 +2,14 @@
# Permissions are hereby granted under the terms of the Apache 2.0 License:
# https://opensource.org/license/apache-2-0.
-from .sentinel1 import make_s1_grd_datatree, make_s1_ocn_datatree
+from .sentinel1 import make_s1_grd_datatree, make_s1_ocn_datatree, make_s1_slc_datatree
from .sentinel2 import make_s2_msi, make_s2_msi_l1c, make_s2_msi_l2a
from .sentinel3 import make_s3_olci_efr, make_s3_slstr_lst, make_s3_slstr_rbt
__all__ = [
"make_s1_grd_datatree",
"make_s1_ocn_datatree",
+ "make_s1_slc_datatree",
"make_s2_msi",
"make_s2_msi_l1c",
"make_s2_msi_l2a",
=====================================
tests/helpers/sentinel1.py
=====================================
@@ -86,6 +86,145 @@ def make_s1_grd_datatree() -> xr.DataTree:
return dt
+def _make_slc_group(
+ azimuth_time: np.ndarray,
+ slant_range_time: np.ndarray,
+ scale: float,
+) -> xr.Dataset:
+ values = np.arange(
+ 1, len(azimuth_time) * len(slant_range_time) + 1, dtype="float32"
+ ).reshape(len(azimuth_time), len(slant_range_time))
+ return xr.Dataset(
+ {
+ "slc": xr.DataArray(
+ values * scale,
+ dims=("azimuth_time", "slant_range_time"),
+ coords={
+ "azimuth_time": azimuth_time,
+ "slant_range_time": slant_range_time,
+ "line": ("azimuth_time", np.arange(len(azimuth_time))),
+ "pixel": ("slant_range_time", np.arange(len(slant_range_time))),
+ },
+ )
+ }
+ )
+
+
+def make_s1_slc_datatree() -> xr.DataTree:
+ azimuth_time_by_swath = {
+ "IW1": np.array(
+ [
+ "2024-01-01T00:00:00",
+ "2024-01-01T00:00:01",
+ "2024-01-01T00:00:02",
+ ],
+ dtype="datetime64[ns]",
+ ),
+ "IW2": np.array(
+ [
+ "2024-01-01T00:00:00",
+ "2024-01-01T00:00:01",
+ "2024-01-01T00:00:02",
+ ],
+ dtype="datetime64[ns]",
+ ),
+ "IW3": np.array(
+ [
+ "2024-01-01T00:00:00",
+ "2024-01-01T00:00:01",
+ "2024-01-01T00:00:02",
+ ],
+ dtype="datetime64[ns]",
+ ),
+ }
+ slant_range_time_by_swath = {
+ "IW1": np.array([0.0, 1.0, 2.0, 3.0], dtype="float32"),
+ "IW2": np.array([1.0, 2.0, 3.0, 4.0], dtype="float32"),
+ "IW3": np.array([2.0, 3.0, 4.0, 5.0], dtype="float32"),
+ }
+ axis = np.array(["x", "y", "z"])
+ polarizations = ("VV", "VH")
+ swaths = ("IW1", "IW2", "IW3")
+
+ beta_nought = xr.Dataset(
+ {
+ "beta_nought": xr.DataArray(
+ np.ones((3, 4), dtype="float32"),
+ dims=("azimuth_time", "slant_range_time"),
+ coords={
+ "azimuth_time": azimuth_time_by_swath["IW1"],
+ "slant_range_time": slant_range_time_by_swath["IW1"],
+ },
+ )
+ }
+ )
+ orbit = xr.Dataset(
+ {
+ "position": xr.DataArray(
+ np.ones((3, 3), dtype="float32"),
+ dims=("azimuth_time", "axis"),
+ coords={
+ "azimuth_time": azimuth_time_by_swath["IW1"],
+ "axis": axis,
+ },
+ )
+ }
+ )
+ burst_info = xr.Dataset(
+ {
+ "first_valid_sample": xr.DataArray(
+ np.array([1, 1, -1], dtype="int32"),
+ dims=("azimuth_time",),
+ coords={"azimuth_time": azimuth_time_by_swath["IW1"]},
+ ),
+ "last_valid_sample": xr.DataArray(
+ np.array([3, 3, -1], dtype="int32"),
+ dims=("azimuth_time",),
+ coords={"azimuth_time": azimuth_time_by_swath["IW1"]},
+ ),
+ }
+ )
+
+ dt_nodes = {}
+ for mode_i, pol in enumerate(polarizations):
+ for swath_i, swath in enumerate(swaths):
+ group_name = f"S1A_IW_SLC_TEST_{pol}_{swath}_0"
+ dt_nodes.update(
+ {
+ f"{group_name}/measurements": _make_slc_group(
+ azimuth_time_by_swath[swath],
+ slant_range_time_by_swath[swath],
+ scale=float(mode_i + swath_i + 1),
+ ),
+ f"{group_name}/quality/calibration": beta_nought.assign_coords(
+ azimuth_time=azimuth_time_by_swath[swath],
+ slant_range_time=slant_range_time_by_swath[swath],
+ ).copy(deep=True),
+ f"{group_name}/conditions/orbit": orbit.assign_coords(
+ azimuth_time=azimuth_time_by_swath[swath]
+ ).copy(deep=True),
+ f"{group_name}/conditions/burst_info": burst_info.assign_coords(
+ azimuth_time=azimuth_time_by_swath[swath]
+ ).copy(deep=True),
+ }
+ )
+
+ dt = xr.DataTree.from_dict(dt_nodes)
+ dt.attrs["stac_discovery"] = {"bbox": [0.0, 50.0, 1.0, 51.0]}
+ dt["S1A_IW_SLC_TEST_VH_IW1_0"].attrs["other_metadata"] = {
+ "image_annotation": {
+ "image_information": {
+ "range_pixel_spacing": 10.0,
+ "slant_range_time": 1.0e-4,
+ "product_first_line_utc_time": "2024-01-01T00:00:00",
+ "azimuth_time_interval": 0.5,
+ "azimuth_pixel_spacing": 20.0,
+ }
+ }
+ }
+ return dt
+
+
def _make_owi_measurements(height: int, width: int):
shape = (height, width)
=====================================
xarray_eopf/amodes/sentinel1.py
=====================================
@@ -2,18 +2,20 @@
# Permissions are hereby granted under the terms of the Apache 2.0 License:
# https://opensource.org/license/apache-2-0.
-
import functools
import os
import re
+import uuid
import warnings
+import weakref
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 dask.array as da
import flox.xarray
+import fsspec
import numpy as np
import pyproj
import pystac_client
@@ -23,9 +25,17 @@ from xcube_resampling import resample_in_space
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
+
+# noinspection PyProtectedMember
+from xcube_resampling.utils import (
+ SourceTileIndexing,
+ _reorganize_tiled_array,
+ reproject_bbox,
+ transform_resolution,
+)
from xarray_eopf.amode import AnalysisMode, AnalysisModeRegistry
+from xarray_eopf.constants import FloatInt
from xarray_eopf.source import get_source_path
from xarray_eopf.utils import NameFilter, assert_arg_has_length, assert_arg_is_instance
@@ -36,66 +46,36 @@ _CRS_ECEF = pyproj.CRS.from_string("EPSG:4978")
_CRS_WGS84 = pyproj.CRS.from_string("EPSG:4326")
_DEM_CHUNKSIZE = dict(lat=1800, lon=1800)
_CHUNKSIZE = (2048, 2048)
+_SLC_SWATHS = ["IW1", "IW2", "IW3"]
+_SENTINEL1_POLARZIATION_MODES = ["VV", "VH", "HV", "HH"]
+_REGISTERED_CACHE_URIS = []
@dataclass(frozen=True)
class GridParams:
- """RTC grid parameters."""
-
- slr0: float
- d_slr: float
- spacing_slr: float
+ """Range and azimuth grid parameters used for geocoding and RTC."""
+
+ range0: float
+ range0_scale: float
+ d_range: float
+ d_range_scale: float
+ spacing_range: float
+ spacing_range_scale: float
az0: np.datetime64
+ az0_scale: 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"}
-
-
- at 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,13 +104,24 @@ class Sen1(AnalysisMode, ABC):
class Sen1GRD(Sen1):
product_type = "GRDH"
+ range_coord = "ground_range"
+ footprint_scale_factor = (3.0, 3.0)
+ cache_fs: fsspec.AbstractFileSystem | None = None
+ cache_uri: str | None = None
def get_applicable_params(self, **kwargs) -> dict[str, Any]:
params = {}
resolution = kwargs.get("resolution")
if resolution is not None:
- assert_arg_is_instance(resolution, "resolution", (float, int))
+ assert_arg_is_instance(resolution, "resolution", (float, int, tuple))
+ if isinstance(resolution, tuple):
+ assert_arg_has_length(resolution, "resolution", 2)
+ if not all(isinstance(v, (float, int)) for v in resolution):
+ raise TypeError(
+ "resolution argument must contain exactly "
+ "two float or int values."
+ )
params.update(resolution=resolution)
bbox = kwargs.get("bbox")
@@ -178,6 +169,11 @@ class Sen1GRD(Sen1):
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(
@@ -185,17 +181,36 @@ class Sen1GRD(Sen1):
datatree: xr.DataTree,
includes: str | Iterable[str] | None = None,
excludes: str | Iterable[str] | None = None,
- resolution: float = None,
+ resolution: FloatInt | tuple[FloatInt, FloatInt] | None = None,
bbox: Sequence[float | int] | None = None,
crs: pyproj.CRS | None = None,
interp_methods: Literal["nearest", "bilinear"] = "bilinear",
- footprint_scale_factor: tuple[float, float] = (3.0, 3.0),
+ footprint_scale_factor: tuple[float, float] | None = None,
dem: xr.DataArray | None = None,
apply_rtc: bool = True,
+ cache_uri: str | None = None,
) -> xr.Dataset:
- # get dem data array
+ if footprint_scale_factor is not None:
+ self.footprint_scale_factor = footprint_scale_factor
+
+ 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
+ _register_cache_uri(self.cache_uri)
+ weakref.finalize(self, _cleanup_registered_cache_uris)
+
+ # Build the DEM if one was not supplied. When no bbox was passed,
+ # fall back to the product footprint from STAC metadata.
if dem is None:
if bbox is None:
+ warnings.warn(
+ f"No bounding box specified. Processing the full Sentinel-1 "
+ f"{self.product_type} product, which may take some time."
+ )
bbox = datatree.attrs["stac_discovery"]["bbox"]
bbox = [
min(bbox[0], bbox[2]),
@@ -205,93 +220,392 @@ class Sen1GRD(Sen1):
]
dem = get_dem(bbox, resolution=resolution, crs=crs)
- # load measurement data
- grd = None
- group = ""
- for mode in ["VV", "VH", "HV", "HH"]:
+ grd = self._open_data(datatree, includes, excludes)
+ return self._terrain_correct(
+ datatree,
+ grd,
+ dem,
+ apply_rtc=apply_rtc,
+ interp_method=interp_methods,
+ )
+
+ @staticmethod
+ def _open_data(
+ datatree: xr.DataTree,
+ includes: str | Iterable[str] | None = None,
+ excludes: str | Iterable[str] | None = None,
+ ) -> xr.Dataset:
+ """Load measurements from the input datatree and apply radiometric
+ calibration using look-up tables provided in the product.
+ """
+
+ dataset = None
+ measurement_group = ""
+ for mode in _SENTINEL1_POLARZIATION_MODES:
children = [x for x in datatree.children if mode in x]
if children:
- group = children[0]
- if grd is None:
- grd = datatree[group].measurements.to_dataset()
- grd = grd.rename({"grd": mode.lower()})
+ if len(children) != 1:
+ raise ValueError(
+ f"Expected exactly one measurement group for polarization "
+ f"'{mode}', but found {len(children)}: {children}. "
+ "The Sentinel-1 Level-1 GRD product appears to be invalid."
+ )
+ measurement_group = children[0]
+ if dataset is None:
+ dataset = datatree[measurement_group].measurements.to_dataset()
+ dataset = dataset.rename({"grd": mode.lower()})
else:
- grd[mode.lower()] = datatree[group].measurements.to_dataset().grd
+ dataset[mode.lower()] = (
+ datatree[measurement_group].measurements.to_dataset().grd
+ )
# filter dataset by variable names
name_filter = NameFilter(includes=includes, excludes=excludes)
- variable_names = [k for k in grd.data_vars if name_filter.accept(str(k))]
+ variable_names = [k for k in dataset.data_vars if name_filter.accept(str(k))]
if not variable_names:
raise ValueError("No valid variable names found in dataset")
- grd = grd[variable_names]
+ dataset = dataset[variable_names]
- # get calibration LUT data
- lut = datatree[group].quality.calibration.beta_nought
- lut_interp = lut.interp(ground_range=grd.ground_range).chunk(
- dict(ground_range=2048)
- )
- lut_interp = lut_interp.interp(azimuth_time=grd.azimuth_time).chunk(
- dict(azimuth_time=2048)
- )
- grd = (grd / lut_interp) ** 2
+ # Mask invalid border pixels (zero-valued fill areas) as NaN
+ dataset = dataset.where(dataset != 0)
+
+ # convert amplitudes to beta nought backscatter.
+ beta_lut = datatree[measurement_group].quality.calibration.beta_nought.values
+ assert len(np.unique(beta_lut)) == 1
+ dataset = (dataset / beta_lut[0, 0]) ** 2
rename_dict = {name: f"beta0_{name}" for name in variable_names}
- grd = grd.rename(rename_dict)
- for var in grd.data_vars:
- grd[var].attrs.update(
+ dataset = dataset.rename(rename_dict)
+ for var in dataset.data_vars:
+ dataset[var].attrs.update(
long_name="beta nought backscatter coefficient",
units="1",
)
+ return dataset
+ def _terrain_correct(
+ self,
+ datatree: xr.DataTree,
+ data: xr.Dataset,
+ dem: xr.DataArray,
+ apply_rtc: bool = True,
+ interp_method: Literal["nearest", "bilinear"] = "nearest",
+ ) -> xr.Dataset:
+ """Geocode the calibrated GRD data and optionally apply
+ radiometric terrain correction (RTC).
+ """
+ group = next(iter(datatree.children))
orbit = datatree[f"{group}/conditions/orbit"].to_dataset()
- sat_position = orbit["position"]
-
- 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)
+ sat_position = orbit["position"].compute()
+
+ if self.range_coord == "ground_range":
+ gcp = datatree[f"{group}/conditions/gcp"].to_dataset()
+ time_slr_gcp = gcp["slant_range_time_gcp"]
+ else:
+ time_slr_gcp = None
+
+ grid_params = _get_grid_parameters(
+ datatree,
+ self.footprint_scale_factor,
+ range_coord=self.range_coord,
+ )
- return terrain_correct(
- grd,
- time_slr_gcp,
- sat_position,
- dem,
+ gm_dem = GridMapping.from_dataset(dem.to_dataset(name="dem"))
+ src_loc = get_source_location(
+ dem=dem,
+ time_slr_gcp=time_slr_gcp,
+ sat_position=sat_position,
grid_params=grid_params,
+ gm_dem_grid=gm_dem,
apply_rtc=apply_rtc,
- interp_method=interp_methods,
+ range_coord=self.range_coord,
)
+ 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, range_coord=self.range_coord
+ )
+
+ 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,
+ range_coord=self.range_coord,
+ )
+ geocoded /= gamma_weights
+ rename_dict = {
+ name: str(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
+
+
+class Sen1SLC(Sen1GRD):
+ product_type = "SLC"
+ range_coord = "slant_range_time"
+ footprint_scale_factor = (3.0, 15.0)
+ cache_fs: fsspec.AbstractFileSystem | None = None
+ cache_uri: str | None = None
+
+ def _open_data(
+ self,
+ datatree: xr.DataTree,
+ includes: str | Iterable[str] | None = None,
+ excludes: str | Iterable[str] | None = None,
+ ) -> xr.Dataset:
+ """Load, calibrate, and merge SLC bursts into a single dataset."""
+
+ children, modes = self._get_groups(datatree)
+
+ # First combine all polarizations within each burst.
+ dss_all = []
+ for swath_i, _ in enumerate(_SLC_SWATHS):
+ dss_swath = []
+ for burst_i, _ in enumerate(children[0][swath_i]):
+ burst_dataset = xr.Dataset()
+ for mode_i, mode in enumerate(modes):
+ child = children[mode_i][swath_i][burst_i]
+ burst = datatree[child]
+ beta0 = self._calibrate_burst(burst)
+ beta0 = self._extract_valid_region(beta0, burst)
+ beta0 = beta0.drop_vars(["line", "pixel"])
+ beta0 = beta0.rename({"slc": f"beta0_{mode.lower()}"})
+ burst_dataset.update(beta0)
+ dss_swath.append(burst_dataset)
+ dss_all.append(dss_swath)
+
+ dss_swaths = []
+ for swath_i, swath in enumerate(_SLC_SWATHS):
+ dss_swaths.append(self._merge_bursts(dss_all[swath_i]))
+ dss_swaths = self._align_azimuth(dss_swaths)
+ merged = self._merge_swaths(dss_swaths)
+
+ def _expand_names(
+ names: str | Iterable[str] | None,
+ ) -> str | list[str] | None:
+ if names is None:
+ return None
+ if isinstance(names, str):
+ names = [names]
+ expanded = []
+ for name in names:
+ expanded.append(name)
+ if not str(name).startswith("beta0_"):
+ expanded.append(f"beta0_{name}")
+ return expanded
+
+ name_filter = NameFilter(
+ includes=_expand_names(includes),
+ excludes=_expand_names(excludes),
+ )
+ variable_names = [k for k in merged.data_vars if name_filter.accept(str(k))]
+ if not variable_names:
+ raise ValueError("No valid variable names found in dataset")
+ merged = merged[variable_names]
+ return merged
@staticmethod
- def _get_grid_parameters(
- dt: xr.DataTree,
- footprint_scale_factor: tuple[float, float],
- ) -> GridParams:
- """Build grid parameters for RTC from Sentinel-1 metadata.
-
- Args:
- dt: Source data tree.
- footprint_scale_factor: Scaling for SAR footprint spacing.
-
- Returns:
- Grid parameters for terrain correction.
- """
+ def _get_groups(datatree: xr.DataTree) -> tuple[list[list[list[str]]], list[str]]:
+ """Return child group names organized by polarization, swath, and burst."""
+ modes_sel = []
+ children = []
+ for mode in _SENTINEL1_POLARZIATION_MODES:
+ mode_children = []
+ for swath_i, swath in enumerate(_SLC_SWATHS):
+ bursts = [x for x in datatree.children if f"_{mode}_{swath}" in x]
+ if bursts:
+ mode_children.append(bursts)
+ modes_sel.append(mode)
+ if mode_children:
+ children.append(mode_children)
+ return children, list(np.unique(modes_sel))
- 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"
+ @staticmethod
+ def _calibrate_burst(burst: xr.DataTree) -> xr.Dataset:
+ """Convert a burst measurement to beta nought."""
+ slc = burst.measurements.to_dataset()
+ beta_lut = burst.quality.calibration.beta_nought
+ assert len(np.unique(beta_lut)) == 1
+ return (abs(slc) / beta_lut[0, 0].values) ** 2
+
+ @staticmethod
+ def _extract_valid_region(dataset: xr.Dataset, burst: xr.DataTree) -> xr.Dataset:
+ """Trim invalid burst samples using the burst information mask."""
+ first = burst.conditions.burst_info.first_valid_sample.values
+ last = burst.conditions.burst_info.last_valid_sample.values
+ first_idxs = np.unique(first)
+ assert len(first_idxs) == 2
+ assert -1 in first_idxs
+ valid = np.where(first != -1)[0]
+ first_idx = first_idxs[1]
+ last_idxs = np.unique(last[valid])
+ assert len(last_idxs) == 1
+ last_idx = last_idxs[0]
+ return dataset.isel(
+ azimuth_time=valid, slant_range_time=slice(first_idx, last_idx)
+ )
+
+ @staticmethod
+ def _merge_bursts(dss: list[xr.Dataset], tolerance: float = 0.01) -> xr.Dataset:
+ """Merge overlapping bursts along azimuth time."""
+ idxs = np.zeros((len(dss), 2), dtype=int)
+ tol = 0.01 * np.diff(dss[0].azimuth_time.values[:2])[0]
+
+ for i, (ds0, ds1) in enumerate(zip(dss[:-1], dss[1:])):
+
+ t0 = ds0.azimuth_time.values
+ t1 = ds1.azimuth_time.values
+
+ mask0 = t0 >= t1[0] - tol
+ mask1 = t1 <= t0[-1] + tol
+ n_overlap = min(mask0.sum(), mask1.sum())
+
+ if n_overlap == 0:
+ raise ValueError(f"No overlap found between swaths {i} and {i + 1}")
+
+ keep0 = n_overlap // 2
+ keep1 = n_overlap - keep0
+
+ idxs[i, 1] = len(t0) - keep0
+ idxs[i + 1, 0] = keep1
+
+ idxs[-1, 1] = dss[-1].sizes["azimuth_time"]
+
+ dss_cut = [
+ ds.isel(azimuth_time=slice(start, stop))
+ for ds, (start, stop) in zip(dss, idxs)
]
- 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
+ out = xr.concat(
+ dss_cut,
+ dim="azimuth_time",
+ join="outer",
+ coords="minimal",
+ )
- return GridParams(
- slr0=attrs["slant_range_time"],
- d_slr=slant_range_time_interval_s,
- spacing_slr=slant_range_spacing_m,
- 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],
+ dt = out.azimuth_time.diff("azimuth_time").values
+ nominal_dt = np.median(dt)
+ relative_deviation = np.abs(dt - nominal_dt) / nominal_dt
+ if any(relative_deviation > tolerance):
+ warnings.warn(
+ f"Azimuth time spacing is not regular. "
+ f"There are intervals which deviate more than "
+ f"{tolerance * 100:.3f}% from the median step "
+ f"({nominal_dt / np.timedelta64(1, 's'):.9f} s). "
+ f"Maximum deviation: "
+ f"{relative_deviation.max() * 100:.3f}%."
+ )
+ out = out.assign_coords(
+ azimuth_time=out.azimuth_time.values[0]
+ + np.arange(out.sizes["azimuth_time"]) * nominal_dt
)
+ return out
+
+ @staticmethod
+ def _align_azimuth(
+ dss: list[xr.Dataset], tolerance: float = 0.01
+ ) -> list[xr.Dataset]:
+ """Align all swaths to a shared azimuth-time axis."""
+ tol = 0.1 * np.diff(dss[0].azimuth_time.values[:2])[0]
+
+ start = max(ds.azimuth_time.values[0] for ds in dss)
+
+ stop = min(ds.azimuth_time.values[-1] for ds in dss)
+
+ dss_align = [ds.sel(azimuth_time=slice(start - tol, stop + tol)) for ds in dss]
+
+ # check number of azimuth lines
+ sizes = [ds.sizes["azimuth_time"] for ds in dss_align]
+ if len(np.unique(sizes)) != 1:
+ warnings.warn(
+ f"Aligned swaths have different azimuth sizes: {sizes}. "
+ "TOPSAR merge may require additional trimming."
+ )
+ min_size = min(sizes)
+ dss_align = [ds.isel(azimuth_time=slice(0, min_size)) for ds in dss_align]
+
+ # check azimuth spacing
+ for i, ds in enumerate(dss_align):
+ dt = ds.azimuth_time.diff("azimuth_time").values
+ nominal_dt = np.median(dt)
+ deviation = np.abs(dt - nominal_dt) / nominal_dt
+
+ if np.any(deviation > tolerance):
+ warnings.warn(
+ f"Azimuth spacing is irregular in swath {i}. "
+ f"Maximum deviation: {deviation.max():.3%} "
+ f"(threshold {tolerance:.3%})."
+ )
+
+ azimuth_time = dss_align[0].azimuth_time
+
+ return [ds.assign_coords(azimuth_time=azimuth_time) for ds in dss_align]
+
+ @staticmethod
+ def _merge_swaths(dss: list[xr.Dataset], tolerance: float = 0.01) -> xr.Dataset:
+ """Merge swaths along slant range."""
+ idxs = np.zeros((len(dss), 2), dtype=int)
+ step_tol = tolerance * np.diff(dss[0].slant_range_time.values[:2])[0]
+
+ for i, (ds0, ds1) in enumerate(zip(dss[:-1], dss[1:])):
+
+ r0 = ds0.slant_range_time.values
+ r1 = ds1.slant_range_time.values
+
+ # overlap region
+ mask0 = r0 >= r1[0] - step_tol
+ mask1 = r1 <= r0[-1] + step_tol
+
+ n_overlap = min(mask0.sum(), mask1.sum())
+
+ if n_overlap == 0:
+ raise ValueError(f"No overlap found between swaths {i} and {i+1}")
+
+ keep0 = n_overlap // 2
+ keep1 = n_overlap - keep0
+
+ idxs[i, 1] = len(r0) - keep0
+ idxs[i + 1, 0] = keep1
+
+ idxs[-1, 1] = dss[-1].sizes["slant_range_time"]
+
+ dss_cut = [
+ ds.isel(slant_range_time=slice(start, stop))
+ for ds, (start, stop) in zip(dss, idxs)
+ ]
+
+ out = xr.concat(
+ dss_cut,
+ dim="slant_range_time",
+ coords="minimal",
+ )
+
+ dr = out.slant_range_time.diff("slant_range_time").values
+ nominal_dr = np.median(dr)
+ deviation = np.abs(dr - nominal_dr) / nominal_dr
+ if np.any(deviation > tolerance):
+ warnings.warn(
+ f"Slant range spacing is irregular. "
+ f"Maximum deviation: {deviation.max():.3%} "
+ f"(threshold {tolerance:.3%})."
+ )
+ return out
+
class Sen1OCN(Sen1):
product_type = "OCN"
@@ -301,7 +615,14 @@ class Sen1OCN(Sen1):
resolution = kwargs.get("resolution")
if resolution is not None:
- assert_arg_is_instance(resolution, "resolution", (float, int))
+ assert_arg_is_instance(resolution, "resolution", (float, int, tuple))
+ if isinstance(resolution, tuple):
+ assert_arg_has_length(resolution, "resolution", 2)
+ if not all(isinstance(v, (float, int)) for v in resolution):
+ raise TypeError(
+ "resolution argument must contain exactly "
+ "two float or int values."
+ )
params.update(resolution=resolution)
bbox = kwargs.get("bbox")
@@ -336,7 +657,7 @@ class Sen1OCN(Sen1):
datatree: xr.DataTree,
includes: str | Iterable[str] | None = None,
excludes: str | Iterable[str] | None = None,
- resolution: float = None,
+ resolution: FloatInt | tuple[FloatInt, FloatInt] | None = None,
bbox: Sequence[float | int] | None = None,
crs: pyproj.CRS | None = None,
interp_methods: SpatialInterpMethods | None = None,
@@ -351,20 +672,20 @@ class Sen1OCN(Sen1):
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,
@@ -396,12 +717,12 @@ class Sen1OCN(Sen1):
bbox = source_gm.xy_bbox
if resolution is None:
if crs and not crs.is_geographic:
- center_lat = (
+ ref_point = (
(source_gm.xy_bbox[0] + source_gm.xy_bbox[2]) / 2,
(source_gm.xy_bbox[1] + source_gm.xy_bbox[3]) / 2,
)
resolution = transform_resolution(
- center_lat, source_gm.xy_res, source_gm.crs, crs
+ ref_point, source_gm.xy_res, source_gm.crs, crs
)
else:
resolution = source_gm.xy_res
@@ -423,14 +744,25 @@ class Sen1OCN(Sen1):
def register(registry: AnalysisModeRegistry):
- """Register Sentinel-1 analysis modes."""
registry.register(Sen1GRD)
+ registry.register(Sen1SLC)
registry.register(Sen1OCN)
+def _register_cache_uri(cache_uri: str) -> None:
+ _REGISTERED_CACHE_URIS.append(cache_uri)
+
+
+def _cleanup_registered_cache_uris() -> None:
+ for cache_uri in _REGISTERED_CACHE_URIS:
+ fs, path = fsspec.url_to_fs(cache_uri)
+ if fs.exists(path):
+ fs.rm(path, recursive=True)
+
+
def get_dem(
bbox: Sequence[float | int],
- resolution: float | None = None,
+ resolution: FloatInt | tuple[FloatInt, FloatInt] | None = None,
crs: pyproj.CRS | None = None,
):
"""Fetch and prepare a DEM for the given area of interest.
@@ -446,6 +778,7 @@ def get_dem(
Raises:
ValueError: If required credentials are missing or resolution is invalid.
"""
+
# check that environment variables are set
missing = [
name
@@ -473,8 +806,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),
@@ -494,64 +826,105 @@ def get_dem(
lon=slice(bbox_wgs84[0], bbox_wgs84[2]),
).chunk(_DEM_CHUNKSIZE)
else:
- if resolution is None:
- raise ValueError("Resolution must be provided if CRS is not None.")
+ dem = dem.to_dataset(name="dem")
if crs is None:
crs = _CRS_WGS84
+ if resolution is None:
+ source_gm = GridMapping.from_dataset(dem)
+ ref_point = (
+ (bbox_wgs84[0] + bbox_wgs84[2]) / 2,
+ (bbox_wgs84[1] + bbox_wgs84[3]) / 2,
+ )
+ resolution = transform_resolution(
+ ref_point, source_gm.xy_res, source_gm.crs, crs
+ )
target_gm = GridMapping.regular_from_bbox(
bbox,
resolution,
crs,
tile_size=(_DEM_CHUNKSIZE["lat"], _DEM_CHUNKSIZE["lon"]),
)
- dem = resample_in_space(dem.to_dataset(name="dem"), target_gm=target_gm).dem
+ dem = resample_in_space(dem, target_gm=target_gm).dem
return dem
-def convert_dem_to_ecef(dem: xr.DataArray, gm_dem: GridMapping) -> xr.DataArray:
+def _get_grid_parameters(
+ datatree: xr.DataTree,
+ footprint_scale_factor: tuple[float, float],
+ range_coord: str = "ground_range",
+) -> GridParams:
+ """Build grid parameters for RTC from Sentinel-1 metadata.
+
+ Args:
+ datatree: Source data tree.
+ footprint_scale_factor: Scaling for SAR footprint spacing.
+
+ Returns:
+ Grid parameters for terrain correction.
+ """
+
+ vh_group = [x for x in datatree.children if "VH" in x][0]
+ attrs = datatree[f"{vh_group}"].attrs["other_metadata"]["image_annotation"][
+ "image_information"
+ ]
+
+ az_scale, range_scale = footprint_scale_factor
+ if range_coord == "ground_range":
+ range0 = 0.0
+ d_range = attrs["range_pixel_spacing"]
+ spacing_range = attrs["range_pixel_spacing"]
+ else:
+ meas = datatree[f"{vh_group}/measurements"]
+ range0 = meas[range_coord].values[0]
+ d_range = meas[range_coord].values[1] - meas[range_coord].values[0]
+ spacing_range = 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(
+ range0=range0,
+ d_range=d_range,
+ spacing_range=spacing_range,
+ d_range_scale=d_range * range_scale,
+ range0_scale=(range0 - (0.5 * d_range) + (0.5 * d_range * range_scale)),
+ spacing_range_scale=spacing_range * range_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,
+ )
+
+
+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:
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.
"""
- x_dim, y_dim = gm_dem.xy_var_names
- transformer = pyproj.Transformer.from_crs(gm_dem.crs, _CRS_ECEF, always_xy=True)
+ x_dim, y_dim = gm_dem_params["xy_var_names"]
+ xx, yy = np.meshgrid(dem[x_dim].values, dem[y_dim].values, indexing="xy")
- 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",
- )
-
- 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,21 +1126,31 @@ 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,
+ range_coord: str = "ground_range",
+ 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:
- 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.
+ range_coord: Name of the coordinate representing the across-track axis.
+ 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.
@@ -775,11 +1158,21 @@ def backward_geocode(
t_shift: Time shift for the secant method.
Returns:
- Orbit time, distance vector, and velocity vector.
+ A dataset containing the optimized ground_range and azimuth time
+ for each target pixel and optionally the gamma area needed for RTC.
Raises:
ValueError: If the method is not supported.
"""
+ assert pos_coeff is not None
+ assert vel_coeff is not None
+ if range_coord == "ground_range":
+ 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,58 +1190,44 @@ 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"])
+ range_values = (
+ get_ground_range(gr_coeff, azimuth_time, time_slr)
+ if range_coord == "ground_range"
+ else time_slr
)
-
+ out = xr.Dataset({"azimuth_time": azimuth_time, range_coord: range_values})
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:
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.
"""
- 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 +1249,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 +1266,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,19 +1283,21 @@ 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.
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:
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,102 +1305,107 @@ 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:
- 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.
"""
- 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:
- 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.
"""
- 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,
+ range_coord: str = "ground_range",
) -> xr.DataArray:
"""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.
+ range_coord: Name of the coordinate representing the across-track axis.
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_scale) / _ONE_SECOND / params.d_az_scale
+ )
+ src_loc["gr_idx"] = (
+ src_loc[range_coord] - params.range0_scale
+ ) / params.d_range_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.spacing_range_scale * params.spacing_az_scale)
def fit_ground_range(time_slr_gcp: xr.DataArray, deg: int = 8) -> xr.DataArray:
@@ -1045,132 +1425,230 @@ 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"].data):
+ 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"],
+ range_coord: str = "ground_range",
) -> xr.Dataset:
"""Geocode data from SAR grid to map coordinates.
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.
+ src_loc: Source location dataset with target coordinates.
+ grid_params: Grid parameters for index conversion.
interp_method: Interpolation method.
+ range_coord: Name of the coordinate representing the across-track axis.
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[range_coord] - grid_params.range0) / grid_params.d_range
+ scr_indexing = _compute_indexing(data, az_idx, gr_idx, range_coord=range_coord)
+ 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_array.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,
- time_slr_gcp: xr.DataArray,
- sat_position: xr.DataArray,
+def get_source_location(
dem: xr.DataArray,
- apply_rtc: bool = True,
- grid_params: GridParams | None = None,
- interp_method: Literal["nearest", "bilinear"] = "nearest",
+ time_slr_gcp: xr.DataArray | None,
+ sat_position: xr.DataArray,
+ grid_params: GridParams,
+ gm_dem_grid: GridMapping,
+ apply_rtc: bool,
+ range_coord: str = "ground_range",
) -> 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.
+ gr_coeff = fit_ground_range(time_slr_gcp) if range_coord == "ground_range" else None
- 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)
-
- 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")
+ azimuth_data = xr.zeros_like(dem, dtype="datetime64[ns]")
+ template = xr.Dataset(
+ {"azimuth_time": azimuth_data, range_coord: 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
+ if "spatial_ref" in template:
+ template = template.drop_vars("spatial_ref")
+ gm_dem_params = {
+ "crs": gm_dem_grid.crs.to_wkt(),
+ "xy_var_names": gm_dem_grid.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,
+ "range_coord": range_coord,
+ "apply_rtc": apply_rtc,
+ "gm_dem_params": gm_dem_params,
+ },
+ template=template,
+ )
+ out.coords["spatial_ref"] = xr.DataArray(0, attrs=gm_dem_grid.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
+
+
+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,
+ range_coord: str = "ground_range",
+) -> 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[range_coord]))),
+ )
+ 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,
+ x_idx: np.ndarray,
+ y_idx: np.ndarray,
+ interp_method: Literal["nearest", "bilinear"] | None = None,
+) -> np.ndarray:
+ """Sample a 3D array at fractional indices (y_idx, x_idx)."""
+ if interp_method == "nearest":
+ 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)
+ 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}'."
+ )
=====================================
xarray_eopf/amodes/sentinel3.py
=====================================
@@ -56,7 +56,14 @@ class Sen3(AnalysisMode, ABC):
resolution = kwargs.get("resolution")
if resolution is not None:
- assert_arg_is_instance(resolution, "resolution", (int, float))
+ assert_arg_is_instance(resolution, "resolution", (float, int, tuple))
+ if isinstance(resolution, tuple):
+ assert_arg_has_length(resolution, "resolution", 2)
+ if not all(isinstance(v, (float, int)) for v in resolution):
+ raise TypeError(
+ "resolution argument must contain exactly "
+ "two float or int values."
+ )
params.update(resolution=resolution)
bbox = kwargs.get("bbox")
=====================================
xarray_eopf/backend.py
=====================================
@@ -19,6 +19,7 @@ from .constants import (
OP_MODE_ANALYSIS,
OP_MODE_NATIVE,
OP_MODES,
+ FloatInt,
OpMode,
Sen1InterpMethods,
)
@@ -122,14 +123,16 @@ class EopfBackend(BackendEntrypoint):
variables: str | Iterable[str] | None = None,
# params for op_mode=analysis
product_type: str | None = None,
- resolution: int | float | None = None,
+ resolution: FloatInt | tuple[FloatInt, FloatInt] | None = None,
bbox: Sequence[int | float] | None = None,
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
@@ -154,9 +157,11 @@ class EopfBackend(BackendEntrypoint):
group_sep: Separator string used to concatenate groups names
to create prefixes for unique variable and dimension names.
Defaults to the underscore character (`"_"`)
- resolution: Target resolution for all spatial
- data variables / bands. For Sentinel-2 products it be one of
- `10`, `20`, or `60`. Only used if `op_mode="analysis"`.
+ resolution: Target spatial resolution for all spatial
+ data variables/bands. The resolution can be specified as a float, integer,
+ or a tuple in the form `(easting, northing)`. For Sentinel-2 products,
+ valid resolutions are `10`, `20`, and `60` meters. This parameter is only
+ used when `op_mode="analysis"`.
bbox: Bounding box [west, south, east, north], used for subsetting.
crs: coordinate reference system of output dataset. Can be provided as a
`str` or a `pyproj.CRS` object. If a string is given, it will be parsed
@@ -207,6 +212,11 @@ class EopfBackend(BackendEntrypoint):
when not provided.
apply_rtc: Whether to apply radiometric terrain correction (RTC) for
Sentinel-1 analysis mode. Defaults to `True`.
+ 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
@@ -266,6 +276,7 @@ class EopfBackend(BackendEntrypoint):
dem=dem,
footprint_scale_factor=footprint_scale_factor,
apply_rtc=apply_rtc,
+ cache_uri=cache_uri,
)
dataset = analysis_mode.convert_datatree(
datatree, includes=variables, **params
=====================================
xarray_eopf/version.py
=====================================
@@ -2,4 +2,4 @@
# Permissions are hereby granted under the terms of the Apache 2.0 License:
# https://opensource.org/license/apache-2-0.
-version = "0.2.9"
+version = "0.3.0"
View it on GitLab: https://salsa.debian.org/debian-gis-team/xarray-eopf/-/commit/df45a2b43a86415bcd338a49fc55028c8e5a053d
--
View it on GitLab: https://salsa.debian.org/debian-gis-team/xarray-eopf/-/commit/df45a2b43a86415bcd338a49fc55028c8e5a053d
You're receiving this email because of your account on salsa.debian.org. Manage all notifications: https://salsa.debian.org/-/profile/notifications | Help: https://salsa.debian.org/help
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://alioth-lists.debian.net/pipermail/pkg-grass-devel/attachments/20260729/abd2e04a/attachment-0001.htm>
More information about the Pkg-grass-devel
mailing list