Focus using external ephemeris (#948)

* Read orbit and attitude data files if provided.

* Add python interface to update Attitude reference epoch.

* Add python interface to update Orbit reference epoch.

* Bind copy constructors for Orbit and Attitude

* Update SLC attitude writer (epoch added in b119fa2)

* Use a consistent time epoch for all objects.

* Use external orbit and attitude in focus.py unit test.

* Remove unnecessary raw string delimiters.

Co-authored-by: Geoffrey M Gunter <[email protected]>

* Add copy dunder methods.

* Add contains() method to Attitude, Orbit, and their python bindings.

* Require orbit covers raw data interval and load earlier so available for output grid init.

* Improve default output grid.

* Do the reskew calculations in a rigorous way and kludge short orbits.

* Move useful junk into isce3 API

* Remove superfluous approx reskew calculation.

* Remove extraneous return parameters.

* Note assumption

* Remove superfluous raw string delimiters.

Co-authored-by: Geoffrey M Gunter <[email protected]>

* Remove shebang

Co-authored-by: Geoffrey M Gunter <[email protected]>

* Remove sh

Co-authored-by: Geoffrey M Gunter <[email protected]>

* Note optional arguments.

Co-authored-by: Geoffrey M Gunter <[email protected]>

* Add blank line

Co-authored-by: Geoffrey M Gunter <[email protected]>

* Add blank line

Co-authored-by: Geoffrey M Gunter <[email protected]>

* Use Orbit.contains method in test.

Co-authored-by: Geoffrey M Gunter <[email protected]>

* Only copy orbit and attitude if necessary.

Co-authored-by: Geoffrey M Gunter <[email protected]>

* Go Vader on ob1

Co-authored-by: Geoffrey M Gunter <[email protected]>

* Add missing newline

* use typing.Optional annotations

* Add type annotations to sar_duration

* Use clearer units for chirp duration.

* More appropriate log levels for grid truncation messages.

* Maintain identity in trivial rdr2rdr case.

* Fix up commit d6c74a8

* Prefer builtin round()

Co-authored-by: Geoffrey M Gunter <[email protected]>

Author: 2 people

cxx/isce3/core/Attitude.h

@@ -59,6 +59,11 @@ class Attitude {
     /** Time of last measurement relative to reference epoch (s) */
     double endTime() const { return _time[size() - 1]; }
 
+    /** Check if time falls in the valid interpolation domain. */
+    bool contains(double time) const {
+        return (startTime() <= time) && (time <= endTime());
+    }
+
     /** UTC time of first measurement */
     DateTime startDateTime() const
     {

cxx/isce3/core/Orbit.h

@@ -76,7 +76,12 @@ class Orbit {
     /** Reference epoch (UTC) */
     const DateTime & referenceEpoch() const { return _reference_epoch; }
 
-    /** Set reference epoch (UTC) */
+    /** Set reference epoch (UTC)
+     *
+     * Also updates relative time tags so that
+     *      referenceEpoch() + TimeDelta(time()[i])
+     * results in the same absolute time tags before and after this call.
+     */
     void referenceEpoch(const DateTime &);
 
     /** Interpolation method */
@@ -103,6 +108,11 @@ class Orbit {
     /** UTC time of last state vector */
     DateTime endDateTime() const { return _reference_epoch + TimeDelta(endTime()); }
 
+    /** Check if time falls in the valid interpolation domain. */
+    bool contains(double time) const {
+        return (startTime() <= time) && (time <= endTime());
+    }
+
     /** Time interval between state vectors (s) */
     double spacing() const { return _time.spacing(); }
 

python/extensions/pybind_isce3/core/Attitude.cpp

@@ -35,6 +35,22 @@ void addbinding(pybind11::class_<Attitude>& pyAttitude)
             "Serialize Attitude to h5py.Group object.",
             py::arg("h5py_group"))
 
+        // add a function rather than just a settable property since
+        // self.time is also modified
+        .def("update_reference_epoch",
+            py::overload_cast<const DateTime&>(&Attitude::referenceEpoch),
+            R"(
+            Set reference epoch and modify time tags so that
+                self.reference_epoch + TimeDelta(self.time[i])
+            remains invariant for all i in range(len(self.time)).)",
+            py::arg("reference_epoch"))
+
+        .def("copy", [](const Attitude& self) { return Attitude(self);})
+        .def("__copy__", [](const Attitude& self) { return Attitude(self);})
+        .def("__deepcopy__", [](const Attitude& self, py::dict) {
+                return Attitude(self);
+            }, py::arg("memo"))
+
         .def_property_readonly("size", &Attitude::size)
         .def_property_readonly("time", &Attitude::time)
         .def_property_readonly("quaternions", &Attitude::quaternions)
@@ -44,5 +60,8 @@ void addbinding(pybind11::class_<Attitude>& pyAttitude)
         .def_property_readonly("end_datetime", &Attitude::endDateTime)
         .def_property_readonly("reference_epoch",
                 py::overload_cast<>(&Attitude::referenceEpoch, py::const_))
+        .def("contains", &Attitude::contains,
+            "Check if time falls in the valid interpolation domain.",
+            py::arg("time"))
         ;
 }

python/extensions/pybind_isce3/core/Orbit.cpp

@@ -34,6 +34,20 @@ void addbinding(py::class_<Orbit> & pyOrbit)
         .def(py::init<std::vector<StateVector>, isce3::core::DateTime>())
         .def_property_readonly("reference_epoch",
                 py::overload_cast<>(&Orbit::referenceEpoch, py::const_))
+        // add a function rather than just a settable property since
+        // self.time is also modified
+        .def("update_reference_epoch",
+            py::overload_cast<const DateTime&>(&Orbit::referenceEpoch),
+            R"(
+            Set reference epoch and modify time tags so that
+                self.reference_epoch + TimeDelta(self.time[i])
+            remains invariant for all i in range(len(self.time)).)",
+            py::arg("reference_epoch"))
+        .def("copy", [](const Orbit& self) { return Orbit(self); })
+        .def("__copy__", [](const Orbit& self) { return Orbit(self); })
+        .def("__deepcopy__", [](const Orbit& self,  py::dict) {
+                return Orbit(self);
+            }, py::arg("memo"))
         .def_property_readonly("time", py::overload_cast<>(&Orbit::time, py::const_))
         .def_property_readonly("position", [](const Orbit & self) {
             return py::roarray(toBuffer(self.position()), py::cast(self));
@@ -79,5 +93,9 @@ void addbinding(py::class_<Orbit> & pyOrbit)
             },
             "Interpolate platform position and velocity",
             py::arg("t"))
+
+        .def("contains", &Orbit::contains,
+            "Check if time falls in the valid interpolation domain.",
+            py::arg("time"))
         ;
 }

python/packages/isce3/focus/__init__.py

@@ -1 +1,2 @@
 from pybind_isce3.focus import *
+from .sar_duration import get_sar_duration

python/packages/isce3/focus/sar_duration.py

@@ -0,0 +1,36 @@
+import isce3
+from isce3.core import Orbit, Ellipsoid
+import numpy as np
+
+
+def get_sar_duration(t: float, r: float, orbit: Orbit, ellipsoid: Ellipsoid,
+                     azres: float, wavelength: float):
+    """
+    Get approximate synthetic aperture duration (coherent processing interval).
+
+    Parameters
+    ---------
+    t : float
+        Azimuth time, in seconds since orbit epoch
+    r : float
+        Range, in meters
+    orbit : isce3.core.Orbit
+    ellipsoid : isce3.core.Ellipsoid
+    azres : float
+        Desired azimuth resolution, in meters
+
+    Returns
+    -------
+    cpi : float
+        Synthetic aperture duration (approximate), s
+    """
+    pos, vel = orbit.interpolate(t)
+    vs = np.linalg.norm(vel)
+    # assume radar and target are close enough that r_dir() is constant
+    lon, lat, h = ellipsoid.xyz_to_lon_lat(pos)
+    hdg = isce3.geometry.heading(lon, lat, vel)
+    a = ellipsoid.r_dir(hdg, lat)
+    # assume circular orbit (no nadir component of velocity)
+    vg = vs / (1 + h / a)
+    # small angle approximations and ignore window factor
+    return wavelength / (2 * azres) * r / vg

python/packages/isce3/geometry/__init__.py

@@ -1 +1,2 @@
 from pybind_isce3.geometry import *
+from .rdr2rdr import rdr2rdr

python/packages/isce3/geometry/rdr2rdr.py

@@ -0,0 +1,72 @@
+from isce3.core import Ellipsoid, LookSide, LUT2d, Orbit
+from isce3.geometry import DEMInterpolator, rdr2geo, geo2rdr
+from typing import Union, Optional
+
+
+def rdr2rdr(t: float,
+            r: float,
+            orbit: Orbit,
+            side: Union[str, LookSide],
+            doppler: LUT2d,
+            wavelength: float,
+            dem: DEMInterpolator = DEMInterpolator(),
+            ellipsoid: Ellipsoid = Ellipsoid(),
+            doppler_out: Optional[LUT2d] = None,
+            orbit_out: Optional[Orbit] = None,
+            rdr2geo_params=dict(),
+            geo2rdr_params=dict()):
+    """
+    Convert coordinates from one radar geometry to another with a different
+    Doppler (reskew) or orbit (motion compensation).
+
+    Parameters
+    ----------
+    t : float
+	    Azimuth time, seconds past orbit epoch
+    r : float
+	    Range, m
+    orbit : isce3.core.Orbit
+	    Orbit defining radar motion on input path
+    side : {"left", "right", isce3.core.LookSide.Left, isce3.core.LookSide.Right}
+	    Flag desribing which side the radar is looking.
+    doppler : isce3.core.LUT2d
+        Doppler look up table vs input range and azimuth time, Hz
+    wavelength : float
+        Wavelength associated with Doppler, m
+    dem : isce3.geometry.DEMInterpolator, optional
+        Digital elevation model, m above ellipsoid.  Defaults to h=0.
+    ellipsoid : isce3.core.Ellipsoid, optional
+        Ellipsoid describing surface.  Defaults to WGS-84.
+    doppler_out : isce3.core.LUT2d, optional
+        Doppler look up table vs output range and azimuth time, Hz
+        Defaults to input `doppler`.
+    orbit_out : isce3.core.Orbit, optional
+        Orbit defining radar motion on output path.  Defaults to input `orbit`.
+    rdr2geo_params : dict, optional
+        Dictionary specifying convergence paramters for rdr2geo solver.
+        Keys among {"threshold", "maxiter", "extraiter"}
+        See isce3.geometry.rdr2geo
+    geo2rdr_params: dict, optional
+        Dictionary specifying convergence paramters for geo2rdr solver.
+        Keys among {"threshold", "maxiter", "delta_range"}
+        See isce3.geometry.geo2rdr
+
+    Returns
+    -------
+    t_out : float
+        Azimuth time in output geometry, seconds past orbit epoch
+    r_out : float
+        Range in output geometry, m
+    """
+    if (orbit_out is None) and (doppler_out is None):
+        return t, r
+    if orbit_out is None:
+        orbit_out = orbit
+    if doppler_out is None:
+        doppler_out = doppler
+    doppler_in = doppler.eval(t, r)
+    llh = rdr2geo(t, r, orbit, side, doppler_in, wavelength, dem, ellipsoid,
+                  **rdr2geo_params)
+    tout, rout = geo2rdr(llh, ellipsoid, orbit_out, doppler_out, wavelength,
+                         side, **geo2rdr_params)
+    return tout, rout

python/packages/nisar/products/writers/SLC.py

@@ -268,7 +268,7 @@ def set_orbit(self, orbit: Orbit, accel=None, type="Custom"):
             " WGS84 G1762 reference frame")
         d.attrs["units"] = np.string_("meters per second squared")
 
-    def set_attitude(self, attitude: Attitude, epoch: DateTime, type="Custom"):
+    def set_attitude(self, attitude: Attitude, type="Custom"):
         log.info("Writing attitude to SLC")
         g = self.root.require_group("metadata/attitude")
         d = g.require_dataset("attitudeType", (), "S10", data=np.string_(type))
@@ -279,7 +279,7 @@ def set_attitude(self, attitude: Attitude, epoch: DateTime, type="Custom"):
         d.attrs["description"] = np.string_("Time vector record. This record"
             " contains the time corresponding to attitude and quaternion"
             " records")
-        d.attrs["units"] = np.string_(time_units(epoch))
+        d.attrs["units"] = np.string_(time_units(attitude.reference_epoch))
         # TODO attitude rates
         n = len(attitude.time)
         qdot = np.zeros((n, 3))