feat: implement envisat gtc (#44)

* feat: implement envisat gtc

* ci: update tox dependencies

* tests: fix test_xarray_backends.py

* Added direct parsing of Envisat files.

This corrects slant range mistake from previous version for terrain correction and it is required for further development.

* fix: small import error

* style: code style improvements

* style: code style improvements

* style: code style improvements

* style: code style improvements

* tests: fix tests

* Fix typing

* A bit of code style changes

* Fix sonarqube issues

---------

Co-authored-by: Anton Perepelenko <Anton.Perepelenko@cgi.com>
Co-authored-by: Priit Pender <priit@alloca.ee>

Author: Achaad

asar_xarray/asar.py

@@ -3,20 +3,21 @@
 import os
 from typing import Dict, Any
 
+import pandas as pd
 import xarray as xr
 import numpy as np
 from osgeo import gdal
 from xarray.backends import AbstractDataStore
 from xarray.core.types import ReadBuffer
-from asar_xarray import reader, utils
+from asar_xarray import reader, utils, envisat_direct
 from loguru import logger
 
 from asar_xarray.derived_subdatasets_metadata import process_derived_subdatasets_metadata
 from asar_xarray.general_metadata import process_general_metadata
 from asar_xarray.records_metadata import process_records_metadata
 
 
-def get_attributes(gdal_dataset: gdal.Dataset) -> Dict[str, Any]:
+def get_metadata(gdal_dataset: gdal.Dataset) -> Dict[str, Any]:
     """
     Build xarray attributes from gdal dataset to be used in xarray.
 
@@ -54,24 +55,110 @@ def open_asar_dataset(filepath: str | os.PathLike[Any] | ReadBuffer[Any] | Abstr
     gdal_dataset: gdal.Dataset = reader.get_gdal_dataset(filepath)
 
     # Extract metadata attributes
-    attributes = get_attributes(gdal_dataset)
+    metadata = get_metadata(gdal_dataset)
 
-    # Read pixel data from the dataset
-    data = gdal_dataset.ReadAsArray()
+    # Duplicate, read directly from file, as gdal does not parse some necessary metadata
+    metadata["direct_parse"] = envisat_direct.parse_direct(filepath)
 
     # Create an xarray Dataset with pixel data and metadata attributes
-    dataset: xr.Dataset = xr.Dataset(
-        data_vars={'pixel_values': (('y', 'x'), data)},
-        coords={
-            'x': np.arange(data.shape[1]),
-            'y': np.arange(data.shape[0])
-        },
-        attrs=attributes
-    )
+    dataset: xr.Dataset = create_dataset(metadata, filepath)
 
     return dataset
 
 
+def create_dataset(metadata: dict[str, Any], filepath: str) -> xr.Dataset:
+    """
+    Create an xarray Dataset from ASAR metadata and file path.
+
+    This function constructs the coordinates, attributes, and data variables
+    for an ASAR product, using the provided metadata and file path.
+
+    :param metadata: Dictionary containing ASAR product metadata.
+    :param filepath: Path to the ASAR dataset file.
+    :return: An xarray Dataset with pixel data, coordinates, and attributes.
+    """
+    number_of_samples = metadata["line_length"]
+    product_first_line_utc_time = metadata["first_line_time"]
+    product_last_line_utc_time = metadata["last_line_time"]
+    print(product_first_line_utc_time)
+
+    number_of_lines = metadata["records"]["main_processing_params"]["num_output_lines"]
+    azimuth_time_interval = 1 / metadata["records"]["main_processing_params"]["image_parameters"]["prf_value"][0]
+    range_sampling_rate = metadata["records"]["main_processing_params"]["range_samp_rate"]
+    image_slant_range_time = metadata["direct_parse"]["slant_time_first"] * 1e-9
+
+    number_of_bursts = 0
+
+    attrs = {
+        "family_name": "Envisat",
+        "number": 1,
+        "mode": 1,
+        "swaths": metadata["swath"],
+        "orbit_number": metadata["abs_orbit"],
+        "relative_orbit_number": metadata["rel_orbit"],
+        "pass": metadata["pass"],
+        "transmitter_receiver_polarisations": metadata["mds1_tx_rx_polar"],
+        "product_type": "SLC",
+        "start_time": product_first_line_utc_time,
+        "stop_time": product_last_line_utc_time,
+
+        "radar_frequency": metadata["records"]["main_processing_params"]["radar_freq"] / 1e9,
+        "ascending_node_time": "",
+        "azimuth_pixel_spacing": metadata["records"]["main_processing_params"]["azimuth_spacing"],
+        "range_pixel_spacing": metadata["records"]["main_processing_params"]["range_samp_rate"],
+        "product_first_line_utc_time": product_first_line_utc_time,
+        "product_last_line_utc_time": product_last_line_utc_time,
+        "azimuth_time_interval": azimuth_time_interval,
+        "image_slant_range_time": image_slant_range_time,
+        "range_sampling_rate": range_sampling_rate,
+        "incidence_angle_mid_swath": metadata["direct_parse"]["incidence_angle_center"],
+        "metadata": metadata
+    }
+
+    azimuth_time = compute_azimuth_time(
+        product_first_line_utc_time, product_last_line_utc_time, number_of_lines
+    )
+
+    if number_of_bursts == 0:
+        swap_dims = {"line": "azimuth_time", "pixel": "slant_range_time"}
+    else:
+        raise NotImplementedError(
+            "Burst processing is not implemented yet."
+        )
+
+    coords: dict[str, Any] = {
+        "pixel": np.arange(0, number_of_samples, dtype=int),
+        "line": np.arange(0, number_of_lines, dtype=int),
+        # set "units" explicitly as CF conventions don't support "nanoseconds".
+        # See: https://github.com/pydata/xarray/issues/4183#issuecomment-685200043
+        "azimuth_time": (
+            "line",
+            azimuth_time,
+            {},
+            {"units": f"microseconds since {azimuth_time[0]}"},
+        ),
+    }
+
+    slant_range_time = np.linspace(
+        image_slant_range_time,
+        image_slant_range_time + (number_of_samples - 1) / range_sampling_rate,
+        number_of_samples,
+    )
+    coords["slant_range_time"] = ("pixel", slant_range_time)
+
+    data = xr.open_dataarray(filepath, engine='rasterio')
+    data.encoding.clear()
+    data = data.squeeze("band").drop_vars(["band", "spatial_ref"])
+    data = data.rename({"y": "line", "x": "pixel"})
+    data = data.assign_coords(coords)
+    data = data.swap_dims(swap_dims)
+
+    data.attrs.update(attrs)
+    data.encoding.update({})
+
+    return xr.Dataset(attrs=attrs, data_vars={"measurements": data})
+
+
 def get_chirp_parameters(dataset: gdal.Dataset) -> dict[str, Any]:
     """
     Parse dataset metadata for chirp parameters.
@@ -134,3 +221,27 @@ def get_chirp_cal_pulse_info(metadata: dict[str, str]) -> list[dict[str, Any]]:
 
     # Convert the dictionary of dictionaries to the list of dictionaries
     return list(cal_info_dict.values())
+
+
+def compute_azimuth_time(
+        product_first_line_utc_time: np.datetime64,
+        product_last_line_utc_time: np.datetime64,
+        number_of_lines: int
+) -> np.ndarray:
+    """
+    Compute an array of azimuth times for each line in the ASAR product.
+
+    This function generates a sequence of evenly spaced time values between the
+    first and last line UTC times, corresponding to the number of lines in the product.
+
+    :param product_first_line_utc_time: UTC time of the first line (as np.datetime64).
+    :param product_last_line_utc_time: UTC time of the last line (as np.datetime64).
+    :param number_of_lines: Total number of lines in the product.
+    :return: Numpy array of azimuth times for each line.
+    """
+    azimuth_time = pd.date_range(
+        start=product_first_line_utc_time,
+        end=product_last_line_utc_time,
+        periods=number_of_lines
+    )
+    return np.asarray(azimuth_time.values, dtype='datetime64[ns]')

asar_xarray/envisat_direct.py

@@ -0,0 +1,92 @@
+"""Module for parsing Envisat direct access data structures and extracting metadata."""
+import struct
+from typing import Any
+
+
+def parse_int(s: str) -> int:
+    """Parse an integer value from a string representation of a field."""
+    s = s.replace("<bytes>", "")
+    s = s[s.index("=") + 1:]
+    # print(s)
+    return int(s)
+
+
+class EnvisatADS:
+    """Class representing an Envisat Annotation Data Set (ADS) descriptor."""
+
+    def __init__(self, buffer: bytes) -> None:
+        """Initialize the ADS descriptor from a buffer of bytes."""
+        str_arr = buffer.decode("ascii").split("\n")
+        name = str_arr[0].replace("DS_NAME=\"", "").replace("\"", "").strip()
+        self.name = name
+        self.num = parse_int(str_arr[5])
+        self.size = parse_int(str_arr[4])
+        self.offset = parse_int(str_arr[3])
+
+    def __str__(self) -> str:
+        """Return a string representation of the ADS descriptor."""
+        return "Envisat ADS: \"{}\" {} {} {}".format(self.name, self.offset, self.size, self.num)
+
+
+def parse_direct(path: str) -> dict[str, Any]:
+    """
+    Parse an Envisat product file and extract relevant metadata fields.
+
+    Args:
+    ----
+        path (str): Path to the Envisat product file.
+
+    returns: Dictionary containing extracted metadata fields.
+    """
+    metadata = {}
+    file_buffer = None
+    with open(path, "rb") as fp:
+        file_buffer = fp.read()
+
+    # read main product header and confirm sph size location
+    mph_size = 1247
+    mph_str = file_buffer[0:mph_size].decode("ascii")
+    assert (mph_str.index("SPH_SIZE") == 1104)
+    sph_size = parse_int(mph_str[1104:1104 + 20])
+
+    sph_buf = file_buffer[mph_size:mph_size + sph_size]
+    dsd_size = 280
+    dsd_num = 18
+    dsd_buf = sph_buf[sph_size - dsd_size * dsd_num:]
+
+    for i in range(dsd_num):
+        ads = EnvisatADS(dsd_buf[i * dsd_size:(i + 1) * dsd_size])
+        if ads.name == "GEOLOCATION GRID ADS":
+            rec_size = 521
+            assert ((ads.size // ads.num) == rec_size)
+            geoloc_buf = file_buffer[ads.offset:ads.offset + ads.size]
+            middle_idx = ads.num // 2
+            geoloc_record = geoloc_buf[middle_idx * rec_size: (middle_idx + 1) * rec_size]
+            # Geolocation Grid ADSRs header
+            header_size = 12 + 1 + 4 + 4 + 4
+            # tiepoints, 11 of big endian floats for each of the following:
+            # samp numbers, slant range times, angles, lats, longs
+            block_size = 11 * 4
+            slant_time_offset = header_size + 1 * block_size
+            incidence_angle_offset = header_size + 2 * block_size
+            # adjust to middle of 11
+            incidence_angle_offset += 5 * 4
+
+            slant_time_first = struct.unpack(">f", geoloc_record[slant_time_offset:slant_time_offset + 4])[0]
+            incidence_angle_middle = \
+                struct.unpack(">f", geoloc_record[incidence_angle_offset:incidence_angle_offset + 4])[0]
+
+            metadata["slant_time_first"] = slant_time_first
+            metadata["incidence_angle_center"] = incidence_angle_middle
+
+        if ads.name == "MAIN PROCESSING PARAMS ADS":
+
+            main_processing_params_buf = file_buffer[ads.offset:ads.offset + ads.size]
+            if len(main_processing_params_buf) == 10069:
+                sigma_buf = main_processing_params_buf[2029:2029 + 4020]
+                gammma_buf = main_processing_params_buf[2029 + 4020:]
+
+                metadata["sigma_calib_vector"] = sigma_buf
+                metadata["gamma_calib_vector"] = gammma_buf
+
+    return metadata

asar_xarray/main.py

@@ -7,4 +7,4 @@
     test_file = ('../tests/resources'
                  '/ASA_IMS_1PNESA20040109_194924_000000182023_00157_09730_0000.N1')
     dataset = xr.open_dataset(test_file, engine='asar')
-    print(dataset)
+    print(dataset.attrs)

environment.yml

@@ -14,5 +14,8 @@ dependencies:
   - lxml>=5.3.0
   - mypy>=1.15.0
   - loguru>=0.7.0
+  - pandas-stubs >=2.0.0
   - pip:
     - pytest-stub>=1.1.0
+    - rasterio>=1.4.0
+    - rioxarray>=0.19.0

pyproject.toml

@@ -17,7 +17,7 @@ classifiers = [
     'Topic :: Scientific/Engineering :: GIS'
 ]
 readme = 'README.md'
-license = { text = 'Apache-2.0' }
+license-files = ['LICENSE']
 keywords = ['ASAR', 'ENVISAT', 'xarray', 'earth-observation', 'remote-sensing', 'radar', 'satellite-imagery', 'sar',
     'synthetic-aperture-radar']
 requires-python = '>=3.12'

tests/test_asar.py

@@ -7,6 +7,6 @@ def test_open_asar_dataset(resources_dir: str) -> None:
     filepath = os.path.join(resources_dir, 'ASA_IMS_1PNESA20040109_194924_000000182023_00157_09730_0000.N1')
     dataset = asar.open_asar_dataset(filepath)
     assert dataset is not None
-    assert dataset.dims['x'] == 5174
-    assert dataset.dims['y'] == 30181
-    assert dataset.data_vars['pixel_values'].data.shape == (30181, 5174)
+    assert dataset.dims['slant_range_time'] == 5174
+    assert dataset.dims['azimuth_time'] == 30181
+    assert dataset.data_vars['measurements'].data.shape == (30181, 5174)

tests/test_xaray_backends.py

@@ -9,6 +9,6 @@ def test_open_dataset(resources_dir: str) -> None:
     file = os.path.join(resources_dir, 'ASA_IMS_1PNESA20040109_194924_000000182023_00157_09730_0000.N1')
     dataset = xr.open_dataset(file, engine='asar')
     assert dataset is not None
-    assert dataset.dims['x'] == 5174
-    assert dataset.dims['y'] == 30181
-    assert dataset.data_vars['pixel_values'].data.shape == (30181, 5174)
+    assert dataset.dims['slant_range_time'] == 5174
+    assert dataset.dims['azimuth_time'] == 30181
+    assert dataset.data_vars['measurements'].data.shape == (30181, 5174)

tox.ini

@@ -19,6 +19,8 @@ deps =
     pytest
     pytest-cov
     mypy
+    rasterio>=1.4.0
+    rioxarray>=0.19.0
 commands =
     coverage run -m pytest -s
     coverage xml