Improve ASTER sensor angle calculation

* make sensor geometry calculation more efficient and more generic. Apprication is extended to ascending satellite mode (before only descending was covered).
* add sensor geometry for backward looking telescope of channel 3.

typhon/cloudmask/aster.py

@@ -486,11 +486,7 @@ def read_coordinates(self, channel="1"):
             ndarray, ndarray: latitude, longitude
         """
 
-        sdsidx = {
-            "VNIR": (0, 1),
-            "SWIR": (6, 7),
-            "TIR": (14, 15)
-        }
+        sdsidx = {"VNIR": (0, 1), "SWIR": (6, 7), "TIR": (14, 15)}
 
         latstr = ":".join(
             (
@@ -595,19 +591,22 @@ def dt_estimate_scanlines(self, sensor="vnir"):
 
         return np.linspace(-0.5, 0.5, scanlines[sensor]) * dtdelta + self.datetime
 
-    def sensor_angles(self, sensor="vnir"):
-        """Calculate a sun reflection angle for a given ASTER image depending on
-        the sensor-sun geometry and the sensor settings.
+    def sensor_angles(self, channel="1"):
+        """Calculate sensor zenith and azimuth angles.
+        
+        Angles are derived for each pixel depending on the channel and the
+        corresponding ASTER subsensor ("VNIR", "SWIR", "TIR"), as well as on the
+        subsensor geometry and settings.
 
         Note:
-            All angular values are given in [°].
+            All angular values are given in degree.
 
         Parameters:
-            sensor (str): ASTER sensor ("vnir", "swir", or "tir").
+            channel (str): ASTER channel number. '1', '2', '3N', '3B', '4',
+                    '5', '6', '7', '8', '9', '10', '11', '12', '13', '14'.
 
         Returns:
-            sensor_angle (ndarray): 2d field of size `cloudmask` of reflection
-            angles for eachannel pixel.
+            ndarray, ndarray: sensor zenith, sensor azimuth angles in degree.
 
         References:
              Kang Yang, Huaguo Zhang, Bin Fu, Gang Zheng, Weibing Guan, Aiqin Shi
@@ -617,79 +616,101 @@ def sensor_angles(self, sensor="vnir"):
              Algorithm Theoretical Basis Document for ASTER Level-1 Data 
              Processing (Ver. 3.0).1996.
         """
+        if channel != "3B":
+            sensor = self.channel2sensor[channel]
+        else:
+            sensor = "VNIRB"
+
         # Angular data from ASTER metadata data.
         S = float(self.meta["MAPORIENTATIONANGLE"])
 
-        if sensor == "vnir":
-            P = float(self.meta["POINTINGANGLE.1"])
-            FOV = 6.09  # Field of view
-            # IFOV = np.rad2deg(21.3e-6) Instantaneous field of view
-            field = self.get_radiance(channel="1")
-        elif sensor == "swir":
-            P = float(self.meta["POINTINGANGLE.2"])
-            # IFOV = np.rad2deg(42.6e-6)
-            FOV = 4.9
-            field = self.get_radiance(channel="4")
-        elif sensor == "tir":
-            P = float(self.meta["POINTINGANGLE.3"])
-            # IFOV = np.rad2deg(127.8e-6)
-            FOV = 4.9
-            field = self.get_radiance(channel="10")
-        else:
-            raise ValueError("ASTER sensors are named: vnir, swir, tir.")
-
-        # Instantaneous FOV (IVOF)
-        # IFOV = FOV / cloudmask.shape[1]
-
-        # Construct n-array indexing pixels n=- right and n=+ left from the image
-        # center and perpendicular to flight direction.
-        n = np.zeros(np.shape(field))
-        nswath = []
-        # assign indices n to pixel, neglegting the frist and last 5 swath rows
-        for i in np.arange(start=5, stop=field.shape[0] - 5, step=1, dtype=int):
-            # for all scan rows: extract first and last not-NaN index
-            ind1 = next(x[0] for x in enumerate(field[i]) if ~np.isnan(x[1]))
-            # print('Index 1: ', ind1)
-            ind2 = field.shape[1] - next(
-                x[0] for x in enumerate(field[i, 5:-5][::-1]) if ~np.isnan(x[1])
-            )
-            ind_mid = (ind1 + ind2) / 2
-            nswath.append(ind2 - ind1 + 1)
-            # create indexing arrays for right and left side of swath
-            right_arr = np.arange(start=-round(ind_mid), stop=0, step=1, dtype=int) * -1
-            left_arr = (
-                np.arange(
-                    start=1, stop=field.shape[1] - round(ind_mid) + 1, step=1, dtype=int
-                )
-                * -1
-            )
-            n[i] = np.asarray(list(right_arr) + list(left_arr))
+        FOV = {"VNIR": 6.09, "VNIRB": 5.19, "SWIR": 4.9, "TIR": 4.9}
 
-        # add the first and last 5 scan lines as dupliates of the 6th and 6th-last row.
-        n[:5] = n[5]  # first
-        n[-5:] = n[-6]  # last
-
-        n[np.isnan(field)] = np.nan
+        P = {
+            "VNIR": float(self.meta["POINTINGANGLE.1"]),
+            "VNIRB": float(self.meta["POINTINGANGLE.1"]),
+            "SWIR": float(self.meta["POINTINGANGLE.2"]),
+            "TIR": float(self.meta["POINTINGANGLE.3"]),
+        }
 
-        # instantaneous field of view (IFOV)
-        FOV_S = FOV / np.cos(np.deg2rad(S))  # correct for map orientation
-        IFOV = FOV_S / np.median(nswath)  # median scan line width
+        # cut overlap area of backward pointing telescope
+        if channel != "3B":
+            field = self.read_digitalnumbers(channel)
+        elif channel == "3B" and self.meta["FLYINGDIRECTION"] == "DE":
+            field = self.read_digitalnumbers(channel)[400:]
+        elif channel == "3B" and self.meta["FLYINGDIRECTION"] == "AE":
+            field = self.read_digitalnumbers(channel)[:400]
+
+        # design n field
+        sidx = np.arange(np.shape(field)[1])
+
+        mid0 = sidx[np.isfinite(field[5, :])][[0, -1]].mean()
+        mid1 = sidx[np.isfinite(field[-5, :])][[0, -1]].mean()
+
+        f = interpolate.interp1d(
+            np.array([5, np.shape(field)[0] - 5]),
+            np.array([mid0, mid1]),
+            kind="linear",
+            fill_value="extrapolate",
+        )
 
-        # Thresholding index for separating left and right from nadir in
-        # azimuth calculation.
-        n_az = n * IFOV + P
+        mids = f(np.arange(np.shape(field)[0]))
+        # costructing an n-array indexing the pixels symmetric to the center of the
+        # swath. If pointing angle is zero, the sensor zenith angle is zero in the
+        # swath center.
+        n = sidx - mids[:, np.newaxis]
+
+        # left and right side of nadir are defined such that the sign follows the
+        # roll angle sign, which is negative on the right and positive on the left
+        # side of the sensor in flying direction (!), NOT in projected image. The
+        # sides therefore depend on the ascending / decending mode defined in the
+        # meta data.
+        flyingdir = self.meta["FLYINGDIRECTION"]
+        if flyingdir is "DE":
+            n *= -1
+
+        swath_widths = np.sum(np.isfinite(field), axis=1)
+        # average swath width, but exluding possible NaN-scanlines at beginning and
+        # end of the image.
+        swath_width = np.mean(swath_widths[swath_widths > 4200])
+
+        n_angles = n * FOV[sensor] / swath_width + P[sensor]
+        azimuth = np.full_like(field, np.nan)
+
+        if channel != "3B":
+            zenith = abs(n_angles)
+
+            if flyingdir is "DE":
+                azimuth[n_angles > 0] = 90 + S
+                azimuth[n_angles <= 0] = 270 + S
+            else:
+                azimuth[n_angles < 0] = 90 + S
+                azimuth[n_angles >= 0] = 270 + S
+        else:
+            h = 705000  # in km above the equator
+            zenith = np.rad2deg(
+                np.arctan(
+                    np.sqrt(
+                        (h * np.tan(np.deg2rad(P[sensor])) + 15 * n) ** 2
+                        + (h * np.tan(np.deg2rad(27.6)) / np.cos(np.deg2rad(P[sensor])))
+                        ** 2
+                    )
+                    / h
+                )
+            )
 
-        # ASTER zenith angle
-        zenith = abs(np.asarray(n) * IFOV + P)
+            x = np.rad2deg(np.arctan(0.6 / np.tan(np.deg2rad(n_angles))))
+            if flyingdir is "DE":
+                azimuth[n_angles > 0] = np.array(90 - x + S)[n_angles > 0]
+                azimuth[n_angles <= 0] = np.array(270 - x + S)[n_angles <= 0]
+            else:
+                azimuth[n_angles < 0] = np.array(90 - x + S)[n_angles < 0]
+                azimuth[n_angles >= 0] = np.array(270 - x + S)[n_angles >= 0]
 
-        # ASTER azimuth angle
-        azimuth = np.ones(field.shape) * np.nan
-        # swath right side
-        azimuth[np.logical_and(n_az < 0, ~np.isnan(n))] = 90 + S
-        # swath left side
-        azimuth[np.logical_and(n_az >= 0, ~np.isnan(n))] = 270 + S
+        zenith[np.isnan(field)] = np.nan
+        azimuth[np.isnan(field)] = np.nan
 
-        return (zenith, azimuth)
+        return zenith, azimuth
 
     def reflection_angles(self):
         """Calculate the reflected sun angle, theta_r, of specular reflection
@@ -957,6 +978,7 @@ def cloudtopheight_IR(bt, cloudmask, latitude, month, method="modis"):
 
     return (bt_clear_avg - bt_cloudy) / lapserate
 
+
 def geocentric2geodetic(latitude):
     """Translate geocentric to geodetic latitudes.