Merge f2df9a7 into 48b6553

.github/PULL_REQUEST_TEMPLATE.md

@@ -5,6 +5,5 @@
  - [ ] Closes #xxxx <!-- remove if there is no corresponding issue, which should only be the case for minor changes -->
  - [ ] Tests added <!-- for all bug fixes or enhancements -->
  - [ ] Tests passed <!-- for all non-documentation changes) -->
- - [ ] Passes ``git diff origin/master **/*py | flake8 --diff`` <!-- remove if you did not edit any Python files -->
+ - [ ] Passes ``flake8 pyorbital`` <!-- remove if you did not edit any Python files -->
  - [ ] Fully documented <!-- remove if this change should not be visible to users, e.g., if it is an internal clean-up, or if this is part of a larger project that will be documented later -->
-1

pyorbital/geoloc.py

@@ -30,22 +30,20 @@
 # - test !!!
 
 from __future__ import print_function
-
 import numpy as np
-from numpy import cos, sin, sqrt
 
 # DIRTY STUFF. Needed the get_lonlatalt function to work on pos directly if
 # we want to print out lonlats in the end.
 from pyorbital import astronomy
-from pyorbital.orbital import *
+from pyorbital.orbital import XKMPER, F
 from pyorbital.orbital import Orbital
 
-a = 6378.137  # km
-b = 6356.75231414  # km, GRS80
-# b = 6356.752314245 # km, WGS84
+A = 6378.137  # WGS84 Equatorial radius (km)
+B = 6356.75231414  # km, GRS80
+# B = 6356.752314245 # km, WGS84
 
 
-def geodetic_lat(point, a=a, b=b):
+def geodetic_lat(point, a=A, b=B):
 
     x, y, z = point
     r = np.sqrt(x * x + y * y)
@@ -55,25 +53,23 @@ def geodetic_lat(point, a=a, b=b):
     e2 = (a * a - b * b) / (a * a)
     while True:
         phi = geod_lat
-        C = 1 / sqrt(1 - e2 * sin(phi)**2)
-        geod_lat = np.arctan2(z + a * C * e2 * sin(phi), r)
+        C = 1 / np.sqrt(1 - e2 * np.sin(phi)**2)
+        geod_lat = np.arctan2(z + a * C * e2 * np.sin(phi), r)
         if np.allclose(geod_lat, phi):
             return geod_lat
 
 
-def subpoint(query_point, a=a, b=b):
-    """Get the point on the ellipsoid under the *query_point*.
-    """
+def subpoint(query_point, a=A, b=B):
+    """Get the point on the ellipsoid under the *query_point*."""
     x, y, z = query_point
-    r = sqrt(x * x + y * y)
 
     lat = geodetic_lat(query_point)
     lon = np.arctan2(y, x)
     e2_ = (a * a - b * b) / (a * a)
-    n__ = a / sqrt(1 - e2_ * sin(lat)**2)
-    nx_ = n__ * cos(lat) * cos(lon)
-    ny_ = n__ * cos(lat) * sin(lon)
-    nz_ = (1 - e2_) * n__ * sin(lat)
+    n__ = a / np.sqrt(1 - e2_ * np.sin(lat)**2)
+    nx_ = n__ * np.cos(lat) * np.cos(lon)
+    ny_ = n__ * np.cos(lat) * np.sin(lon)
+    nz_ = (1 - e2_) * n__ * np.sin(lat)
 
     return np.stack([nx_, ny_, nz_], axis=0)
 
@@ -108,7 +104,7 @@ def vectors(self, pos, vel, roll=0.0, pitch=0.0, yaw=0.0):
         # looking down at the centre of the ellipsoid and the surface of the
         # ellipsoid. Nadir on the other hand is the point which vertical goes
         # through the satellite...
-        #nadir = -pos / vnorm(pos)
+        # nadir = -pos / vnorm(pos)
 
         nadir = subpoint(-pos)
         nadir /= vnorm(nadir)
@@ -128,8 +124,8 @@ def vectors(self, pos, vel, roll=0.0, pitch=0.0, yaw=0.0):
         return qrotate(xy_rotated, nadir, yaw)
 
     def times(self, start_of_scan):
-        #tds = [timedelta(seconds=i) for i in self._times]
-        tds = self._times.astype('timedelta64[us]')
+        # tds = [timedelta(seconds=i) for i in self._times]
+        # tds = self._times.astype('timedelta64[us]')
         try:
             return np.array(self._times) + np.datetime64(start_of_scan)
         except ValueError:
@@ -168,30 +164,28 @@ def qrotate(vector, axis, angle):
     This function uses quaternion rotation.
     """
     n_axis = axis / vnorm(axis)
-    sin_angle = np.expand_dims(sin(angle / 2), 0)
+    sin_angle = np.expand_dims(np.sin(angle / 2), 0)
     if np.ndim(n_axis) == 1:
         n_axis = np.expand_dims(n_axis, 1)
         p__ = np.dot(n_axis, sin_angle)[:, np.newaxis]
     else:
         p__ = n_axis * sin_angle
 
-    q__ = Quaternion(cos(angle / 2), p__)
+    q__ = Quaternion(np.cos(angle / 2), p__)
     shape = vector.shape
     return np.einsum("kj, ikj->ij",
                      vector.reshape((3, -1)),
                      q__.rotation_matrix()[:3, :3]).reshape(shape)
 
 
 def get_lonlatalt(pos, utc_time):
-    """Calculate sublon, sublat and altitude of satellite, considering the
-    earth an ellipsoid.
+    """Calculate sublon, sublat and altitude of satellite, considering the earth an ellipsoid.
 
     http://celestrak.com/columns/v02n03/
+
     """
     (pos_x, pos_y, pos_z) = pos / XKMPER
-    lon = ((np.arctan2(pos_y * XKMPER, pos_x * XKMPER) - astronomy.gmst(utc_time))
-           % (2 * np.pi))
-
+    lon = ((np.arctan2(pos_y * XKMPER, pos_x * XKMPER) - astronomy.gmst(utc_time)) % (2 * np.pi))
     lon = np.where(lon > np.pi, lon - np.pi * 2, lon)
     lon = np.where(lon <= -np.pi, lon + np.pi * 2, lon)
 
@@ -273,6 +267,7 @@ def hnorm(m):
     """
     return np.sqrt((m**2).sum(1))
 
+
 if __name__ == '__main__':
     # NOAA 18 (from the 2011-10-12, 16:55 utc)
     # 1 28654U 05018A   11284.35271227  .00000478  00000-0  28778-3 0  9246

pyorbital/geoloc_example.py

@@ -26,6 +26,8 @@
 import numpy as np
 from datetime import datetime
 from pyorbital.geoloc import ScanGeometry, compute_pixels, get_lonlatalt
+from mpl_toolkits.basemap import Basemap
+import matplotlib.pyplot as plt
 
 tle1 = "1 33591U 09005A   12345.45213434  .00000391  00000-0  24004-3 0  6113"
 tle2 = "2 33591 098.8821 283.2036 0013384 242.4835 117.4960 14.11432063197875"
@@ -61,19 +63,18 @@
 
 print(pos_time)
 
-
 # Plot the result
-from mpl_toolkits.basemap import Basemap
-import matplotlib.pyplot as plt
-
-m = Basemap(projection='stere', llcrnrlat=24, urcrnrlat=70, llcrnrlon=-25, urcrnrlon=120, lat_ts=58, lat_0=58, lon_0=14, resolution='l')
-
+m = Basemap(projection='stere', llcrnrlat=24, urcrnrlat=70, llcrnrlon=-25, urcrnrlon=120,
+            lat_ts=58, lat_0=58, lon_0=14, resolution='l')
 
 # convert and plot the predicted pixels in red
 x, y = m(pos_time[0], pos_time[1])
-p1 = m.plot(x,y, marker='+', color='red', markerfacecolor='red', markeredgecolor='red', markersize=1, markevery=1, zorder=4, linewidth=0.0)
+p1 = m.plot(x, y, marker='+', color='red', markerfacecolor='red', markeredgecolor='red', markersize=1, markevery=1,
+            zorder=4, linewidth=0.0)
 m.fillcontinents(color='0.85', lake_color=None, zorder=3)
-m.drawparallels(np.arange(-90.,90.,5.), labels=[1,0,1,0],fontsize=10, dashes=[1, 0], color=[0.8,0.8,0.8], zorder=1)
-m.drawmeridians(np.arange(-180.,180.,5.), labels=[0,1,0,1],fontsize=10, dashes=[1, 0], color=[0.8,0.8,0.8], zorder=2)
+m.drawparallels(np.arange(-90., 90., 5.), labels=[1, 0, 1, 0], fontsize=10, dashes=[1, 0],
+                color=[0.8, 0.8, 0.8], zorder=1)
+m.drawmeridians(np.arange(-180., 180., 5.), labels=[0, 1, 0, 1], fontsize=10, dashes=[1, 0],
+                color=[0.8, 0.8, 0.8], zorder=2)
 
 plt.show()

pyorbital/geoloc_instrument_definitions.py

@@ -263,9 +263,12 @@ def amsua_edge_geom(scans_nb):
 
 def mhs(scans_nb, edges_only=False):
     """ Describe MHS instrument geometry
+
     See:
+
     - https://www.eumetsat.int/website/home/Satellites/CurrentSatellites/Metop/MetopDesign/MHS/index.html
-    - https://www1.ncdc.noaa.gov/pub/data/satellite/publications/podguides/N-15%20thru%20N-19/pdf/0.0%20NOAA%20KLM%20Users%20Guide.pdf
+    - https://www1.ncdc.noaa.gov/pub/data/satellite/publications/podguides/
+          N-15%20thru%20N-19/pdf/0.0%20NOAA%20KLM%20Users%20Guide.pdf
       (NOAA KLM Users Guide –August 2014 Revision)
 
     Parameters:
@@ -316,10 +319,12 @@ def mhs_edge_geom(scans_nb):
 ################################################################
 
 def hirs4(scans_nb, edges_only=False):
-    """ Describe HIRS/4 instrument geometry
+    """Describe HIRS/4 instrument geometry.
+
     See:
     - https://www.eumetsat.int/website/home/Satellites/CurrentSatellites/Metop/MetopDesign/HIRS/index.html
-    - https://www1.ncdc.noaa.gov/pub/data/satellite/publications/podguides/N-15%20thru%20N-19/pdf/0.0%20NOAA%20KLM%20Users%20Guide.pdf
+    - https://www1.ncdc.noaa.gov/pub/data/satellite/publications/podguides/
+          N-15%20thru%20N-19/pdf/0.0%20NOAA%20KLM%20Users%20Guide.pdf
       (NOAA KLM Users Guide –August 2014 Revision)
 
     Parameters:
@@ -372,8 +377,12 @@ def hirs4_edge_geom(scans_nb):
 def atms(scans_nb, edges_only=False):
     """ Describe MHS instrument geometry
     See:
-    https://dtcenter.org/com-GSI/users/docs/presentations/2013_workshop/Garrett_GSI_2013.pdf (Assimilation of Suomi-NPP ATMS, Kevin Garrett et al., August 8, 2013)
-    https://www.star.nesdis.noaa.gov/star/documents/meetings/2016JPSSAnnual/S4/S4_13_JPSSScience2016_session4Part2_ATMS_Scan_Reversal_HYANG.pdf (Suomi NPP ATMS Scan Reversal Study, Hu (Tiger) Yang, NOAA/STAR ATMS SDR Working Group)
+
+    - https://dtcenter.org/com-GSI/users/docs/presentations/2013_workshop/
+          Garrett_GSI_2013.pdf (Assimilation of Suomi-NPP ATMS, Kevin Garrett et al., August 8, 2013)
+    - https://www.star.nesdis.noaa.gov/star/documents/meetings/2016JPSSAnnual/
+          S4/S4_13_JPSSScience2016_session4Part2_ATMS_Scan_Reversal_HYANG.pdf
+          (Suomi NPP ATMS Scan Reversal Study, Hu (Tiger) Yang, NOAA/STAR ATMS SDR Working Group)
 
     Parameters:
        scans_nb | int -  number of scan lines

pyorbital/orbital.py

@@ -66,7 +66,7 @@
 XKE = 0.743669161e-1
 XKMPER = 6378.135
 XMNPDA = 1440.0
-#MFACTOR = 7.292115E-5
+# MFACTOR = 7.292115E-5
 AE = 1.0
 SECDAY = 8.6400E4
 
@@ -194,7 +194,7 @@ def get_last_an_time(self, utc_time):
 
         # Bisect to z within 1 km
         while np.abs(pos1[2]) > 1:
-            pos0, vel0 = pos1, vel1
+            # pos0, vel0 = pos1, vel1
             dt = (t_old - t_new) / 2
             t_mid = t_old - dt
             pos1, vel1 = self.get_position(t_mid, normalize=False)
@@ -463,15 +463,15 @@ def fprime(timex):
         tx0 = utc_time - timedelta(seconds=1.0)
         tx1 = utc_time
         idx = 0
-        #eps = 500.
+        # eps = 500.
         eps = 100.
         while abs(tx1 - tx0) > precision and idx < nmax_iter:
             tx0 = tx1
             fpr = fprime(tx0)
             # When the elevation is high the scale is high, and when
             # the elevation is low the scale is low
-            #var_scale = np.abs(np.sin(fpr[0] * np.pi/180.))
-            #var_scale = np.sqrt(var_scale)
+            # var_scale = np.abs(np.sin(fpr[0] * np.pi/180.))
+            # var_scale = np.sqrt(var_scale)
             var_scale = np.abs(fpr[0])
             tx1 = tx0 - timedelta(seconds=(eps * var_scale * fpr[1]))
             idx = idx + 1
@@ -547,20 +547,20 @@ class _SGDP4(object):
     def __init__(self, orbit_elements):
         self.mode = None
 
-        perigee = orbit_elements.perigee
+        # perigee = orbit_elements.perigee
         self.eo = orbit_elements.excentricity
         self.xincl = orbit_elements.inclination
         self.xno = orbit_elements.original_mean_motion
-        k_2 = CK2
-        k_4 = CK4
-        k_e = XKE
+        # k_2 = CK2
+        # k_4 = CK4
+        # k_e = XKE
         self.bstar = orbit_elements.bstar
         self.omegao = orbit_elements.arg_perigee
         self.xmo = orbit_elements.mean_anomaly
         self.xnodeo = orbit_elements.right_ascension
         self.t_0 = orbit_elements.epoch
         self.xn_0 = orbit_elements.mean_motion
-        A30 = -XJ3 * AE**3
+        # A30 = -XJ3 * AE**3
 
         if not(0 < self.eo < ECC_LIMIT_HIGH):
             raise OrbitalError('Eccentricity out of range: %e' % self.eo)
@@ -632,12 +632,11 @@ def __init__(self, orbit_elements):
 
         self.c1 = self.bstar * self.c2
 
-        self.c4 = (2.0 * self.xnodp * coef_1 * self.aodp * betao2 * (self.eta *
-                                                                     (2.0 + 0.5 * etasq) + self.eo * (0.5 + 2.0 *
-                                                                                                      etasq) - (2.0 * CK2) * tsi / (self.aodp * psisq) * (-3.0 *
-                                                                                                                                                          self.x3thm1 * (1.0 - 2.0 * eeta + etasq *
-                                                                                                                                                                         (1.5 - 0.5 * eeta)) + 0.75 * self.x1mth2 * (2.0 *
-                                                                                                                                                                                                                     etasq - eeta * (1.0 + etasq)) * np.cos(2.0 * self.omegao))))
+        self.c4 = (2.0 * self.xnodp * coef_1 * self.aodp * betao2 * (
+            self.eta * (2.0 + 0.5 * etasq) + self.eo * (0.5 + 2.0 * etasq) - (2.0 * CK2) * tsi /
+            (self.aodp * psisq) * (-3.0 * self.x3thm1 * (1.0 - 2.0 * eeta + etasq * (1.5 - 0.5 * eeta)) +
+                                   0.75 * self.x1mth2 * (2.0 * etasq - eeta * (1.0 + etasq)) *
+                                   np.cos(2.0 * self.omegao))))
 
         self.c5, self.c3, self.omgcof = 0.0, 0.0, 0.0
 
@@ -710,7 +709,7 @@ def propagate(self, utc_time):
         kep = {}
 
         # get the time delta in minutes
-        #ts = astronomy._days(utc_time - self.t_0) * XMNPDA
+        # ts = astronomy._days(utc_time - self.t_0) * XMNPDA
         # print utc_time.shape
         # print self.t_0
         utc_time = dt2np(utc_time)
@@ -878,6 +877,7 @@ def kep2xyz(kep):
 
     return np.array((x, y, z)), np.array((v_x, v_y, v_z))
 
+
 if __name__ == "__main__":
     obs_lon, obs_lat = np.deg2rad((12.4143, 55.9065))
     obs_alt = 0.02

pyorbital/tests/__init__.py

@@ -19,28 +19,26 @@
 
 # You should have received a copy of the GNU General Public License
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
-
-"""The tests package.
-"""
+"""The tests package."""
 
 from pyorbital.tests import (test_aiaa, test_tlefile, test_orbital,
                              test_astronomy, test_geoloc)
 import unittest
 
+
 def suite():
-    """The global test suite.
-    """
+    """The global test suite."""
     mysuite = unittest.TestSuite()
     # Test the documentation strings
-    #mysuite.addTests(doctest.DocTestSuite(image))
+    # mysuite.addTests(doctest.DocTestSuite(image))
     # Use the unittests also
     mysuite.addTests(test_aiaa.suite())
     mysuite.addTests(test_tlefile.suite())
     mysuite.addTests(test_orbital.suite())
     mysuite.addTests(test_astronomy.suite())
     mysuite.addTests(test_geoloc.suite())
-
     return mysuite
 
+
 if __name__ == '__main__':
     unittest.TextTestRunner(verbosity=2).run(suite())

pyorbital/tests/test_aiaa.py

@@ -28,7 +28,7 @@
 
 import os
 import unittest
-from datetime import datetime, timedelta
+from datetime import datetime
 
 import numpy as np
 
@@ -107,10 +107,10 @@ def test_aiaa(self):
 
                     try:
                         o = LineOrbital("unknown", line1, line2)
-                    except NotImplementedError as e:
+                    except NotImplementedError:
                         test_line = f__.readline()
                         continue
-                    except ChecksumError as e:
+                    except ChecksumError:
                         self.assertTrue(test_line.split()[1] in [
                                         "33333", "33334", "33335"])
                     for delay in times: