Fix issue with improper geolocation of IMP products.

Fix antenna gain with IMS mode.

Author: priit111

asar_xarray/asar.py

@@ -89,10 +89,16 @@ def create_dataset(metadata: dict[str, Any], filepath: str) -> xr.Dataset:
     range_sampling_rate = metadata["records"]["main_processing_params"]["range_samp_rate"]
     image_slant_range_time = metadata["direct_parse"]["slant_time_first"] * 1e-9
 
-    if metadata["sph_descriptor"] == "Image Mode SLC Image":
+    range_pixel_spacing = metadata["records"]["main_processing_params"]["range_spacing"]
+
+    product_str = metadata["sph_descriptor"]
+
+    if product_str == "Image Mode SLC Image":
         product_type = "SLC"
+    elif product_str == "Image Mode Precision Image":
+        product_type = "GRD"
     else:
-        raise RuntimeError("Only Image mode SLC files(IMS) supported for now")
+        raise RuntimeError("Only Image mode files supported(IMS & IMP) supported for now")
 
     attrs = {
         "family_name": "Envisat",
@@ -106,8 +112,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_pixel_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"],
@@ -116,15 +122,14 @@ 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] = {
@@ -140,12 +145,34 @@ def create_dataset(metadata: dict[str, Any], filepath: str) -> xr.Dataset:
         ),
     }
 
-    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)
+    if product_type == "SLC":
+        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)
+    else:
+
+        # generate slant range times from ground ranges for ground range products
+        # this is easier due to the Envisat only having GRSR metadata unlike Sentinel1
+        # however xarray/Sarsen handle GRD conversion differently for S1
+        ground_range = np.linspace(
+            0,
+            range_pixel_spacing * (number_of_samples - 1),
+            number_of_samples,
+        )
+
+        # numpy polyval expects the polynomial top be highest ranked first
+        coeffs = list(reversed(metadata["direct_parse"]["grsr_coeffs"]))
+
+        slant_ranges = np.polyval(coeffs, ground_range)
+        slant_ranges *= 2
+
+        c = 299792458
+        slant_range_times = slant_ranges / c
+
+        coords["slant_range_time"] = ("pixel", slant_range_times)
 
     data = xr.open_dataarray(filepath, engine='rasterio')
     data.encoding.clear()
@@ -180,7 +207,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['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/envisat_direct.py

@@ -175,9 +175,10 @@ def parse_direct(path: str, gdal_metadata) -> dict[str, Any]:
 
             r = struct.unpack(">ff5f", srgr_buf[13:41])
 
-
+            # Envisat file specification says it is SR/GR Conversion ADSR,
+            # however the polynomial is for GR to SR...
             srgr_coeffs = list(r[2:])
-            metadata["srgr_coeffs"] = srgr_coeffs
+            metadata["grsr_coeffs"] = srgr_coeffs
 
 
 
@@ -187,59 +188,69 @@ def parse_direct(path: str, gdal_metadata) -> dict[str, Any]:
     # 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 * slant_time_first * 1e-9 / 2
-
-    osv = gdal_metadata["records"]["main_processing_params"]["orbit_state_vectors"]
 
-    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
 
     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
+    if gdal_metadata["sample_type"] == "DETECTED":
+        gain_arr = np.ones(n_samp)
+        spreading_loss = np.ones(n_samp)
+    else:
+        # antenna gain
+        c = 299792458
+        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 * slant_time_first * 1e-9 / 2
+
+        osv = gdal_metadata["records"]["main_processing_params"]["orbit_state_vectors"]
+
+        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
+
+        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
+
+            sar_dis = math.sqrt(sat_x**2 + sat_y**2 + sat_z**2)
+
+            distance = calc_distance(lat, lon)
 
-        sar_dis = math.sqrt(sat_x**2 + sat_y**2 + sat_z**2)
+            # 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))
 
-        distance = calc_distance(lat, lon)
+            # 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)
 
-        # 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))
+            gain_arr.append(1/gain)
 
-        # 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)
 
-        gain_arr.append(gain)
 
 
 
 
+        # calculate spreading loss compensation
 
+        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)
 
-    # calculate spreading loss compensation
 
-    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"]