Merge branch 'develop' into 'master'

Bump to version 0.4

See merge request !20

drama/geo/derived_geo.py

 import numpy as np
 import scipy.interpolate as interpolate
+import scipy.constants
 
 from drama.geo import geometry as geom
 from drama.geo.swath_geo import SingleSwath
 
 
 class BistaticRadarGeometry(object):
-    """This class calculates useful geometric relations.
+    """
+    This class calculates useful geometric relations.
 
     Parameters
     ----------
-    orbit_height :
+    orbit_height : float
         orbit height above surface
-    r_planet :
+    r_planet : float
         radious of planet, default's to Earth's
-    gm_planet :
+    gm_planet : float
         mass of planet, default's to Earths
-    degrees :
+    degrees : bool
         if set, everything is in degrees, defaults to False.
-
-    Returns
-    -------
-
     """
     def __init__(self, par_file, companion_delay=None):
         self.swth_t = SingleSwath(par_file=par_file)
@@ -67,6 +65,11 @@ class BistaticRadarGeometry(object):
         self.v_r[..., 0] = np.sum(self.swth_r.v_ecef.reshape(ashp) * self.swth_t.local_x, axis=-1)
         self.v_r[..., 1] = np.sum(self.swth_r.v_ecef.reshape(ashp) * self.swth_t.local_y, axis=-1)
         self.v_r[..., 2] = np.sum(self.swth_r.v_ecef.reshape(ashp) * self.swth_t.local_z, axis=-1)
+        # Now an additional line of sight for an antenna separated by a certain Baseline
+        r_r2 = self.r_r + self.antaxes_local[...,1]  # this adds baseline along antenna axis_
+        r_r2_n = np.linalg.norm(r_r2, axis=2)
+        self.r_v_r2 = r_r2 / r_r2_n[:,:, np.newaxis]
+
         # Magnitude (norm) of r_t and r_r
         self.r_n_t = np.linalg.norm(self.r_t, axis=2)
         self.r_n_r = np.linalg.norm(self.r_r, axis=2)
@@ -105,16 +108,6 @@ class BistaticRadarGeometry(object):
         self.set_lat(lat)
 
     def set_lat(self, lat):
-        """
-
-        Parameters
-        ----------
-        lat :
-
-        Returns
-        -------
-
-        """
         self._lat = lat
         mid_range = int(self.swth_t.lats.shape[1] / 2)
         lats = self.swth_t.lats
@@ -466,21 +459,30 @@ class BistaticRadarGeometry(object):
             return np.stack([grad_x, grad_y], axis=-1)
 
     def train_off_track_displacement(self, inc, dotr_s=1, ascending=True):
-        """ Calculates displacement of target in SAR image due to a motion induced
-        derivative of the bistatic range
+        """
+        Calculates displacement of target in SAR image due to a motion induced
+        derivative of the bistatic range.
 
         Parameters
         ----------
-        inc :
-            incident angle, in radians
-        dotr_s:
-            time derivative of bistatic range due to target motion
-        ascending :
-             (Default value = True)
+        inc : float | np.ndarray
+            Incident angle, in radians.
+        dotr_s: float
+            Time derivative of bistatic range due to target motion.
+        ascending : bool
+             If true, ascending orbit, else descending (Default value = True).
 
         Returns
         -------
-
+        tuple
+            Tuple containing the temporal offset, the offset in the azimuth
+            (cross-track) direction, and a tuple containing parameters related
+            to the derivatives. The latter contains three variables: the ratio
+            of the first derivative of the bistatic range wrt to the azimuth
+            distance to the second derivative, the derivate of the bistatic
+            range wrt to time, the 2nd derivate of the bistatic range wrt to and
+            an alternative expression for the temporal offset.
+            time
         """
         # azimuth (time) shift
         c1 = self.inc2sgrad_r_b(inc, ascending=ascending)[:, 0] / self.inc2sgrad_dot_r_b(inc, ascending=ascending)[:, 0]
@@ -498,8 +500,8 @@ class BistaticRadarGeometry(object):
     def k_ground(self, df, dta, f0=None, ascending=True):
         if f0 is None:
             f0 = self.conf.sar.f0
-        k0 = 2 * np.pi * f0 / 3e8
-        dk = df * 2 * np.pi / 3e8
+        k0 = 2 * np.pi * f0 / scipy.constants.c
+        dk = df * 2 * np.pi / scipy.constants.c
         if ascending:
             orbind = self.__asc_latind
         else:
@@ -508,13 +510,165 @@ class BistaticRadarGeometry(object):
         v_r = self.v_r[orbind]
         r_i = self.r_v_t[orbind]
         r_s = self.r_v_r[orbind]
-        dr_i = (np.sum(r_i * v_t, axis=-1)[..., np.newaxis] * r_i - v_t) / bsgeo.r_n_t[orbind, ...,np.newaxis]
-        dr_s = (np.sum(r_s * v_r, axis=-1)[..., np.newaxis] * r_s - v_r) / bsgeo.r_n_r[orbind, ...,np.newaxis]
+        dr_i = (np.sum(r_i * v_t, axis=-1)[..., np.newaxis] * r_i - v_t) / self.r_n_t[orbind, ...,np.newaxis]
+        dr_s = (np.sum(r_s * v_r, axis=-1)[..., np.newaxis] * r_s - v_r) / self.r_n_r[orbind, ...,np.newaxis]
         kv = (k0 + dk) * (r_i + r_s) + k0 * (dr_i + dr_s) * dta
         return kv
 
+    def dk_ground_dt(self, f0=None, ascending=True):
+        """
+        Derivative of bistatic wavenumber vector at the ground w.r.t. azimuth time.
+
+        Parameters
+        ----------
+        f0 : float
+            Carrier frequency, defaults to value in conf file passed to initialize function.
+        ascending : bool
+            True for ascending.
+
+        Returns
+        -------
+        np.ndarray
+            The derivative w.r.t azimuth time.
+        """
+        if f0 is None:
+            f0 = self.conf.sar.f0
+        k0 = 2 * np.pi * f0 / scipy.constants.c
+        if ascending:
+            orbind = self.__asc_latind
+        else:
+            orbind = self.__dsc_latind
+        v_t = self.v_t[orbind]
+        v_r = self.v_r[orbind]
+        r_i = self.r_v_t[orbind]
+        r_s = self.r_v_r[orbind]
+        dr_i = (np.sum(r_i * v_t, axis=-1)[..., np.newaxis] * r_i - v_t) / self.r_n_t[orbind, ...,np.newaxis]
+        dr_s = (np.sum(r_s * v_r, axis=-1)[..., np.newaxis] * r_s - v_r) / self.r_n_r[orbind, ...,np.newaxis]
+        dkvdt = k0 * (dr_i + dr_s)
+        return dkvdt
+
+    def dk_ground_df(self, ascending=True):
+        """
+        Derivative of bistatic wavenumber vector at the ground w.r.t. frequency.
+
+        Parameters
+        ----------
+        ascending : bool
+            True for ascending.
+
+        Returns
+        -------
+        np.ndarray
+            The derivative w.r.t frequency.
+        """
+        if ascending:
+            orbind = self.__asc_latind
+        else:
+            orbind = self.__dsc_latind
+        r_i = self.r_v_t[orbind]
+        r_s = self.r_v_r[orbind]
+        dkvdf = 2 * np.pi / scipy.constants.c * (r_i + r_s)
+        return dkvdf
+
+    def dk_ground_dbati(self, f0=None, ascending=True):
+        """
+        Derivative of bistatic wavenumber vector at the ground w.r.t. baseline
+        (or along-track position on receive antenna).
+
+        Parameters
+        ----------
+        f0 : float
+            Carrier frequency, defaults to value in conf file passed to initialize function.
+        ascending : bool
+            True for ascending.
+
+        Returns
+        -------
+        np.ndarray
+            The derivative w.r.t the baseline.
+        """
+        if f0 is None:
+            f0 = self.conf.sar.f0
+        k0 = 2 * np.pi * f0 / scipy.constants.c
+        if ascending:
+            orbind = self.__asc_latind
+        else:
+            orbind = self.__dsc_latind
+        #v_t = self.v_t[orbind]
+        #v_r = self.v_r[orbind]
+        #r_i = self.r_v_t[orbind]
+        r_s = self.r_v_r[orbind]
+        r_s2 = self.r_v_r2[orbind]
+        #dr_i = (np.sum(r_i * v_t, axis=-1)[..., np.newaxis] * r_i - v_t) / self.r_n_t[orbind, ...,np.newaxis]
+        #dr_s = (np.sum(r_s * v_r, axis=-1)[..., np.newaxis] * r_s - v_r) / self.r_n_r[orbind, ...,np.newaxis]
+        #dr_s2 = (np.sum(r_s2 * v_r, axis=-1)[..., np.newaxis] * r_s - v_r) / self.r_n_r[orbind, ...,np.newaxis]
+        dkvdbati = k0 * (r_s2 - r_s)
+        return dkvdbati
+
+    def bati2insarpar(self, inc, f0=None, ascending=True):
+        """
+        Computes short-baseline InSAR paramters for a 1 m baseline along antenna
+        axis.
+
+        Parameters
+        ----------
+        inc: float, np.ndarray
+            angle of incidence of incident wave, in radians
+        f0 : float
+            Carrier frequency, defaults to value in conf file passed to initialize function.
+        ascending : bool
+            True for ascending.
+
+        Returns
+        -------
+        dict
+            A dictionary with:
+            'dfdb': The frequency (spectral) shift, in Hz, per 1 m baseline
+            'dtdb': The ati lag per 1 m baseline
+            'dkzdb': The delta kz (i.e. 2\pi/h_amb) per 1 m baseline
+            'dkzdb_rar': The delta kz for the real-aperture interferometer,
+            i.e. without coregistration, which only applies if also no SAR
+            processing is done.
+        """
+        if ascending:
+            orbind = self.__asc_latind
+        else:
+            orbind = self.__dsc_latind
+        dkdt = self.dk_ground_dt(f0=f0, ascending=ascending)
+        dkdf = self.dk_ground_df(ascending=ascending)
+        dkdbati = self.dk_ground_dbati(f0=f0, ascending=ascending)
+        A = np.zeros((dkdt.shape[0],2,2))
+        A[:, :, 0] = dkdf[:,0:2]
+        A[:, :, 1] = dkdt[:,0:2]
+        invA = np.linalg.inv(A)
+        b = - dkdbati[:, 0:2]
+        df_dt = np.einsum("...ij,...j->...i", invA, b)
+        # dkz of second antenna after alignment
+        dkz = dkdbati[:,2] + dkdt[:, 2] * df_dt[:,1] + dkdf[:, 2] * df_dt[:,0]
+        inc_m = self.swth_t.master_inc[orbind, :]
+        inc2dfdb = interpolate.interp1d(inc_m, df_dt[:,0],
+                                        "quadratic",
+                                        bounds_error=False,
+                                        fill_value=np.NaN)
+        inc2dtdb = interpolate.interp1d(inc_m, df_dt[:,1],
+                                        "quadratic",
+                                        bounds_error=False,
+                                        fill_value=np.NaN)
+        inc2dkzdb = interpolate.interp1d(inc_m, dkz,
+                                         "quadratic",
+                                         bounds_error=False,
+                                         fill_value=np.NaN)
+        inc2dkz0db = interpolate.interp1d(inc_m, dkdbati[:,2],
+                                         "quadratic",
+                                         bounds_error=False,
+                                         fill_value=np.NaN)
+        return {"dfdb": inc2dfdb(inc),
+                "dtdb": inc2dtdb(inc),
+                "dkzdb":inc2dkzdb(inc),
+                "dkzdb_rar":inc2dkz0db(inc)}
+
     def inc2specsup(self, inc, mode='IWS', ascending=True, f0=None):
-        """Short summary.
+        """Calculate the spectral support from the incidence angle.
 
         Parameters
         ----------
@@ -523,7 +677,7 @@ class BistaticRadarGeometry(object):
         mode : string
             'IWS' (default).
         ascending : bool
-            Self-explainatory.
+            True for ascending orbit, False for descending.
         f0 : float
             carrier frequency to superseed value read from config file.
 
@@ -531,7 +685,6 @@ class BistaticRadarGeometry(object):
         -------
         dictionary
             Corners of spectral support region and other goodies.
-
         """
         if ascending:
             orbind = self.__asc_latind
@@ -549,7 +702,7 @@ class BistaticRadarGeometry(object):
             sw = swath[0]
         proc_bw = self.conf.IWS.proc_bw[sw]
         pulse_bw = self.conf.IWS.pulse_bw[sw]
-        dta =proc_bw/(self.inc2ddot_r_b((inc))/(3e8/f0))
+        dta =proc_bw/(self.inc2ddot_r_b((inc))/(scipy.constants.c/f0))
         v11 = self.k_ground(-pulse_bw/2, -dta/2, f0=f0, ascending=ascending)[inc_ind]
         v21 = self.k_ground(pulse_bw/2, -dta/2, f0=f0, ascending=ascending)[inc_ind]
         v12 = self.k_ground(-pulse_bw/2, dta/2, f0=f0, ascending=ascending)[inc_ind]
@@ -562,7 +715,9 @@ class BistaticRadarGeometry(object):
 
     def spectral_mask(self, inc, daz, mode='IWS', ascending=True, f0=None, drg=None,
                       baseband=False, Nx=128, Ny=128):
-        """Short summary.
+        """
+        Find the positions of the spectral support region for a set of incidence
+        angles and return the corresponding mask.
 
         Parameters
         ----------
@@ -590,9 +745,7 @@ class BistaticRadarGeometry(object):
         -------
         dictionary
             Spectral support mask and kx and ky vectors.
-
         """
-        #from matplotlib.collections import PatchCollection
         from matplotlib.path import Path
         support = self.inc2specsup(inc, mode=mode, ascending=ascending, f0=f0)
         verts = np.zeros((4 , 2))
@@ -602,7 +755,7 @@ class BistaticRadarGeometry(object):
         verts[3,:] = support['specsup'][2][0:2]
 
         if drg is None:
-            drg = 3e8/2/support['pbw']/1.1/np.sin(support['inc'])
+            drg = scipy.constants.c/2/support['pbw']/1.1/np.sin(support['inc'])
         if baseband:
             kx = 2*np.pi * np.fft.fftshift(np.fft.fftfreq(256, drg))
             ky = 2*np.pi * np.fft.fftshift(np.fft.fftfreq(256, daz))
@@ -613,7 +766,7 @@ class BistaticRadarGeometry(object):
             ky = 2*np.pi * np.fft.fftshift(np.fft.fftfreq(Ny, daz)) + support['center'][1]
         kkx, kky = np.meshgrid(kx,ky)
         poly_path=Path(verts)
-        mask = 1.0* poly_path.contains_points(np.stack((kkx.flatten(),kky.flatten()), axis=-1)).reshape(kkx.shape)
+        mask = 1.0 * poly_path.contains_points(np.stack((kkx.flatten(),kky.flatten()), axis=-1)).reshape(kkx.shape)
         for amb in [-2,-1,1,2]:
             mask_a = poly_path.contains_points(np.stack((kkx.flatten(),kky.flatten() + amb * 2*np.pi/daz), axis=-1)).reshape(kkx.shape)
             mask = mask + 1.0 * mask_a
@@ -682,3 +835,28 @@ if __name__ == '__main__':
     plt.figure()
     plt.imshow(np.fft.fftshift(np.real(np.fft.ifft2(np.fft.fft2(specmask["mask"])**2))), origin='lower')
     plt.colorbar()
+#%%
+    print(bsgeo.antaxes_sat[1000])
+    # print(bsgeo.antaxes_sat[3000])
+    print(bsgeo.antaxes_local[1500,200])
+    plt.figure()
+    plt.plot(bsgeo.antaxes_sat[:,0,1],label='x')
+    plt.plot(bsgeo.antaxes_sat[:,1,1],label='y')
+    plt.plot(bsgeo.antaxes_sat[:,2,1],label='z')
+    plt.legend()
+    plt.figure()
+    plt.plot(bsgeo.antaxes_local[3000,:,0,1],label='x')
+    plt.plot(bsgeo.antaxes_local[3000,:,1,1],label='y')
+    plt.plot(bsgeo.antaxes_local[3000,:,2,1],label='z')
+    plt.legend()
+    bsgeo.r_r.shape
+#%%
+    inc = np.linspace(30,45)
+    atipar = bsgeo.bati2insarpar(np.radians(inc))
+    plt.figure()
+    plt.plot(inc,atipar["dtdb"]*9)
+    plt.figure()
+    plt.plot(inc,atipar["dfdb"]*9)
+    plt.figure()
+    plt.plot(inc, 2*np.pi/atipar["dkzdb"]/9)
+    plt.plot(inc, 2*np.pi/atipar["dkzdb_rar"]/9,"--")

setup.cfg

 [metadata]
 name = drama
-version = 0.3
+version = 0.4
 author = Paco López-Dekker
 author_email = F.LopezDekker@tudelft.nl
 description = Delft Radar Modelling and perfornance Analysis (DRaMA)