Skip to content

Commit

Permalink
Merge remote-tracking branch 'origin/main' into declare-typed
Browse files Browse the repository at this point in the history
  • Loading branch information
alexamici committed Apr 25, 2022
2 parents 5fa8bb4 + 09af28d commit 51add7c
Show file tree
Hide file tree
Showing 3 changed files with 62 additions and 21 deletions.
22 changes: 15 additions & 7 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -58,6 +58,11 @@ Due to the inherent complexity and redundancy of the SAFE format *xarray-sentine
maps it to a tree of *groups* where every *group* may be opened as a `Dataset`,
but it may also contain *subgroups*, that are listed in the `subgroups` attribute.

The following sections show some example of xarray-sentinel usage.
In the `notebooks` folder you
can also find notebooks, one for each supported product, that allow you to explore the
data in more detail using the xarray-sentinel functions.

### The root dataset

For example let's explore the Sentinel-1 SLC Stripmap product in the local folder
Expand Down Expand Up @@ -137,6 +142,9 @@ associated with the image line
and `slant_range_time` is an `np.float64` coordinate that contains the two-way range time interval
in seconds associated with the image pixel.

Since Sentinel-1 IPF version 3.40, a unique identifier for bursts has been added to the SLC product metadata.
For these products, the list of the burst ids is stored the `burst_ids` dataset attribute.

### Metadata datasets

The measurement group contains several subgroups with metadata associated with the image. Currently,
Expand Down Expand Up @@ -230,14 +238,14 @@ but also because the measurement array is a collage of sub-images called *bursts

*xarray-sentinel* provides a helper function that crops a burst out of a measurement dataset for you.

You need to first open the desired measurement dataset, for example, the VH polarisation
of the first IW swath of the `S1B_IW_SLC__1SDV_20210401T052622_20210401T052650_026269_032297_EFA4`
You need to first open the desired measurement dataset, for example, the HH polarisation
of the first IW swath of the `S1A_IW_SLC__1SDH_20220414T102209_20220414T102236_042768_051AA4_E677.SAFE`
product, in the current folder:

```python-repl
>>> slc_iw_v340_path = "tests/data/S1A_IW_SLC__1SDH_20220414T102209_20220414T102236_042768_051AA4_E677.SAFE"
>>> slc_iw1_hh = xr.open_dataset(slc_iw_v340_path, group="IW1/HH", engine="sentinel-1")
>>> slc_iw1_hh
>>> slc_iw1_v340_hh = xr.open_dataset(slc_iw_v340_path, group="IW1/HH", engine="sentinel-1")
>>> slc_iw1_v340_hh
<xarray.Dataset>
Dimensions: (pixel: 21169, line: 13500)
Coordinates:
Expand Down Expand Up @@ -271,7 +279,7 @@ Now the 9th burst out of 9 can be cropped from the swath data using `burst_index

```python-repl
>>> import xarray_sentinel
>>> xarray_sentinel.crop_burst_dataset(slc_iw1_hh, burst_index=8)
>>> xarray_sentinel.crop_burst_dataset(slc_iw1_v340_hh, burst_index=8)
<xarray.Dataset>
Dimensions: (slant_range_time: 21169, azimuth_time: 1500)
Coordinates:
Expand All @@ -298,11 +306,11 @@ Attributes: ...
```

For products processed with processor versions 3.40 or higher, it is also possible to select the burst
If IPF processor version is 3.40 or higher, it is also possible to select the burst
to be cropped using the `burst_id` key:

```python-repl
>>> xarray_sentinel.crop_burst_dataset(slc_iw1_hh, burst_id=365923)
>>> xarray_sentinel.crop_burst_dataset(slc_iw1_v340_hh, burst_id=365923)
<xarray.Dataset>
Dimensions: (slant_range_time: 21169, azimuth_time: 1500)
Coordinates:
Expand Down
6 changes: 3 additions & 3 deletions tests/test_30_xarray_backends.py
Original file line number Diff line number Diff line change
Expand Up @@ -199,15 +199,15 @@ def test_burst_id_attribute() -> None:
)

res = xr.open_dataset(product_path, engine="sentinel-1", group="IW1/HH") # type: ignore
assert "bursts_ids" in res.attrs
assert len(res.attrs["bursts_ids"]) == res.attrs["number_of_bursts"]
assert "burst_ids" in res.attrs
assert len(res.attrs["burst_ids"]) == res.attrs["number_of_bursts"]

product_path = (
DATA_FOLDER
/ "S1B_IW_SLC__1SDV_20210401T052622_20210401T052650_026269_032297_EFA4.SAFE"
)
res = xr.open_dataset(product_path, engine="sentinel-1", group="IW1/VV") # type: ignore
assert "bursts_ids" not in res.attrs
assert "burst_ids" not in res.attrs


def test_open_pol_dataset_preferred_chunks() -> None:
Expand Down
55 changes: 44 additions & 11 deletions xarray_sentinel/sentinel1.py
Original file line number Diff line number Diff line change
Expand Up @@ -465,10 +465,10 @@ def open_pol_dataset(
swap_dims = {"line": "azimuth_time", "pixel": "slant_range_time"}
else:
if "burstId" in swath_timing["burstList"]["burst"][0]:
bursts_ids = []
burst_ids = []
for burst in swath_timing["burstList"]["burst"]:
bursts_ids.append(burst["burstId"]["$"])
attrs["bursts_ids"] = bursts_ids
burst_ids.append(burst["burstId"]["$"])
attrs["burst_ids"] = burst_ids
lines_per_burst = swath_timing["linesPerBurst"]
attrs.update(
{
Expand Down Expand Up @@ -582,6 +582,17 @@ def crop_burst_dataset(
use_center: bool = False,
burst_id: T.Optional[int] = None,
) -> DataArrayOrDataset:
"""
Returns the measurement dataset cropped to the selected burst.
Only one keyword between 'burst_index' and 'azimuth_anx_time' and 'burst_id' must be defined.
:param xr.Dataset pol_dataset: measurement dataset
:param int burst_index: burst index can take values from 1 to the number of bursts
:param float azimuth_anx_time: azimuth anx time of first line of the bursts
To use the center instead of the first line, set `use_center=True`
:param bool use_center: If `true`, it uses as reference the azimuth anx time of the burst center instead of the first line
:param int burst_id: for product processed with Sentinel-1 IPF version 3.40 or higher,
the burst can be selected using the relative burst id.
"""
burst_definitions = (
(burst_index is not None)
+ (azimuth_anx_time is not None)
Expand All @@ -598,16 +609,16 @@ def crop_burst_dataset(
pol_dataset, azimuth_anx_time, use_center=use_center
)
elif burst_id is not None:
bursts_ids = pol_dataset.attrs.get("bursts_ids")
if bursts_ids is None:
burst_ids = pol_dataset.attrs.get("burst_ids")
if burst_ids is None:
raise TypeError(
"'bursts_ids' list can't be found in product attributes, "
"'burst_ids' list can't be found in product attributes, "
"probably Sentinel-1 IPF processor version is older than 3.40"
)
try:
burst_index = bursts_ids.index(burst_id)
burst_index = burst_ids.index(burst_id)
except ValueError:
raise KeyError(f"{burst_id=} not found in product {bursts_ids=}")
raise KeyError(f"{burst_id=} not found in product {burst_ids=}")
else:
raise TypeError(
"one keyword between 'burst_index' and 'azimuth_anx_time' must be defined"
Expand All @@ -628,9 +639,9 @@ def crop_burst_dataset(
ds.attrs["azimuth_anx_time"] = burst_azimuth_anx_times.values[0] / ONE_SECOND
ds = ds.swap_dims({"line": "azimuth_time", "pixel": "slant_range_time"})
ds.attrs["burst_index"] = burst_index
if "bursts_ids" in ds.attrs:
ds.attrs["burst_id"] = ds.attrs["bursts_ids"][burst_index]
_ = ds.attrs.pop("bursts_ids")
if "burst_ids" in ds.attrs:
ds.attrs["burst_id"] = ds.attrs["burst_ids"][burst_index]
_ = ds.attrs.pop("burst_ids")
_ = ds.attrs.pop("subgroups", None)
return ds

Expand All @@ -648,6 +659,11 @@ def mosaic_slc_iw(
def calibrate_amplitude(
digital_number: xr.DataArray, calibration_lut: xr.DataArray
) -> xr.DataArray:
"""Returns the calibrated amplitude. The calibration is done using the calibration LUT in the product metadata.
:param digital_number: digital numbers to be calibrated
:param calibration_lut: calibration LUT (sigmaNought, betaNought or gamma).
The LUT can be opened using the measurement sub-group `calibration`
"""
calibration = calibration_lut.interp(
line=digital_number.line,
pixel=digital_number.pixel,
Expand All @@ -669,6 +685,14 @@ def calibrate_intensity(
as_db: bool = False,
min_db: T.Optional[float] = -40.0,
) -> xr.DataArray:
"""
Returns the calibrated intensity. The calibration is done using the calibration LUT in the product metadata.
:param digital_number: digital numbers to be calibrated
:param calibration_lut: calibration LUT (sigmaNought, betaNought or gamma).
The LUT can be opened using the measurement sub-group `calibration`.
:param as_db: if True, returns the data in db
:param min_db: minimal value in db, to avoid infinity values.
"""
amplitude = calibrate_amplitude(digital_number, calibration_lut)
intensity = abs(amplitude) ** 2
if as_db:
Expand All @@ -693,6 +717,15 @@ def slant_range_time_to_ground_range(
slant_range_time: xr.DataArray,
coordinate_conversion: xr.Dataset,
) -> xr.DataArray:
"""
Convert the slant range time coordinates to ground range coordinates using the coordinate conversion `sr0`
and `srgrCoefficients` product metadata
:param azimuth_time: azimuth time coordinates
:param slant_range_time: slant range time
:param coordinate_conversion: coordinate conversion dataset.
The coordinate conversion dataset can be opened using the measurement sub-groub `coordinate_conversion`
:return:
"""
slant_range = SPEED_OF_LIGHT / 2.0 * slant_range_time
cc = coordinate_conversion.interp(azimuth_time=azimuth_time)
x = slant_range - cc.sr0
Expand Down

0 comments on commit 51add7c

Please sign in to comment.