Added ERS1 & ERS2 support with the external file.

Calibration vector now properly uses the antenna gain information.

Author: priit111

asar_xarray/asar.py

@@ -35,8 +35,7 @@ def get_metadata(gdal_dataset: gdal.Dataset) -> Dict[str, Any]:
     return attributes
 
 
-def open_asar_dataset(
-        filepath: str | os.PathLike[Any] | ReadBuffer[Any] | bytes | memoryview | AbstractDataStore) -> xr.Dataset:
+def open_asar_dataset(filepath: str | os.PathLike[Any] | ReadBuffer[Any] | AbstractDataStore) -> xr.Dataset:
     """
     Open an ASAR dataset and converts it into an xarray Dataset.
 
@@ -107,8 +106,8 @@ def create_dataset(metadata: dict[str, Any], filepath: str) -> xr.Dataset:
         "product_type": product_type,
         "start_time": product_first_line_utc_time,
         "stop_time": product_last_line_utc_time,
-        "range_pixel_spacing": metadata["range_spacing"],
-        "azimuth_spacing": metadata["azimuth_spacing"],
+        "range_pixel_spacing" : metadata["range_spacing"],
+        "azimuth_pixel_spacing" : metadata["azimuth_spacing"],
         "radar_frequency": metadata["records"]["main_processing_params"]["radar_freq"] / 1e9,
         "ascending_node_time": "",
         "azimuth_pixel_spacing": metadata["records"]["main_processing_params"]["azimuth_spacing"],
@@ -117,14 +116,15 @@ def create_dataset(metadata: dict[str, Any], filepath: str) -> xr.Dataset:
         "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"] * 2 * math.pi / 360,
+        "incidence_angle_mid_swath": metadata["direct_parse"]["incidence_angle_center"] * 2 * math.pi / 360 ,
         "metadata": metadata
     }
 
     azimuth_time = compute_azimuth_time(
         product_first_line_utc_time, product_last_line_utc_time, number_of_lines
     )
 
+
     swap_dims = {"line": "azimuth_time", "pixel": "slant_range_time"}
 
     coords: dict[str, Any] = {
@@ -180,6 +180,7 @@ def get_chirp_parameters(dataset: gdal.Dataset) -> dict[str, Any]:
             new_key = key.replace('CHIRP_PARAMS_ADS_CHIRP_', '').lower()
             params['chirp'][new_key] = float(value)
 
+    #params['elev_corr_factor'] = float(metadata.get('CHIRP_PARAMS_ADS_ELEV_CORR_FACTOR'))
     params['cal_pulse_info'] = get_chirp_cal_pulse_info(metadata)
 
     return params

asar_xarray/aux/ERS1/ER1_XCA_AXVXXX20100209_000000_19910717_000000_19980224_000000

@@ -4,6 +4,7 @@
 import os
 import math
 import pathlib
+import numpy as np
 
 
 def parse_int(s: str) -> int:
@@ -14,6 +15,35 @@ def parse_int(s: str) -> int:
     return int(s)
 
 
+
+def calc_distance(latitude, longitude, altitude=0.0):
+    """Lat/Lon to distance from earth center"""
+    DTOR = math.pi / 180.0
+    RTOD = 180 / math.pi
+    A = 6378137.0
+    B = 6356752.3142451794975639665996337
+    FLAT_EARTH_COEF = 1.0 / ((A - B) / A)
+    E2 = 2.0 / FLAT_EARTH_COEF - 1.0 / (FLAT_EARTH_COEF * FLAT_EARTH_COEF)
+    EP2 = E2 / (1 - E2)
+
+    lat = latitude * DTOR
+    lon = longitude * DTOR
+
+    sin_lat = math.sin(lat)
+    cos_lat = math.cos(lat)
+
+    N = (A / math.sqrt(1.0 - E2 * sin_lat * sin_lat))
+    NcosLat = (N + altitude) * cos_lat
+
+    sin_lon = math.sin(lon)
+    cos_lon = math.cos(lon)
+
+    x_pos = NcosLat * cos_lon
+    y_pos = NcosLat * sin_lon
+    z_pos = (N + altitude - E2 * N) * sin_lat
+    return math.sqrt(x_pos**2 + y_pos**2 + z_pos**2)
+
+
 class EnvisatADS:
     """Class representing an Envisat Annotation Data Set (ADS) descriptor."""
 
@@ -32,7 +62,7 @@ def __str__(self) -> str:
         return "Envisat ADS: \"{}\" {} {} {}".format(self.name, self.offset, self.size, self.num)
 
 
-def parse_direct(path: str, gdal_metadata: dict[str, Any]) -> dict[str, Any]:
+def parse_direct(path: str, gdal_metadata) -> dict[str, Any]:
     """
     Parse an Envisat product file and extract relevant metadata fields.
 
@@ -42,7 +72,8 @@ def parse_direct(path: str, gdal_metadata: dict[str, Any]) -> dict[str, Any]:
 
     returns: Dictionary containing extracted metadata fields.
     """
-    metadata: dict[str, Any] = {}
+    metadata = {}
+    file_buffer = None
     with open(path, "rb") as fp:
         file_buffer = fp.read()
 
@@ -59,145 +90,161 @@ def parse_direct(path: str, gdal_metadata: dict[str, Any]) -> dict[str, Any]:
 
     for i in range(dsd_num):
         ads = EnvisatADS(dsd_buf[i * dsd_size:(i + 1) * dsd_size])
+        #print(ads.name)
         if ads.name == "GEOLOCATION GRID ADS":
-            extract_geolocation_metadata(ads, file_buffer, metadata)
+            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
+
+            lats_buffer = geoloc_buf[header_size + 11 * 4 * 3 : header_size + 11 * 4 * 4]
+            lons_buffer = geoloc_buf[header_size + 11 * 4 * 4: header_size + 11 * 4 * 5]
+
+
+            lats = list(struct.unpack(">11i", lats_buffer))
+            lons = list(struct.unpack(">11i", lons_buffer))
+            lats = [e * 1e-6 for e in lats]
+            lons = [e * 1e-6 for e in lons]
+
+
+            # 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
+
+
+
+        ext_cal_buf = None
         if ads.name == "EXTERNAL CALIBRATION":
-            extract_external_calibration_metadata(ads, gdal_metadata, metadata)
 
-        if ads.name == "SR GR ADS" and ads.size > 0:
-            srgr_coeffs = extract_srgr_coeffs(ads, file_buffer)
+            auxfolder = pathlib.Path(os.path.abspath(__file__)).parent
+            prod_name = gdal_metadata["product"]
+            auxfolder /= "aux"
+            if ".N1" in prod_name:
+                auxfolder /= "ASAR_Auxiliary_Files"
+                auxfolder /= "ASA_XCA_AX"
+            elif ".E1" in prod_name:
+                auxfolder /= "ERS1"
+            elif ".E2" in prod_name:
+                auxfolder /= "ERS2"
+
+            for p in os.scandir(auxfolder):
+                if ads.filename == p.name:
+                    with open(p.path, "rb") as fp:
+                        ext_cal_buf = fp.read()
+                        pol = gdal_metadata["mds1_tx_rx_polar"]
+                        swath = gdal_metadata["swath"]
+                        swath_offset = ord(swath[2]) - ord("1")
+                        pol_offset = {"H/H" : 0, "V/V" :1, "H/V" : 2, "V/H":3}[pol]
+
+                        # Envisat_Product_Spec_Vol8.pdf
+                        # 8.6.2 External Calibration Data
+                        offset = 1247 + 378
+
+                        offset += 12
+                        offset += 4
+                        offset += 26 * 28 + 4 * 4
+
+                        mid_angles = struct.unpack(">7f", ext_cal_buf[offset:offset+4*7])
+
+                        offset += 8 * 4
+
+                        offset += 804 * 4 * swath_offset
+                        offset += 201 * 4 * pol_offset
+
+                        antenna_gains = struct.unpack(">201f", ext_cal_buf[offset:offset + 4 * 201])
+                        metadata["antenna_ref_elev_angle"] = mid_angles[swath_offset]
+
+        if ads.name == "SR GR ADS"  and ads.size > 0:
+            srgr_buf = file_buffer[ads.offset:ads.offset + ads.size]
+
+            r = struct.unpack(">ff5f", srgr_buf[13:41])
+
+
+            srgr_coeffs = list(r[2:])
             metadata["srgr_coeffs"] = srgr_coeffs
 
-    # calculate spreading loss compensation
-    c = 299792458
 
+
+
+    
+
+    # antenna gain
+    c = 299792458
     n_samp = gdal_metadata["line_length"]
+    R_first = c * slant_time_first * 1e-9 / 2
     range_spacing = c / (2 * gdal_metadata["records"]["main_processing_params"]["range_samp_rate"])
     range_ref = gdal_metadata["records"]["main_processing_params"]["range_ref"]
-    r_first = c * metadata["slant_time_first"] * 1e-9 / 2
+    R_first = c * slant_time_first * 1e-9 / 2
 
-    spreading_loss = []
-    for n in range(n_samp):
-        r = r_first + n * range_spacing
-        factor = math.sqrt((range_ref / r) ** 3)
-        spreading_loss.append(1 / factor)
+    osv = gdal_metadata["records"]["main_processing_params"]["orbit_state_vectors"]
 
-    cal_factor = gdal_metadata["records"]["main_processing_params"]["calibration_factors"][0]["ext_cal_fact"]
+    sat_x = osv[0]["x_pos_1"] * 1e-2
+    sat_y = osv[0]["y_pos_1"] * 1e-2
+    sat_z = osv[0]["z_pos_1"] * 1e-2
 
-    metadata["cal_factor"] = cal_factor
-    metadata["cal_vector"] = spreading_loss
+    gain_arr = []
+    for n in range(n_samp):
+        # https://github.com/senbox-org/microwave-toolbox/blob/master/sar-op-calibration/src/main/java/eu/esa/sar/calibration/gpf/calibrators/ASARCalibrator.java
+        R = R_first + n * range_spacing
+        n_geogrid = len(lats)
+        geo_interp_idx = (n / n_samp) * (len(lats) - 1)
+        idx = int(geo_interp_idx)
+        fract = geo_interp_idx % 1
+        #interpolate geogrid to match first line geogrid to first OSV
+        lat = lats[idx] + (lats[idx+1] - lats[idx]) * fract
+        lon = lons[idx] + (lons[idx+1] - lons[idx]) * fract
 
-    return metadata
+        sar_dis = math.sqrt(sat_x**2 + sat_y**2 + sat_z**2)
 
+        distance = calc_distance(lat, lon)
 
-def extract_geolocation_metadata(ads: EnvisatADS, file_buffer: bytes, metadata: dict[str, Any]) -> None:
-    """
-    Extract geolocation metadata from the GEOLOCATION GRID ADS record.
+        # elevation angle, cosine law from three sides, slant range, earth center to satellite, earth center to target
+        angle_cos = (R * R + sar_dis * sar_dis - distance * distance) / (2 * R * sar_dis)
+        elev_angle = np.rad2deg(math.acos(angle_cos))
 
-    Args:
-    ----
-        ads: EnvisatADS object containing ADS descriptor information.
-        file_buffer: The full file buffer as bytes.
-        metadata: Dictionary to store extracted metadata.
+        # find the gain from LUT, with 201 gains against reference elevation angle with 0.05 degree steps
+        elev_idx = int((elev_angle - metadata["antenna_ref_elev_angle"]) / 0.05)
+        elev_idx += 100
+        gain = 1.0
+        if elev_idx >= 0 and elev_idx <= len(antenna_gains):
+            # dB -> linear
+            gain = math.pow(10, antenna_gains[elev_idx] / 10)
 
-    Populates the metadata dictionary with:
-        - slant_time_first: First slant range time (float).
-        - incidence_angle_center: Incidence angle at the center (float).
-    """
-    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
-
-    lats_buffer = geoloc_buf[header_size + 11 * 4 * 3: header_size + 11 * 4 * 4]
-    lons_buffer = geoloc_buf[header_size + 11 * 4 * 4: header_size + 11 * 4 * 5]
-
-    lats = list(struct.unpack(">11i", lats_buffer))
-    lons = list(struct.unpack(">11i", lons_buffer))
-    lats = [e * 1e-6 for e in lats]
-    lons = [e * 1e-6 for e in lons]
-
-    # 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
-
-
-def extract_external_calibration_metadata(ads: EnvisatADS, gdal_metadata: dict[str, Any],
-                                          metadata: dict[str, Any]) -> None:
-    """
-    Extract external calibration metadata from the EXTERNAL CALIBRATION ADS record.
+        gain_arr.append(gain)
 
-    Args:
-    ----
-        ads: EnvisatADS object containing ADS descriptor information.
-        gdal_metadata: Dictionary with GDAL metadata required for extraction.
-        metadata: Dictionary to store extracted metadata.
-
-    Populates the metadata dictionary with:
-        - antenna_ref_elev_angle: Reference elevation angle for the antenna (float).
-        - antenna_elev_gains: List of antenna elevation gain values (tuple of 201 floats).
-    """
-    aux_folder = pathlib.Path(os.path.abspath(__file__)).parent
-    aux_folder /= "aux"
-    aux_folder /= "ASAR_Auxiliary_Files"
-    aux_folder /= "ASA_XCA_AX"
 
-    for p in os.scandir(aux_folder):
-        if ads.filename == p.name:
-            with open(p.path, "rb") as fp:
-                ext_cal_buf = fp.read()
-                pol = gdal_metadata["mds1_tx_rx_polar"]
-                swath = gdal_metadata["swath"]
-                swath_offset = ord(swath[2]) - ord("1")
-                pol_offset = {"H/H": 0, "V/V": 1, "H/V": 2, "V/H": 3}[pol]
 
-                # Envisat_Product_Spec_Vol8.pdf
-                # 8.6.2 External Calibration Data
-                offset = 1247 + 378
 
-                offset += 12
-                offset += 4
-                offset += 26 * 28 + 4 * 4
 
-                mid_angles = struct.unpack(">7f", ext_cal_buf[offset:offset + 4 * 7])
 
-                offset += 8 * 4
-
-                offset += 804 * 4 * swath_offset
-                offset += 201 * 4 * pol_offset
+    # calculate spreading loss compensation
 
-                antenna_gains = struct.unpack(">201f", ext_cal_buf[offset:offset + 4 * 201])
-                metadata["antenna_ref_elev_angle"] = mid_angles[swath_offset]
-                metadata["antenna_elev_gains"] = antenna_gains
+    spreading_loss = []
+    for n in range(n_samp):
+        R = R_first + n * range_spacing
+        factor = math.sqrt((range_ref / R) ** 3)
+        spreading_loss.append(1/factor)
 
+    cal_factor = gdal_metadata["records"]["main_processing_params"]["calibration_factors"][0]["ext_cal_fact"]
 
-def extract_srgr_coeffs(ads: EnvisatADS, file_buffer: bytes) -> list[float]:
-    """
-    Extract SRGR (Slant Range to Ground Range) coefficients from the SR GR ADS record.
+    metadata["cal_factor"] = cal_factor
+    metadata["cal_vector"] = np.array(spreading_loss) * np.array(gain_arr)
 
-    Args:
-    ----
-        ads: EnvisatADS object containing ADS descriptor information.
-        file_buffer: The full file buffer as bytes.
 
-    Returns: List of SRGR coefficients (list of 5 floats).
-    """
-    srgr_buf = file_buffer[ads.offset:ads.offset + ads.size]
-    r = struct.unpack(">ff5f", srgr_buf[13:41])
-    srgr_coeffs = list(r[2:])
-    return srgr_coeffs
+    return metadata