From cada77fd86fe44ab89b59ec738fbb00229844003 Mon Sep 17 00:00:00 2001
From: Liang Yu <liangjyu@gmail.com>
Date: Tue, 10 May 2022 17:03:23 -0700
Subject: [PATCH 1/4] initial running but unverified rdr2geo test

---
 tests/prep.py         | 180 ++++++++++++++++++++++++++++++++++++++++++
 tests/test_rdr2geo.py |  13 +++
 2 files changed, 193 insertions(+)
 create mode 100644 tests/prep.py
 create mode 100644 tests/test_rdr2geo.py

diff --git a/tests/prep.py b/tests/prep.py
new file mode 100644
index 00000000..614371eb
--- /dev/null
+++ b/tests/prep.py
@@ -0,0 +1,180 @@
+import datetime
+import pytest
+from types import SimpleNamespace
+
+import isce3
+import numpy as np
+from osgeo import gdal
+from s1reader.s1_burst_slc import Doppler, Sentinel1BurstSlc
+
+def zero_dem(params):
+    lon0 = params.lon0 - params.omega / 2
+    lat0 = params.lat0 - params.omega / 2
+    geo_trans = [lon0, params.omega, 0.0, lat0, 0.0, params.omega]
+
+    length = 1
+    width = params.n_sv + 1
+
+    # create DEM file
+    dem_raster = isce3.io.Raster(params.dem, width, length, 1,
+                                 gdal.GDT_Float32, 'GTiff')
+    dem_raster.set_geotransform(geo_trans)
+    dem_raster.set_epsg(params.epsg)
+    del dem_raster
+
+    # write zero to DEM file
+    data = np.zeros((length, width))
+    ds = gdal.Open(params.dem, gdal.GA_Update)
+    ds.GetRasterBand(1).WriteArray(data)
+    ds = None
+
+
+def solve_geocentric_lat(params):
+    temp = 1 + params.h_sat / params.ellipsoid.a
+    temp1 = rng / params.ellipsoid.a
+    temp2 = rng / (params.ellipsoid.a + params.h_sat)
+
+    cosang = 0.5 * (temp + (1.0/temp) - temp1 * temp2)
+    angdiff = np.arccos(cosang);
+
+    if (params.look_side * satomega > 0):
+        x = satlat0 + angdiff
+    else:
+        x = satlat0 - angdiff
+
+    return x
+
+
+def compute_orbit(params):
+    # Setup orbit
+    statevecs = []
+    clat = np.cos(params.lat0)
+    slat = np.sin(params.lat0)
+    sat_h = params.ellipsoid.a + params.h_sat
+    for ii in range(params.n_sv):
+        total_delta_t = ii * params.orbit_dt
+        t = params.sensing_start + datetime.timedelta(seconds=total_delta_t)
+
+        lon = params.lon0 + params.omega * total_delta_t
+
+        pos = [sat_h * clat * np.cos(lon),
+                sat_h * clat * np.sin(lon),
+                sat_h * slat]
+
+        vel = [-params.omega * pos[1],
+                params.omega * pos[0],
+                0.0]
+
+        sv = isce3.core.StateVector(isce3.core.DateTime(t), pos, vel)
+        statevecs.append(sv)
+
+    # use list of stateVectors to init and return isce3.core.Orbit
+    time_delta = datetime.timedelta(days=2)
+    ref_epoch = isce3.core.DateTime(params.sensing_start - time_delta)
+
+    return isce3.core.Orbit(statevecs, ref_epoch)
+
+
+def common_params(lat):
+    params = SimpleNamespace()
+    sensing_start = "2017-02-12T01:12:30.0"
+    fmt = "%Y-%m-%dT%H:%M:%S.%f"
+    params.sensing_start = datetime.datetime.strptime(sensing_start, fmt)
+    params.lon0 = 0.0
+    params.lat0 = np.radians(lat)
+    params.omega = np.radians(0.1)
+    params.n_sv = 10
+    params.n_samples = 20
+    params.h_sat = 700000.0
+    params.orbit_dt = 10.0
+    params.epsg = 4326
+    params.look_side = isce3.core.LookSide.Right
+    params.ellipsoid = isce3.core.make_projection(params.epsg).ellipsoid
+    params.orbit = compute_orbit(params)
+    params.dem = "zero_dem.tif"
+    return params
+
+
+'''
+def create_burst(ref_epoch_str="2017-02-12T01:12:30.0", az_time_interval):
+    fmt = "%Y-%m-%dT%H:%M:%S.%f"
+'''
+def create_burst(params):
+
+    # place holder value
+    dont_matter_for_now = 0
+
+    radar_freq = dont_matter_for_now
+    wavelength = 0.24
+    azimuth_steer_rate = 0.024389943375862838
+    azimuth_time_interval = 2.0
+    slant_range_time = dont_matter_for_now
+    starting_range = 800000.0
+    iw2_mid_range = dont_matter_for_now
+    range_sampling_rate = dont_matter_for_now
+    range_pxl_spacing = 10.0
+    shape = (1, params.n_samples)
+    az_fm_rate = dont_matter_for_now
+    doppler = Doppler(isce3.core.Poly1d([1]), isce3.core.LUT2d())
+    rng_processing_bandwidth = dont_matter_for_now
+    pol = 'vv'
+    burst_id = 't00_iw1_b000'
+    platform_id = dont_matter_for_now
+    center_pts = ()
+    boundary_pts = []
+    orbit_dir = str(dont_matter_for_now)
+    tiff_path = dont_matter_for_now
+    i_burst = dont_matter_for_now
+    first_valid_sample = dont_matter_for_now
+    last_sample = dont_matter_for_now
+    first_valid_line = dont_matter_for_now
+    last_line = dont_matter_for_now
+    range_window_type = dont_matter_for_now
+    range_window_coeff = dont_matter_for_now
+    i_burst = dont_matter_for_now
+    range_window_type = str(dont_matter_for_now)
+    range_window_coeff = dont_matter_for_now
+    rank = int(dont_matter_for_now)
+    prf_raw_data = dont_matter_for_now
+
+    burst = Sentinel1BurstSlc(params.sensing_start, radar_freq, wavelength,
+                              azimuth_steer_rate, azimuth_time_interval,
+                              slant_range_time, starting_range, iw2_mid_range,
+                              range_sampling_rate, range_pxl_spacing,
+                              shape, az_fm_rate, doppler,
+                              rng_processing_bandwidth, pol, burst_id,
+                              platform_id, center_pts,
+                              boundary_pts, params.orbit, orbit_dir,
+                              tiff_path, i_burst, first_valid_sample,
+                              last_sample, first_valid_line, last_line,
+                              range_window_type, range_window_coeff,
+                              rank, prf_raw_data)
+
+    return burst
+
+
+def cfg_45_lat():
+    cfg = SimpleNamespace()
+
+    params = common_params(45)
+
+    # create and set DEM raster
+    zero_dem(params)
+
+    # set rdr2geo params
+    rdr2geo_params = SimpleNamespace()
+    rdr2geo_params.threshold = 1e-8
+    rdr2geo_params.numiter= 50
+    rdr2geo_params.lines_per_block = 100
+    rdr2geo_params.extraiter = 25
+    rdr2geo_params.compute_mask = True
+    cfg.rdr2geo_params = rdr2geo_params
+
+    # make list with single burst
+    cfg.bursts = [create_burst(params)]
+    cfg.dem = params.dem
+    cfg.gpu_enabled = False
+    cfg.gpu_id = 0
+    cfg.product_path = 'test_out'
+
+    return cfg
diff --git a/tests/test_rdr2geo.py b/tests/test_rdr2geo.py
new file mode 100644
index 00000000..1e646e73
--- /dev/null
+++ b/tests/test_rdr2geo.py
@@ -0,0 +1,13 @@
+#!/usr/bin/env python3
+import pytest
+
+from compass import s1_rdr2geo
+import prep
+
+def test_45_lat_rdr2geo():
+    cfg = prep.cfg_45_lat()
+    s1_rdr2geo.run(cfg)
+
+
+if __name__ == '__main__':
+    test_45_lat_rdr2geo()

From 5c8f857141f5b808362bef9627de09a53c724792 Mon Sep 17 00:00:00 2001
From: Liang Yu <liangjyu@gmail.com>
Date: Tue, 17 May 2022 13:00:40 -0700
Subject: [PATCH 2/4] add unit test comparison (broken) merge with main

---
 src/compass/s1_cslc.py             |  47 +++++++----
 src/compass/s1_geo2rdr.py          |   3 +-
 src/compass/s1_geocode_slc.py      | 127 +++++++++++++++++++++++++++++
 src/compass/s1_rdr2geo.py          |   3 +-
 src/compass/s1_resample_burst.py   |   3 +-
 src/compass/utils/yaml_argparse.py |  11 ++-
 tests/prep.py                      |  53 ++++++++----
 tests/test_rdr2geo.py              |  17 ++++
 8 files changed, 225 insertions(+), 39 deletions(-)
 create mode 100644 src/compass/s1_geocode_slc.py

diff --git a/src/compass/s1_cslc.py b/src/compass/s1_cslc.py
index 812632da..38698cd1 100644
--- a/src/compass/s1_cslc.py
+++ b/src/compass/s1_cslc.py
@@ -1,25 +1,40 @@
-from compass import s1_rdr2geo, s1_geo2rdr, s1_resample_burst
+from compass import s1_rdr2geo, s1_geo2rdr, s1_resample, s1_geocode_slc
+from compass.utils.geo_runconfig import GeoRunConfig
 from compass.utils.runconfig import RunConfig
 from compass.utils.yaml_argparse import YamlArgparse
 
 
-def run(cfg):
-    # If boolean flag "is_reference" is true
-    # Run rdr2geo and archive reference burst
-    if cfg.is_reference:
-        s1_rdr2geo.run(cfg)
-    else:
-        s1_geo2rdr.run(cfg)
-        s1_resample_burst.run(cfg)
+def main(run_config_path, grid_type):
+    if grid_type == 'radar':
+        # CSLC workflow in radar coordinates
+        # get a runconfig dict from command line args
+        cfg = RunConfig.load_from_yaml(parser.run_config_path, 's1_cslc_radar')
+
+        if cfg.is_reference:
+            # reference burst - run rdr2geo and archive it
+            s1_rdr2geo.run(cfg)
+
+        else:
+            # secondary burst - run geo2rdr + resample
+            s1_geo2rdr.run(cfg)
+            s1_resample.run(cfg)
+
+    elif proc_steps_geo.issubset(proc_steps):
+        # CSLC workflow in geo-coordinates
+        # get a runconfig dict from command line argumens
+        cfg = GeoRunConfig.load_from_yaml(parser.run_config_path, 's1_cslc_geo')
+
+        # run geocode_slc
+        s1_geocode_slc.run(cfg)
 
 
 if __name__ == "__main__":
-    '''run rdr2geo from command line'''
+    '''run s1_cslc from command line'''
     # load command line args
-    arg_parser = YamlArgparse()
-
-    # get a runconfig dict from command line args
-    runconfig = RunConfig.load_from_yaml(arg_parser.run_config_path, 'cslc_s1')
+    parser = YamlArgparse()
+    parser.parser.add_argument('--grid-type', dest='grid_type', type=str,
+                               choices=['geo', 'radar'], default='geo',
+                               help='Grid type to perform CSLC processing in\n')
+    parser.parse()
 
-    # run rdr2geo
-    run(runconfig)
+    main(parser.run_config_path, parser.args.grid_type)
diff --git a/src/compass/s1_geo2rdr.py b/src/compass/s1_geo2rdr.py
index 1b74ff3d..ac05e323 100644
--- a/src/compass/s1_geo2rdr.py
+++ b/src/compass/s1_geo2rdr.py
@@ -93,7 +93,8 @@ def run(cfg: dict):
 
 if __name__ == "__main__":
     """Run geo2rdr from command line"""
-    geo2rdr_parser = YamlArgparse()
+    parser = YamlArgparse()
+    parser.parse()
 
     # Get a runconfig dict from command line arguments
     geo2rdr_runconfig = RunConfig.load_from_yaml(
diff --git a/src/compass/s1_geocode_slc.py b/src/compass/s1_geocode_slc.py
new file mode 100644
index 00000000..deae228a
--- /dev/null
+++ b/src/compass/s1_geocode_slc.py
@@ -0,0 +1,127 @@
+from datetime import timedelta
+import os
+import time
+
+import isce3
+import journal
+import numpy as np
+from osgeo import gdal
+
+from compass.utils.geo_metadata import GeoCslcMetadata
+from compass.utils.geo_runconfig import GeoRunConfig
+from compass.utils.helpers import get_module_name
+from compass.utils.range_split_spectrum import range_split_spectrum
+from compass.utils.yaml_argparse import YamlArgparse
+
+
+def run(cfg):
+    '''
+    Run geocode burst workflow with user-defined
+    args stored in dictionary runconfig *cfg*
+
+    Parameters
+    ---------
+    cfg: dict
+        Dictionary with user runconfig options
+    '''
+    module_name = get_module_name(__file__)
+    info_channel = journal.info(f"{module_name}.run")
+    info_channel.log(f"Starting {module_name} burst")
+
+    # Start tracking processing time
+    t_start = time.time()
+
+    # Common initializations
+    dem_raster = isce3.io.Raster(cfg.dem)
+    epsg = dem_raster.get_epsg()
+    proj = isce3.core.make_projection(epsg)
+    ellipsoid = proj.ellipsoid
+    image_grid_doppler = isce3.core.LUT2d()
+    threshold = cfg.geo2rdr_params.threshold
+    iters = cfg.geo2rdr_params.numiter
+    blocksize = cfg.geo2rdr_params.lines_per_block
+    dem_margin = cfg.geocoding_params.dem_margin
+    flatten = cfg.geocoding_params.flatten
+
+    # process one burst only
+    burst = cfg.bursts[0]
+    date_str = burst.sensing_start.strftime("%Y%m%d")
+    burst_id = burst.burst_id
+    pol = burst.polarization
+    geo_grid = cfg.geogrids[burst_id]
+
+    os.makedirs(cfg.output_dir, exist_ok=True)
+
+    scratch_path = f'{cfg.scratch_path}/{burst_id}/{date_str}'
+    os.makedirs(scratch_path, exist_ok=True)
+
+    radar_grid = burst.as_isce3_radargrid()
+    native_doppler = burst.doppler.lut2d
+    orbit = burst.orbit
+
+    # Get azimuth polynomial coefficients for this burst
+    az_carrier_poly2d = burst.get_az_carrier_poly()
+
+    # Split the range bandwidth of the burst, if required
+    if cfg.split_spectrum_params.enabled:
+        rdr_burst_raster = range_split_spectrum(burst,
+                                                cfg.split_spectrum_params,
+                                                scratch_path)
+    else:
+        temp_slc_path = f'{scratch_path}/{burst_id}_{pol}_temp.vrt'
+        burst.slc_to_vrt_file(temp_slc_path)
+        rdr_burst_raster = isce3.io.Raster(temp_slc_path)
+
+    # Generate output geocoded burst raster
+    geo_burst_raster = isce3.io.Raster(
+        f'{cfg.output_dir}/{burst_id}_{date_str}_{pol}.slc',
+        geo_grid.width, geo_grid.length,
+        rdr_burst_raster.num_bands, gdal.GDT_CFloat32,
+        cfg.geocoding_params.output_format)
+
+    # Extract burst boundaries
+    b_bounds = np.s_[burst.first_valid_line:burst.last_valid_line,
+               burst.first_valid_sample:burst.last_valid_sample]
+
+    # Create sliced radar grid representing valid region of the burst
+    sliced_radar_grid = burst.as_isce3_radargrid()[b_bounds]
+
+    # Geocode
+    isce3.geocode.geocode_slc(geo_burst_raster, rdr_burst_raster,
+                              dem_raster,
+                              radar_grid, sliced_radar_grid,
+                              geo_grid, orbit,
+                              native_doppler,
+                              image_grid_doppler, ellipsoid, threshold,
+                              iters,
+                              blocksize, dem_margin, flatten,
+                              azimuth_carrier=az_carrier_poly2d)
+
+    # Set geo transformation
+    geotransform = [geo_grid.start_x, geo_grid.spacing_x, 0,
+                    geo_grid.start_y, 0, geo_grid.spacing_y]
+    geo_burst_raster.set_geotransform(geotransform)
+    geo_burst_raster.set_epsg(epsg)
+    del geo_burst_raster
+
+    # Save burst metadata
+    metadata = GeoCslcMetadata.from_georunconfig(cfg)
+    json_path = f'{cfg.output_dir}/{burst_id}_{date_str}_{pol}.json'
+    with open(json_path, 'w') as f_json:
+        metadata.to_file(f_json, 'json')
+
+    dt = str(timedelta(seconds=time.time() - t_start)).split(".")[0]
+    info_channel.log(f"{module_name} burst successfully ran in {dt} (hr:min:sec)")
+
+
+if __name__ == "__main__":
+    '''Run geocode cslc workflow from command line'''
+    # load arguments from command line
+    parser = YamlArgparse()
+    parser.parse()
+
+    # Get a runconfig dict from command line argumens
+    cfg = GeoRunConfig.load_from_yaml(parser.run_config_path, 's1_cslc_geo')
+
+    # Run geocode burst workflow
+    run(cfg)
diff --git a/src/compass/s1_rdr2geo.py b/src/compass/s1_rdr2geo.py
index 9d706817..08a985a5 100644
--- a/src/compass/s1_rdr2geo.py
+++ b/src/compass/s1_rdr2geo.py
@@ -115,7 +115,8 @@ def run(cfg):
 if __name__ == "__main__":
     '''run rdr2geo from command line'''
     # load command line args
-    rdr2geo_parser = YamlArgparse()
+    parser = YamlArgparse()
+    parser.parse()
 
     # get a runconfig dict from command line args
     rdr2geo_runconfig = RunConfig.load_from_yaml(rdr2geo_parser.args.run_config_path, 'rdr2geo')
diff --git a/src/compass/s1_resample_burst.py b/src/compass/s1_resample_burst.py
index c7532654..8d55cfb3 100644
--- a/src/compass/s1_resample_burst.py
+++ b/src/compass/s1_resample_burst.py
@@ -99,7 +99,8 @@ def run(cfg: dict):
 
 if __name__ == "__main__":
     """Run resample burst from command line"""
-    resample_parser = YamlArgparse()
+    parser = YamlArgparse()
+    parser.parse()
 
     # Get a runconfig dict from command line arguments
     resample_runconfig = RunConfig.load_from_yaml(
diff --git a/src/compass/utils/yaml_argparse.py b/src/compass/utils/yaml_argparse.py
index 28a992db..6cc0bcd2 100644
--- a/src/compass/utils/yaml_argparse.py
+++ b/src/compass/utils/yaml_argparse.py
@@ -2,10 +2,13 @@
 
 class YamlArgparse():
     def __init__(self):
-        '''Initialize YamlArgparse class and parse CLI arguments for COMPASS'''
-        parser = argparse.ArgumentParser(description='', formatter_class=argparse.ArgumentDefaultsHelpFormatter)
-        parser.add_argument('run_config_path', type=str, nargs='?', default=None, help='Path to run config file')
-        self.args = parser.parse_args()
+        '''Initialize YamlArgparse class with basics'''
+        self.parser = argparse.ArgumentParser(description='', formatter_class=argparse.ArgumentDefaultsHelpFormatter)
+        self.parser.add_argument('run_config_path', type=str, nargs='?', default=None, help='Path to run config file')
+
+    def parse():
+        '''Parse CLI arguments for COMPASS'''
+        self.args = self.parser.parse_args()
 
     @property
     def run_config_path(self) -> str:
diff --git a/tests/prep.py b/tests/prep.py
index 614371eb..10051539 100644
--- a/tests/prep.py
+++ b/tests/prep.py
@@ -29,18 +29,20 @@ def zero_dem(params):
     ds = None
 
 
-def solve_geocentric_lat(params):
+def solve_geocentric_lat(slant_range, params):
     temp = 1 + params.h_sat / params.ellipsoid.a
-    temp1 = rng / params.ellipsoid.a
-    temp2 = rng / (params.ellipsoid.a + params.h_sat)
+    temp1 = slant_range / params.ellipsoid.a
+    temp2 = slant_range / (params.ellipsoid.a + params.h_sat)
 
     cosang = 0.5 * (temp + (1.0/temp) - temp1 * temp2)
     angdiff = np.arccos(cosang);
 
-    if (params.look_side * satomega > 0):
-        x = satlat0 + angdiff
+    look_side_int = 1 if params.look_side == isce3.core.LookSide.Left \
+        else -1
+    if (look_side_int * params.omega > 0):
+        x = params.lat0 + angdiff
     else:
-        x = satlat0 - angdiff
+        x = params.lat0 - angdiff
 
     return x
 
@@ -86,6 +88,10 @@ def common_params(lat):
     params.n_sv = 10
     params.n_samples = 20
     params.h_sat = 700000.0
+    params.slant_range0 = 80000.0
+    params.d_slant_range = 10.0
+    params.az_time0 = 5.0
+    params.d_az_time = 2.0
     params.orbit_dt = 10.0
     params.epsg = 4326
     params.look_side = isce3.core.LookSide.Right
@@ -107,12 +113,9 @@ def create_burst(params):
     radar_freq = dont_matter_for_now
     wavelength = 0.24
     azimuth_steer_rate = 0.024389943375862838
-    azimuth_time_interval = 2.0
     slant_range_time = dont_matter_for_now
-    starting_range = 800000.0
     iw2_mid_range = dont_matter_for_now
     range_sampling_rate = dont_matter_for_now
-    range_pxl_spacing = 10.0
     shape = (1, params.n_samples)
     az_fm_rate = dont_matter_for_now
     doppler = Doppler(isce3.core.Poly1d([1]), isce3.core.LUT2d())
@@ -138,9 +141,9 @@ def create_burst(params):
     prf_raw_data = dont_matter_for_now
 
     burst = Sentinel1BurstSlc(params.sensing_start, radar_freq, wavelength,
-                              azimuth_steer_rate, azimuth_time_interval,
-                              slant_range_time, starting_range, iw2_mid_range,
-                              range_sampling_rate, range_pxl_spacing,
+                              azimuth_steer_rate, params.d_az_time,
+                              slant_range_time, params.slant_range0, iw2_mid_range,
+                              range_sampling_rate, params.d_slant_range,
                               shape, az_fm_rate, doppler,
                               rng_processing_bandwidth, pol, burst_id,
                               platform_id, center_pts,
@@ -153,13 +156,31 @@ def create_burst(params):
     return burst
 
 
+def compute_expected_llh(test_params):
+    i = np.arange(test_params.n_samples)
+    ellipsoid = test_params.ellipsoid
+    az_time = test_params.az_time0 + i * test_params.d_az_time
+    slant_range = test_params.slant_range0 + i * test_params.d_slant_range
+    lon = test_params.lon0 + test_params.omega * az_time
+    c_lon = np.cos(lon)
+    s_lon = np.sin(lon)
+    lat = np.arccos(solve_geocentric_lat(slant_range, test_params))
+    c_lat = np.cos(lat)
+    s_lat = np.sin(lat)
+    xyz_vec = np.vstack([ellipsoid.a * c_lat * c_lon,
+                         ellipsoid.a * c_lat * s_lon,
+                         ellipsoid.b * s_lat])
+    llh = [ellipsoid.xyz_to_lon_lat(xyz_pt) for xyz_pt in xyz_vec]
+    return llh
+
+
 def cfg_45_lat():
     cfg = SimpleNamespace()
 
-    params = common_params(45)
+    cfg.test_params = common_params(45)
 
     # create and set DEM raster
-    zero_dem(params)
+    zero_dem(cfg.test_params)
 
     # set rdr2geo params
     rdr2geo_params = SimpleNamespace()
@@ -171,8 +192,8 @@ def cfg_45_lat():
     cfg.rdr2geo_params = rdr2geo_params
 
     # make list with single burst
-    cfg.bursts = [create_burst(params)]
-    cfg.dem = params.dem
+    cfg.bursts = [create_burst(cfg.test_params)]
+    cfg.dem = cfg.test_params.dem
     cfg.gpu_enabled = False
     cfg.gpu_id = 0
     cfg.product_path = 'test_out'
diff --git a/tests/test_rdr2geo.py b/tests/test_rdr2geo.py
index 1e646e73..860843d7 100644
--- a/tests/test_rdr2geo.py
+++ b/tests/test_rdr2geo.py
@@ -2,11 +2,28 @@
 import pytest
 
 from compass import s1_rdr2geo
+import numpy as np
+from osgeo import gdal
+
 import prep
 
+def get_rdr2geo_output(cfg):
+    burst = cfg.bursts[0]
+    burst_id = burst.burst_id
+    date_str = str(burst.sensing_start.date())
+    topo_path = f'{cfg.product_path}/{burst_id}/{date_str}/topo.vrt'
+    print(topo_path)
+    ds = gdal.Open(topo_path)
+    llh = np.vstack([ds.GetRasterBand(1).ReadAsArray() for i in range(3)])
+    return llh
+
 def test_45_lat_rdr2geo():
     cfg = prep.cfg_45_lat()
+
+    expected_llh = prep.compute_expected_llh(cfg.test_params)
+
     s1_rdr2geo.run(cfg)
+    computed_llh = get_rdr2geo_output(cfg)
 
 
 if __name__ == '__main__':

From 312faea69e8c1f4715b5c692c5fdc9cbed9c119c Mon Sep 17 00:00:00 2001
From: Liang Yu <liangjyu@gmail.com>
Date: Wed, 8 Jun 2022 21:04:04 -0700
Subject: [PATCH 3/4] fix slant range typo differentiate start times removed
 double arccos fixed raster band index

---
 tests/prep.py         | 23 ++++++++++++++---------
 tests/test_rdr2geo.py |  4 +++-
 2 files changed, 17 insertions(+), 10 deletions(-)

diff --git a/tests/prep.py b/tests/prep.py
index 10051539..caf28899 100644
--- a/tests/prep.py
+++ b/tests/prep.py
@@ -55,13 +55,13 @@ def compute_orbit(params):
     sat_h = params.ellipsoid.a + params.h_sat
     for ii in range(params.n_sv):
         total_delta_t = ii * params.orbit_dt
-        t = params.sensing_start + datetime.timedelta(seconds=total_delta_t)
+        t = params.orbit_start + datetime.timedelta(seconds=total_delta_t)
 
         lon = params.lon0 + params.omega * total_delta_t
 
         pos = [sat_h * clat * np.cos(lon),
-                sat_h * clat * np.sin(lon),
-                sat_h * slat]
+               sat_h * clat * np.sin(lon),
+               sat_h * slat]
 
         vel = [-params.omega * pos[1],
                 params.omega * pos[0],
@@ -79,16 +79,18 @@ def compute_orbit(params):
 
 def common_params(lat):
     params = SimpleNamespace()
-    sensing_start = "2017-02-12T01:12:30.0"
+    sensing_start = "2017-02-12T01:12:35.0"
+    orbit_start = "2017-02-12T01:12:30.0"
     fmt = "%Y-%m-%dT%H:%M:%S.%f"
     params.sensing_start = datetime.datetime.strptime(sensing_start, fmt)
+    params.orbit_start = datetime.datetime.strptime(orbit_start, fmt)
     params.lon0 = 0.0
     params.lat0 = np.radians(lat)
     params.omega = np.radians(0.1)
     params.n_sv = 10
     params.n_samples = 20
     params.h_sat = 700000.0
-    params.slant_range0 = 80000.0
+    params.slant_range0 = 800000.0
     params.d_slant_range = 10.0
     params.az_time0 = 5.0
     params.d_az_time = 2.0
@@ -139,6 +141,7 @@ def create_burst(params):
     range_window_coeff = dont_matter_for_now
     rank = int(dont_matter_for_now)
     prf_raw_data = dont_matter_for_now
+    range_chirp_ramp_rate = dont_matter_for_now
 
     burst = Sentinel1BurstSlc(params.sensing_start, radar_freq, wavelength,
                               azimuth_steer_rate, params.d_az_time,
@@ -151,7 +154,7 @@ def create_burst(params):
                               tiff_path, i_burst, first_valid_sample,
                               last_sample, first_valid_line, last_line,
                               range_window_type, range_window_coeff,
-                              rank, prf_raw_data)
+                              rank, prf_raw_data, range_chirp_ramp_rate)
 
     return burst
 
@@ -164,13 +167,14 @@ def compute_expected_llh(test_params):
     lon = test_params.lon0 + test_params.omega * az_time
     c_lon = np.cos(lon)
     s_lon = np.sin(lon)
-    lat = np.arccos(solve_geocentric_lat(slant_range, test_params))
+    lat = solve_geocentric_lat(slant_range, test_params)
     c_lat = np.cos(lat)
     s_lat = np.sin(lat)
     xyz_vec = np.vstack([ellipsoid.a * c_lat * c_lon,
                          ellipsoid.a * c_lat * s_lon,
-                         ellipsoid.b * s_lat])
-    llh = [ellipsoid.xyz_to_lon_lat(xyz_pt) for xyz_pt in xyz_vec]
+                         ellipsoid.a * s_lat])
+    llh = [ellipsoid.xyz_to_lon_lat(xyz_vec[:,i])
+           for i in range(test_params.n_samples)]
     return llh
 
 
@@ -193,6 +197,7 @@ def cfg_45_lat():
 
     # make list with single burst
     cfg.bursts = [create_burst(cfg.test_params)]
+    import ipdb; ipdb.set_trace()
     cfg.dem = cfg.test_params.dem
     cfg.gpu_enabled = False
     cfg.gpu_id = 0
diff --git a/tests/test_rdr2geo.py b/tests/test_rdr2geo.py
index 860843d7..bb5f8205 100644
--- a/tests/test_rdr2geo.py
+++ b/tests/test_rdr2geo.py
@@ -14,7 +14,7 @@ def get_rdr2geo_output(cfg):
     topo_path = f'{cfg.product_path}/{burst_id}/{date_str}/topo.vrt'
     print(topo_path)
     ds = gdal.Open(topo_path)
-    llh = np.vstack([ds.GetRasterBand(1).ReadAsArray() for i in range(3)])
+    llh = np.vstack([ds.GetRasterBand(i+1).ReadAsArray() for i in range(3)])
     return llh
 
 def test_45_lat_rdr2geo():
@@ -24,6 +24,8 @@ def test_45_lat_rdr2geo():
 
     s1_rdr2geo.run(cfg)
     computed_llh = get_rdr2geo_output(cfg)
+    print(computed_llh)
+    print(expected_llh)
 
 
 if __name__ == '__main__':

From 45c807c57d4705957821dafffff2608c2eb021bf Mon Sep 17 00:00:00 2001
From: Liang Yu <liangjyu@gmail.com>
Date: Thu, 7 Jul 2022 21:53:51 -0700
Subject: [PATCH 4/4] attempt to use expected height DEM omega 0.1 -> 0.001
 save expected height in cfg add more comments

---
 tests/prep.py         | 40 ++++++++++++++++++++++------------------
 tests/test_rdr2geo.py | 15 +++++++++------
 2 files changed, 31 insertions(+), 24 deletions(-)

diff --git a/tests/prep.py b/tests/prep.py
index caf28899..a9d8eaa1 100644
--- a/tests/prep.py
+++ b/tests/prep.py
@@ -7,23 +7,27 @@
 from osgeo import gdal
 from s1reader.s1_burst_slc import Doppler, Sentinel1BurstSlc
 
-def zero_dem(params):
-    lon0 = params.lon0 - params.omega / 2
-    lat0 = params.lat0 - params.omega / 2
-    geo_trans = [lon0, params.omega, 0.0, lat0, 0.0, params.omega]
+def expected_height_dem(cfg):
+    params = cfg.test_params
+
+    # define geogrid in degrees
+    lon0 = np.degrees(params.lon0) #- params.omega / 2
+    lat0 = np.degrees(params.lat0) #- params.omega / 2
+    omega = np.degrees(params.omega)
+    geo_trans = [lon0, omega, 0.0, lat0, 0.0, omega]
 
     length = 1
-    width = params.n_sv + 1
+    width = params.n_samples#+ 1
 
-    # create DEM file
+    # create DEM empty file
     dem_raster = isce3.io.Raster(params.dem, width, length, 1,
                                  gdal.GDT_Float32, 'GTiff')
     dem_raster.set_geotransform(geo_trans)
     dem_raster.set_epsg(params.epsg)
     del dem_raster
 
-    # write zero to DEM file
-    data = np.zeros((length, width))
+    # populate empty DEM with expected LLH values
+    data = np.array(cfg.expected_llh[:,2])[np.newaxis, ...]
     ds = gdal.Open(params.dem, gdal.GA_Update)
     ds.GetRasterBand(1).WriteArray(data)
     ds = None
@@ -86,7 +90,7 @@ def common_params(lat):
     params.orbit_start = datetime.datetime.strptime(orbit_start, fmt)
     params.lon0 = 0.0
     params.lat0 = np.radians(lat)
-    params.omega = np.radians(0.1)
+    params.omega = np.radians(0.001)
     params.n_sv = 10
     params.n_samples = 20
     params.h_sat = 700000.0
@@ -99,14 +103,10 @@ def common_params(lat):
     params.look_side = isce3.core.LookSide.Right
     params.ellipsoid = isce3.core.make_projection(params.epsg).ellipsoid
     params.orbit = compute_orbit(params)
-    params.dem = "zero_dem.tif"
+    params.dem = "expected_height_dem.tif"
     return params
 
 
-'''
-def create_burst(ref_epoch_str="2017-02-12T01:12:30.0", az_time_interval):
-    fmt = "%Y-%m-%dT%H:%M:%S.%f"
-'''
 def create_burst(params):
 
     # place holder value
@@ -173,18 +173,23 @@ def compute_expected_llh(test_params):
     xyz_vec = np.vstack([ellipsoid.a * c_lat * c_lon,
                          ellipsoid.a * c_lat * s_lon,
                          ellipsoid.a * s_lat])
-    llh = [ellipsoid.xyz_to_lon_lat(xyz_vec[:,i])
-           for i in range(test_params.n_samples)]
+    llh = np.array([ellipsoid.xyz_to_lon_lat(xyz_vec[:,i])
+                    for i in range(test_params.n_samples)])
     return llh
 
 
 def cfg_45_lat():
+    # manually populate namespace identical with attribs expected by s1_*
     cfg = SimpleNamespace()
 
+    # retrieve test-defining parameters
     cfg.test_params = common_params(45)
 
+    # compute expected LLH for (1) test comparison (2) DEM population
+    cfg.expected_llh = compute_expected_llh(cfg.test_params)
+
     # create and set DEM raster
-    zero_dem(cfg.test_params)
+    expected_height_dem(cfg)
 
     # set rdr2geo params
     rdr2geo_params = SimpleNamespace()
@@ -197,7 +202,6 @@ def cfg_45_lat():
 
     # make list with single burst
     cfg.bursts = [create_burst(cfg.test_params)]
-    import ipdb; ipdb.set_trace()
     cfg.dem = cfg.test_params.dem
     cfg.gpu_enabled = False
     cfg.gpu_id = 0
diff --git a/tests/test_rdr2geo.py b/tests/test_rdr2geo.py
index bb5f8205..b8ba6bf1 100644
--- a/tests/test_rdr2geo.py
+++ b/tests/test_rdr2geo.py
@@ -12,20 +12,23 @@ def get_rdr2geo_output(cfg):
     burst_id = burst.burst_id
     date_str = str(burst.sensing_start.date())
     topo_path = f'{cfg.product_path}/{burst_id}/{date_str}/topo.vrt'
-    print(topo_path)
     ds = gdal.Open(topo_path)
     llh = np.vstack([ds.GetRasterBand(i+1).ReadAsArray() for i in range(3)])
-    return llh
+    return llh.transpose()
 
 def test_45_lat_rdr2geo():
+    # set up test configuration for 45 degree latitude
     cfg = prep.cfg_45_lat()
 
-    expected_llh = prep.compute_expected_llh(cfg.test_params)
-
+    # run s1_rdr2geo with configuration from above
     s1_rdr2geo.run(cfg)
+
+    # retrieve and print rdr2geo output
     computed_llh = get_rdr2geo_output(cfg)
-    print(computed_llh)
-    print(expected_llh)
+    print(computed_llh, computed_llh.shape)
+
+    # print analytical solution
+    print(cfg.expected_llh, cfg.expected_llh.shape)
 
 
 if __name__ == '__main__':