.bumpversion.cfg

 [bumpversion]
-current_version = 1.1.1
+current_version = 1.2.0
 commit = True
 tag = True
 

.github/ISSUE_TEMPLATE.md

+#### Code Sample, a minimal, complete, and verifiable piece of code
+
+```python
+# Your code here
+
+```
+#### Problem description
+
+[this should also explain **why** the current behaviour is a problem and why the 
+expected output is a better solution.]
+
+#### Expected Output
+
+#### Actual Result, Traceback if applicable
+
+#### Versions of Python, package at hand and relevant dependencies
+
+Thank you for reporting an issue !

.github/PULL_REQUEST_TEMPLATE.md

+<!-- Please make the PR against the `develop` branch. -->
+
+<!-- Describe what your PR does, and why -->
+
+ - [ ] 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/develop **/*py | flake8 --diff`` <!-- 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 -->

changelog.rst

 Changelog
 =========
 
+v1.2.0 (2018-02-13)
+-------------------
+
+- Update changelog. [Adam.Dybbroe]
+
+- Bump version: 1.1.1 → 1.2.0. [Adam.Dybbroe]
+
+- Merge branch 'develop' into new_release. [Adam.Dybbroe]
+
+- Add for NOAA-20. [Adam.Dybbroe]
+
+  Signed-off-by: Adam.Dybbroe <adam.dybbroe@smhi.se>
+
+
+- Add github issue and PR templates. [Adam.Dybbroe]
+
+  Signed-off-by: Adam.Dybbroe <adam.dybbroe@smhi.se>
+
+
+- Provide units in the get_observer_look documentation. [Martin Raspaud]
+
+  Fixes #17
+
+- Check that there are required amount of items in the row. [Panu
+  Lahtinen]
+
+- Improve documentation: Emphasize the importance of fresh/actual TLEs.
+  [Adam.Dybbroe]
+
+  Signed-off-by: Adam.Dybbroe <adam.dybbroe@smhi.se>
+
+
+- Add GOES-16 and Himawari-9. [Adam.Dybbroe]
+
+  Signed-off-by: Adam.Dybbroe <adam.dybbroe@smhi.se>
+
+
+- Merge pull request #16 from jordanlui/newTLEurl. [Martin Raspaud]
+
+  New TLE URLs
+
+- Additional TLE URLs from CelesTrak added. [Jordan Lui]
+
+- Added additional TLE paths to TLE_URL. [Jordan Lui]
+
+- Replace numpy.rank with numpy.ndim to make compatible with future
+  numpy versions. [Adam.Dybbroe]
+
+  Currently at 1.13.0 a VisibleDeprecationWarning is raised using numpy.rank
+
+  Signed-off-by: Adam.Dybbroe <a000680@c20671.ad.smhi.se>
+
+
+- Numpy 1.12 compatible - accepting number of scans to be a float.
+  [Adam.Dybbroe]
+
+  Signed-off-by: Adam.Dybbroe <a000680@c20671.ad.smhi.se>
+
+
+- Fix avhrr_gac instrument definition. [Martin Raspaud]
+
+- Fix conversion to datetime64. [Martin Raspaud]
+
+- Fix geoloc tests. [Martin Raspaud]
+
+- Add support for 2d time arrays in geoloc. [Martin Raspaud]
+
+- Merge branch 'develop' of github.com:pytroll/pyorbital into develop.
+  [Martin Raspaud]
+
+  Conflicts:
+  	pyorbital/tests/test_aiaa.py
+
+
+- Merge pull request #13 from pytroll/feature_np_datetime64. [Martin
+  Raspaud]
+
+  Feature np datetime64
+
+- Merge branch 'develop' into feature_np_datetime64. [Martin Raspaud]
+
+- Finish conversion to np datetime64. [Martin Raspaud]
+
+- Add conversion func between datetime and np.datetime64. [Martin
+  Raspaud]
+
+- Do not crash when start_of_scan already is datetime64. [Martin
+  Raspaud]
+
+  Signed-off-by: Martin Raspaud <martin.raspaud@smhi.se>
+
+
+- Convert input times to datetime64. [Martin Raspaud]
+
+  Signed-off-by: Martin Raspaud <martin.raspaud@smhi.se>
+
+
+- Allow offset application to be turned off (avhrr) [Martin Raspaud]
+
+  Signed-off-by: Martin Raspaud <martin.raspaud@smhi.se>
+
+
+- Adapt to datetime64. [Martin Raspaud]
+
+  Signed-off-by: Martin Raspaud <martin.raspaud@smhi.se>
+
+
+- Cleanup style. [Martin Raspaud]
+
+- Fix indexing. [Martin Raspaud]
+
+- Merge pull request #19 from howff/patch-1. [Adam Dybbroe]
+
+  Update platforms.txt with NOAA-20, MSG 4, GOES-16
+
+- Update platforms.txt with NOAA-20, MSG 4, GOES-16. [howff]
+
+- Merge pull request #14 from kconkas/master. [Martin Raspaud]
+
+  Python3 fixes for fetch()
+
+- Python3 fixes for fetch() [Kristijan Conkas]
+
 v1.1.1 (2017-01-10)
 -------------------
 

debian/changelog

-pyorbital (1.1.1-2) UNRELEASED; urgency=medium
+pyorbital (1.2.0-1) unstable; urgency=medium
 
+  [ Bas Couwenberg ]
   * Team upload.
   * Add python3-sphinx to build dependencies.
   * Strip trailing whitespace from changelog.
 
- -- Bas Couwenberg <sebastic@debian.org>  Wed, 27 Sep 2017 15:00:00 +0200
+  [ Antonio Valentino ]
+  * New upstream release
+  * Standard version bumped to 4.1.3 (no change)
+  * debian/control:
+    - remove un-necessary Testsuite field
+    - declare the *-doc package Multi-Arch foreign
+  * Set compat to 11
+  * debian/rules
+    - remove unnecessary dh argument: --parallel
+
+ -- Antonio Valentino <antonio.valentino@tiscali.it>  Sun, 18 Feb 2018 09:58:31 +0000
 
 pyorbital (1.1.1-1) unstable; urgency=medium
 

debian/compat

-9
+11

debian/control

@@ -2,9 +2,8 @@ Source: pyorbital
 Maintainer: Debian GIS Project <pkg-grass-devel@lists.alioth.debian.org>
 Uploaders: Antonio Valentino <antonio.valentino@tiscali.it>
 Section: python
-Testsuite: autopkgtest
 Priority: optional
-Build-Depends: debhelper (>= 9.0.0),
+Build-Depends: debhelper (>= 11),
                dh-python,
                python-setuptools,
                python3-setuptools,
@@ -14,7 +13,7 @@ Build-Depends: debhelper (>= 9.0.0),
                python3-numpy,
                python-sphinx,
                python3-sphinx
-Standards-Version: 3.9.8
+Standards-Version: 4.1.3
 Vcs-Browser: https://anonscm.debian.org/cgit/pkg-grass/pyorbital.git
 Vcs-Git: https://anonscm.debian.org/git/pkg-grass/pyorbital.git
 Homepage: https://github.com/mraspaud/pyorbital
@@ -53,6 +52,7 @@ Description: Orbital and astronomy computations in Python 3
 
 Package: python-pyorbital-doc
 Architecture: all
+Multi-Arch: foreign
 Section: doc
 Depends: ${misc:Depends},
          ${sphinxdoc:Depends}

debian/rules

@@ -9,7 +9,7 @@ export PYBUILD_NAME=pyorbital
 
 
 %:
-	dh $@ --parallel --with python2,python3,sphinxdoc --buildsystem=pybuild
+	dh $@ --with python2,python3,sphinxdoc --buildsystem=pybuild
 
 
 override_dh_auto_build:

doc/source/index.rst

@@ -27,15 +27,43 @@ The orbital module enables computation of satellite position and velocity at a s
 
     >>> from pyorbital.orbital import Orbital
     >>> from datetime import datetime
-    >>> orb = Orbital("noaa 18")
+    >>> # Use current TLEs from the internet:
+    >>> orb = Orbital("Suomi NPP")
     >>> now = datetime.utcnow()
     >>> # Get normalized position and velocity of the satellite:
     >>> orb.get_position(now)
-    ([0.57529384846822862, 0.77384005228105424, 0.59301408257897559],
-    [0.031846489698768146, 0.021287993461926374, -0.05854106186659274])
+    (array([-0.20015267,  0.09001458,  1.10686756]),
+     array([ 0.06148495,  0.03234914,  0.00846805]))
     >>> # Get longitude, latitude and altitude of the satellite:
     >>> orb.get_lonlatalt(now)
-    (-1.1625895579622014, 0.55402132517640568, 847.89381184656702)
+    (40.374855865574951, 78.849923885700363, 839.62504115338368)
+ 
+
+Use actual TLEs to increase accuracy
+------------------------------------
+
+    >>> from pyorbital.orbital import Orbital
+    >>> from datetime import datetime
+    >>> orb = Orbital("Suomi NPP")
+    >>> dtobj = datetime(2015,2,7,3,0)
+    >>> orb.get_lonlatalt(dtobj)
+    (152.11564698762811, 20.475251739329622, 829.37355785502211)
+
+But since we are interesting knowing the position of the Suomi-NPP more than
+two and half years from now (September 26, 2017) we can not rely on the current
+TLEs, but rather need a TLE closer to the time of interest:
+
+    >>> snpp = Orbital('Suomi NPP', tle_file='/data/lang/satellit/polar/orbital_elements/TLE/201502/tle-20150207.txt')
+    >>> snpp.get_lonlatalt(dtobj)
+    (105.37373804512762, 79.160752404540133, 838.94605490133154)
+
+If we take a TLE from one week earlier we get a slightly different result:
+
+    >>> snpp = Orbital('Suomi NPP', tle_file='/data/lang/satellit/polar/orbital_elements/TLE/201501/tle-20150131.txt')
+    >>> snpp.get_lonlatalt(dtobj)
+    (104.1539184988462, 79.328272480878141, 838.81555967963391)
+
+
 
 Computing astronomical parameters
 ---------------------------------

etc/platforms.txt

@@ -25,12 +25,15 @@ FY-2G 40367
 FY-3A 32958
 FY-3B 37214
 FY-3C 39260
+FY-3D 43010
 GOES-13 29155
 GOES-14 35491
 GOES-15 36411
+GOES-16 41866
 Himawari-6 28622
 Himawari-7 28937
 Himawari-8 40267
+Himawari-9 41836
 INSAT-3A 27714
 INSAT-3C 27298
 INSAT-3D 39216
@@ -42,6 +45,7 @@ Meteosat-7 24932
 Meteosat-8 27509
 Meteosat-9 28912
 Meteosat-10 38552
+Meteosat-11 40732
 Metop-A 29499
 Metop-B 38771
 NOAA-15 25338
@@ -49,6 +53,7 @@ NOAA-16 26536
 NOAA-17 27453
 NOAA-18 28654
 NOAA-19 33591
+NOAA-20 43013
 RadarSat-2 32382
 Sentinel-1A 39634
 SMOS 36036

pyorbital/__init__.py

+# -*- coding: utf-8 -*-
+
+# Copyright (c) 2017
+
+# Author(s):
+
+#   Martin Raspaud <martin.raspaud@smhi.se>
+
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+import numpy as np
+
+
+def dt2np(utc_time):
+    try:
+        return np.datetime64(utc_time)
+    except ValueError:
+        return utc_time.astype('datetime64[ns]')

pyorbital/astronomy.py

@@ -24,18 +24,19 @@
 Parts taken from http://www.geoastro.de/elevaz/basics/index.htm
 """
 
-import datetime
 import numpy as np
 
+from pyorbital import dt2np
 
 F = 1 / 298.257223563  # Earth flattening WGS-84
 A = 6378.137  # WGS84 Equatorial radius
 MFACTOR = 7.292115E-5
 
+
 def jdays2000(utc_time):
     """Get the days since year 2000.
     """
-    return _days(utc_time - datetime.datetime(2000, 1, 1, 12, 0))
+    return _days(dt2np(utc_time) - np.datetime64('2000-01-01T12:00'))
 
 
 def jdays(utc_time):
@@ -43,22 +44,11 @@ def jdays(utc_time):
     """
     return jdays2000(utc_time) + 2451545
 
-def _fdays(dt):
-    """Get the days (floating point) from *d_t*.
-    """
-    return (dt.days +
-            (dt.seconds +
-             dt.microseconds / (1000000.0)) / (24 * 3600.0))
-
-_vdays = np.vectorize(_fdays)
 
 def _days(dt):
     """Get the days (floating point) from *d_t*.
     """
-    try:
-        return _fdays(dt)
-    except AttributeError:
-        return _vdays(dt)
+    return dt / np.timedelta64(1, 'D')
 
 
 def gmst(utc_time):
@@ -97,6 +87,7 @@ def sun_ecliptic_longitude(utc_time):
     l__ = l_0 + d_l
     return np.deg2rad(l__)
 
+
 def sun_ra_dec(utc_time):
     """Right ascension and declination of the sun at *utc_time*.
     """
@@ -115,6 +106,7 @@ def sun_ra_dec(utc_time):
     right_ascension = 2 * np.arctan2(y__, (x__ + r__))
     return right_ascension, declination
 
+
 def _local_hour_angle(utc_time, longitude, right_ascension):
     """Hour angle at *utc_time* for the given *longitude* and
     *right_ascension*
@@ -122,6 +114,7 @@ def _local_hour_angle(utc_time, longitude, right_ascension):
     """
     return _lmst(utc_time, longitude) - right_ascension
 
+
 def get_alt_az(utc_time, lon, lat):
     """Return sun altitude and azimuth from *utc_time*, *lon*, and *lat*.
     lon,lat in degrees
@@ -137,6 +130,7 @@ def get_alt_az(utc_time, lon, lat):
             np.arctan2(-np.sin(h__), (np.cos(lat) * np.tan(dec) -
                                       np.sin(lat) * np.cos(h__))))
 
+
 def cos_zen(utc_time, lon, lat):
     """Cosine of the sun-zenith angle for *lon*, *lat* at *utc_time*.
     utc_time: datetime.datetime instance of the UTC time
@@ -149,6 +143,7 @@ def cos_zen(utc_time, lon, lat):
     h__ = _local_hour_angle(utc_time, lon, r_a)
     return (np.sin(lat) * np.sin(dec) + np.cos(lat) * np.cos(dec) * np.cos(h__))
 
+
 def sun_zenith_angle(utc_time, lon, lat):
     """Sun-zenith angle for *lon*, *lat* at *utc_time*.
     lon,lat in degrees.
@@ -156,6 +151,7 @@ def sun_zenith_angle(utc_time, lon, lat):
     """
     return np.rad2deg(np.arccos(cos_zen(utc_time, lon, lat)))
 
+
 def sun_earth_distance_correction(utc_time):
     """Calculate the sun earth distance correction, relative to 1 AU.
     """
@@ -174,6 +170,7 @@ def sun_earth_distance_correction(utc_time):
 
     return corr
 
+
 def observer_position(time, lon, lat, alt):
     """Calculate observer ECI position.
 
@@ -197,4 +194,3 @@ def observer_position(time, lon, lat, alt):
     vz = 0
 
     return (x, y, z), (vx, vy, vz)
-    

pyorbital/geoloc.py

 #!/usr/bin/env python
 # -*- coding: utf-8 -*-
 
-# Copyright (c) 2011, 2012, 2013, 2014, 2015.
+# Copyright (c) 2011, 2012, 2013, 2014, 2015, 2017.
 
 # Author(s):
 
@@ -33,7 +33,11 @@ from __future__ import print_function
 
 import numpy as np
 from numpy import cos, sin, sqrt
-from datetime import timedelta
+
+# 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 Orbital
 
 a = 6378.137  # km
@@ -71,7 +75,7 @@ def subpoint(query_point, a=a, b=b):
     ny_ = n__ * cos(lat) * sin(lon)
     nz_ = (1 - e2_) * n__ * sin(lat)
 
-    return np.vstack([nx_, ny_, nz_])
+    return np.stack([nx_, ny_, nz_], axis=0)
 
 
 class ScanGeometry(object):
@@ -89,7 +93,7 @@ class ScanGeometry(object):
                  times,
                  attitude=(0, 0, 0)):
         self.fovs = np.array(fovs)
-        self._times = np.array(times)
+        self._times = np.array(times) * np.timedelta64(1000000000, 'ns')
         self.attitude = attitude
 
     def vectors(self, pos, vel, roll=0.0, pitch=0.0, yaw=0.0):
@@ -117,22 +121,26 @@ class ScanGeometry(object):
         y /= vnorm(y)
 
         # rotate first around x
-        x_rotated = qrotate(nadir, x, self.fovs[:, 0] + roll)
+        x_rotated = qrotate(nadir, x, self.fovs[0] + roll)
         # then around y
-        xy_rotated = qrotate(x_rotated, y, self.fovs[:, 1] + pitch)
+        xy_rotated = qrotate(x_rotated, y, self.fovs[1] + pitch)
         # then around z
         return qrotate(xy_rotated, nadir, yaw)
 
     def times(self, start_of_scan):
-        tds = [timedelta(seconds=i) for i in self._times]
-        return np.array(tds) + start_of_scan
+        #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:
+            return np.array(self._times) + start_of_scan
 
 
 class Quaternion(object):
 
     def __init__(self, scalar, vector):
-        self.__x, self.__y, self.__z = vector
-        self.__w = scalar
+        self.__x, self.__y, self.__z = vector.reshape((3, -1))
+        self.__w = scalar.ravel()
 
     def rotation_matrix(self):
         x, y, z, w = self.__x, self.__y, self.__z, self.__w
@@ -161,22 +169,17 @@ def qrotate(vector, axis, angle):
     """
     n_axis = axis / vnorm(axis)
     sin_angle = np.expand_dims(sin(angle / 2), 0)
-    if np.rank(n_axis) == 1:
+    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__)
+    shape = vector.shape
     return np.einsum("kj, ikj->ij",
-                     vector,
-                     q__.rotation_matrix()[:3, :3])
-
-
-# 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 *
+                     vector.reshape((3, -1)),
+                     q__.rotation_matrix()[:3, :3]).reshape(shape)
 
 
 def get_lonlatalt(pos, utc_time):
@@ -209,10 +212,11 @@ def get_lonlatalt(pos, utc_time):
 # END OF DIRTY STUFF
 
 
-def compute_pixels(tle, sgeom, times, rpy=(0.0, 0.0, 0.0)):
+def compute_pixels(orb, sgeom, times, rpy=(0.0, 0.0, 0.0)):
     """Compute cartesian coordinates of the pixels in instrument scan.
     """
-    (tle1, tle2) = tle
+    if isinstance(orb, (list, tuple)):
+        tle1, tle2 = orb
         orb = Orbital("mysatellite", line1=tle1, line2=tle2)
 
     # get position and velocity for each time of each pixel
@@ -232,8 +236,10 @@ def compute_pixels(tle, sgeom, times, rpy=(0.0, 0.0, 0.0)):
     # b__ = 6356.75231414 # km, GRS80
     b__ = 6356.752314245  # km, WGS84
     radius = np.array([[1 / a__, 1 / a__, 1 / b__]]).T
-    xr_ = vectors * radius
-    cr_ = centre * radius
+    shape = vectors.shape
+
+    xr_ = vectors.reshape([3, -1]) * radius
+    cr_ = centre.reshape([3, -1]) * radius
     ldotc = np.einsum("ij,ij->j", xr_, cr_)
     lsq = np.einsum("ij,ij->j", xr_, xr_)
     csq = np.einsum("ij,ij->j", cr_, cr_)
@@ -241,7 +247,7 @@ def compute_pixels(tle, sgeom, times, rpy=(0.0, 0.0, 0.0)):
     d1_ = (ldotc - np.sqrt(ldotc ** 2 - csq * lsq + lsq)) / lsq
 
     # return the actual pixel positions
-    return vectors * d1_ - centre
+    return vectors * d1_.reshape(shape[1:]) - centre
 
 
 def norm(v):
@@ -252,7 +258,7 @@ def mnorm(m, axis=None):
     """norm of a matrix of vectors stacked along the *axis* dimension.
     """
     if axis is None:
-        axis = np.rank(m) - 1
+        axis = np.ndim(m) - 1
     return np.sqrt((m**2).sum(axis))
 
 

pyorbital/geoloc_instrument_definitions.py

 #!/usr/bin/env python
 # -*- coding: utf-8 -*-
 
-# Copyright (c) 2013, 2014, 2015 Martin Raspaud
+# Copyright (c) 2013, 2014, 2015, 2017 Martin Raspaud
 
 # Author(s):
 
@@ -58,19 +58,21 @@ def avhrr(scans_nb, scan_points,
     # build the avhrr instrument (scan angles)
     avhrr_inst = np.vstack(((scan_points / 1023.5 - 1)
                             * np.deg2rad(-scan_angle),
-                            np.zeros((len(scan_points),)))).transpose()
-    avhrr_inst = np.tile(avhrr_inst, [scans_nb, 1])
+                            np.zeros((len(scan_points),))))
+
+    avhrr_inst = np.tile(
+        avhrr_inst[:, np.newaxis, :], [1, np.int(scans_nb), 1])
 
     # building the corresponding times array
     # times = (np.tile(scan_points * 0.000025 + 0.0025415, [scans_nb, 1])
     #         + np.expand_dims(offset, 1))
 
-    times = np.tile(scan_points * 0.000025, [scans_nb, 1])
+    times = np.tile(scan_points * 0.000025, [np.int(scans_nb), 1])
     if apply_offset:
-        offset = np.arange(scans_nb) * frequency
+        offset = np.arange(np.int(scans_nb)) * frequency
         times += np.expand_dims(offset, 1)
 
-    return ScanGeometry(avhrr_inst, times.ravel())
+    return ScanGeometry(avhrr_inst, times)
 
 
 def avhrr_gac(scan_times, scan_points,
@@ -86,16 +88,17 @@ def avhrr_gac(scan_times, scan_points,
     except TypeError:
         offset = np.arange(scan_times) * frequency
     scans_nb = len(offset)
-    # build the avhrr instrument (scan angles)
+
     avhrr_inst = np.vstack(((scan_points / 1023.5 - 1)
                             * np.deg2rad(-scan_angle),
-                            np.zeros((len(scan_points),)))).transpose()
-    avhrr_inst = np.tile(avhrr_inst, [scans_nb, 1])
+                            np.zeros((len(scan_points),))))
 
+    avhrr_inst = np.tile(
+        avhrr_inst[:, np.newaxis, :], [1, np.int(scans_nb), 1])
     # building the corresponding times array
     times = (np.tile(scan_points * 0.000025, [scans_nb, 1])
              + np.expand_dims(offset, 1))
-    return ScanGeometry(avhrr_inst, times.ravel())
+    return ScanGeometry(avhrr_inst, times)
 
 ################################################################
 # avhrr, all pixels

pyorbital/orbital.py

@@ -30,7 +30,7 @@ from datetime import datetime, timedelta
 
 import numpy as np
 
-from pyorbital import astronomy, tlefile
+from pyorbital import astronomy, dt2np, tlefile
 
 ECC_EPS = 1.0e-6  # Too low for computing further drops.
 ECC_LIMIT_LOW = -1.0e-3
@@ -77,9 +77,9 @@ def get_observer_look(sat_lon, sat_lat, sat_alt, utc_time, lon, lat, alt):
     http://celestrak.com/columns/v02n02/
 
     utc_time: Observation time (datetime object)
-    lon: Longitude of observer position on ground
-    lat: Latitude of observer position on ground
-    alt: Altitude above sea-level (geoid) of observer position on ground
+    lon: Longitude of observer position on ground in degrees east
+    lat: Latitude of observer position on ground in degrees north
+    alt: Altitude above sea-level (geoid) of observer position on ground in km
 
     Return: (Azimuth, Elevation)
     """
@@ -146,7 +146,7 @@ class Orbital(object):
         """
 
         # Propagate backwards to ascending node
-        dt = timedelta(minutes=10)
+        dt = np.timedelta64(10, 'm')
         t_old = utc_time
         t_new = t_old - dt
         pos0, vel0 = self.get_position(t_old, normalize=False)
@@ -226,12 +226,14 @@ class Orbital(object):
         http://celestrak.com/columns/v02n02/
 
         utc_time: Observation time (datetime object)
-        lon: Longitude of observer position on ground
-        lat: Latitude of observer position on ground
-        alt: Altitude above sea-level (geoid) of observer position on ground
+        lon: Longitude of observer position on ground in degrees east
+        lat: Latitude of observer position on ground in degrees north
+        alt: Altitude above sea-level (geoid) of observer position on ground in km
 
         Return: (Azimuth, Elevation)
         """
+
+        utc_time = dt2np(utc_time)
         (pos_x, pos_y, pos_z), (vel_x, vel_y, vel_z) = self.get_position(
             utc_time, normalize=False)
         (opos_x, opos_y, opos_z), (ovel_x, ovel_y, ovel_z) = \
@@ -271,6 +273,7 @@ class Orbital(object):
         """Calculate orbit number at specified time.
         Optionally use TBUS-style orbit numbering (TLE orbit number + 1)
         """
+        utc_time = np.datetime64(utc_time)
         try:
             dt = astronomy._days(utc_time - self.orbit_elements.an_time)
             orbit_period = astronomy._days(self.orbit_elements.an_period)
@@ -287,7 +290,7 @@ class Orbital(object):
 
             self.orbit_elements.an_period = self.orbit_elements.an_time - \
                 self.get_last_an_time(self.orbit_elements.an_time
-                                      - timedelta(minutes=10))
+                                      - np.timedelta64(10, 'm'))
 
             dt = astronomy._days(utc_time - self.orbit_elements.an_time)
             orbit_period = astronomy._days(self.orbit_elements.an_period)
@@ -683,7 +686,12 @@ class _SGDP4(object):
     def propagate(self, utc_time):
         kep = {}
 
-        ts = astronomy._days(utc_time - self.t_0) * XMNPDA
+        # get the time delta in minutes
+        #ts = astronomy._days(utc_time - self.t_0) * XMNPDA
+        # print utc_time.shape
+        # print self.t_0
+        utc_time = dt2np(utc_time)
+        ts = (utc_time - self.t_0) / np.timedelta64(1, 'm')
 
         em = self.eo
         xinc = self.xincl
@@ -796,7 +804,7 @@ class _SGDP4(object):
         xinck = xinc + 1.5 * temp2 * self.cosIO * self.sinIO * cos2u
 
         if np.any(rk < 1):
-            raise Exception('Satellite crased at time %s', utc_time)
+            raise Exception('Satellite crashed at time %s', utc_time)
 
         temp0 = np.sqrt(a)
         temp2 = XKE / (a * temp0)

pyorbital/tests/test_aiaa.py

@@ -24,20 +24,23 @@
 """
 
 # TODO: right formal unit tests.
-from __future__ import with_statement, print_function
+from __future__ import print_function, with_statement
 
 import os
+import unittest
+from datetime import datetime, timedelta
 
-from pyorbital.orbital import Orbital, OrbitElements, _SGDP4
-from pyorbital.tlefile import ChecksumError
-from pyorbital import tlefile, astronomy
 import numpy as np
-from datetime import timedelta, datetime
-import unittest
+
+from pyorbital import astronomy, tlefile
+from pyorbital.orbital import _SGDP4, Orbital, OrbitElements
+from pyorbital.tlefile import ChecksumError
+
 
 class LineOrbital(Orbital):
     """Read TLE lines instead of file.
     """
+
     def __init__(self, satellite, line1, line2):
         satellite = satellite.upper()
         self.satellite_name = satellite
@@ -65,6 +68,7 @@ def get_results(satnumber, delay):
                     utc_time = utc_time.replace(year=int(sline[-4]),
                                                 month=int(sline[-3]),
                                                 day=int(sline[-2]))
+                    utc_time = np.datetime64(utc_time)
                 return (float(sline[1]),
                         float(sline[2]),
                         float(sline[3]),
@@ -74,6 +78,7 @@ def get_results(satnumber, delay):
                         utc_time)
             line = f_2.readline()
 
+
 class AIAAIntegrationTest(unittest.TestCase):
     """Test against the AIAA test cases.
     """
@@ -106,12 +111,15 @@ class AIAAIntegrationTest(unittest.TestCase):
                         test_line = f__.readline()
                         continue
                     except ChecksumError as e:
-                        self.assertTrue(test_line.split()[1] in ["33333", "33334", "33335"])
+                        self.assertTrue(test_line.split()[1] in [
+                                        "33333", "33334", "33335"])
                     for delay in times:
                         try:
-                            test_time = timedelta(minutes=delay) + o.tle.epoch
+                            test_time = delay.astype(
+                                'timedelta64[m]') + o.tle.epoch
                             pos, vel = o.get_position(test_time, False)
-                            res = get_results(int(o.tle.satnumber), float(delay))
+                            res = get_results(
+                                int(o.tle.satnumber), float(delay))
                         except NotImplementedError:
                             # Skipping deep-space
                             break
@@ -144,4 +152,3 @@ def suite():
     mysuite.addTest(loader.loadTestsFromTestCase(AIAAIntegrationTest))
 
     return mysuite
-

pyorbital/tests/test_geoloc.py

 #!/usr/bin/env python
 # -*- coding: utf-8 -*-
 
-# Copyright (c) 2014 Martin Raspaud
+# Copyright (c) 2014, 2017 Martin Raspaud
 
 # Author(s):
 
@@ -24,13 +24,16 @@
 """
 
 import unittest
-from pyorbital.geoloc_instrument_definitions import avhrr
-from pyorbital.geoloc import (ScanGeometry, qrotate,
-                              subpoint, geodetic_lat)
-import numpy as np
 from datetime import datetime, timedelta
 
+import numpy as np
+
+from pyorbital.geoloc import ScanGeometry, geodetic_lat, qrotate, subpoint
+from pyorbital.geoloc_instrument_definitions import avhrr
+
+
 class TestQuaternion(unittest.TestCase):
+
     """Test the quaternion rotation.
     """
 
@@ -62,45 +65,59 @@ class TestQuaternion(unittest.TestCase):
 
 
 class TestGeoloc(unittest.TestCase):
+
     """Test for the core computing part.
     """
 
     def test_scan_geometry(self):
         """Test the ScanGeometry object.
         """
-        instrument = ScanGeometry(np.deg2rad(np.array([[10, 0],
-                                                       [0, 0],
-                                                       [-10, 0]])),
-                                  np.array([-0.1, 0, 0.1]))
-        self.assertTrue(np.allclose(np.rad2deg(instrument.fovs[:, 0]),
-                                    np.array([10, 0, -10])))
+        scans_nb = 1
+
+        xy = np.vstack((np.deg2rad(np.array([10, 0, -10])),
+                        np.array([0, 0, 0])))
+        xy = np.tile(xy[:, np.newaxis, :], [1, np.int(scans_nb), 1])
+
+        times = np.tile([-0.1, 0, 0.1], [np.int(scans_nb), 1])
+
+        instrument = ScanGeometry(xy, times)
+
+        self.assertTrue(np.allclose(np.rad2deg(instrument.fovs[0]),
+                                    np.array([[10, 0, -10]])))
 
         # Test vectors
 
-        pos = np.array([[0, 0, 7000]]).T
-        vel = np.array([[1, 0, 0]]).T
+        pos = np.rollaxis(np.tile(np.array([0, 0, 7000]), [3, 1, 1]), 2)
+        vel = np.rollaxis(np.tile(np.array([1, 0, 0]), [3, 1, 1]), 2)
+        pos = np.stack([np.array([0, 0, 7000])] * 3, 1)[:, np.newaxis, :]
+        vel = np.stack([np.array([1, 0, 0])] * 3, 1)[:, np.newaxis, :]
+
         vec = instrument.vectors(pos, vel)
+
         self.assertTrue(np.allclose(np.array([[0, 0, -1]]),
-                                    vec[:, 1]))
+                                    vec[:, 0, 1]))
 
         # minus sin because we use trigonometrical direction of angles
 
         self.assertTrue(np.allclose(np.array([[0,
                                                -np.sin(np.deg2rad(10)),
                                                -np.cos(np.deg2rad(10))]]),
-                                    vec[:, 0]))
+                                    vec[:, 0, 0]))
         self.assertTrue(np.allclose(np.array([[0,
                                                -np.sin(np.deg2rad(-10)),
                                                -np.cos(np.deg2rad(-10))]]),
-                                    vec[:, 2]))
+                                    vec[:, 0, 2]))
 
         # Test times
 
-        start_of_scan = datetime(2014, 1, 8, 11, 30)
+        start_of_scan = np.datetime64(datetime(2014, 1, 8, 11, 30))
         times = instrument.times(start_of_scan)
-        self.assertEqual(times[1], start_of_scan)
-        self.assertEqual(times[0], start_of_scan - timedelta(seconds=0.1))
-        self.assertEqual(times[2], start_of_scan + timedelta(seconds=0.1))
+
+        self.assertEquals(times[0, 1], start_of_scan)
+        self.assertEquals(times[0, 0], start_of_scan -
+                          np.timedelta64(100, 'ms'))
+        self.assertEquals(times[0, 2], start_of_scan +
+                          np.timedelta64(100, 'ms'))
 
     def test_geodetic_lat(self):
         """Test the determination of the geodetic latitude.
@@ -117,7 +134,6 @@ class TestGeoloc(unittest.TestCase):
                                     np.array([0.78755832699854733,
                                               0.78755832699854733])))
 
-        
     def test_subpoint(self):
         """Test nadir determination.
         """
@@ -125,39 +141,50 @@ class TestGeoloc(unittest.TestCase):
         b = 6356.75231414  # km, GRS80
         point = np.array([0, 0, 7000])
         nadir = subpoint(point, a, b)
-        self.assertTrue(np.allclose(nadir, np.array([[0, 0, b]]).T))
+        self.assertTrue(np.allclose(nadir, np.array([[0, 0, b]])))
 
         point = np.array([7000, 0, 7000])
         nadir = subpoint(point, a, b)
         self.assertTrue(np.allclose(nadir,
                                     np.array([[4507.85431429,
                                                0,
-                                               4497.06396339]]).T))
+                                               4497.06396339]])))
         points = np.array([[7000, 0, 7000],
                            [7000, 0, 7000]]).T
         nadir = subpoint(points, a, b)
-        self.assertTrue(np.allclose(nadir,
+        self.assertTrue(np.allclose(nadir[:, 0],
                                     np.array([[4507.85431429,
                                                0,
-                                               4497.06396339]]).T))
-
-
+                                               4497.06396339]])))
+        self.assertTrue(np.allclose(nadir[:, 1],
+                                    np.array([[4507.85431429,
+                                               0,
+                                               4497.06396339]])))
 
 
 class TestGeolocDefs(unittest.TestCase):
+
     """Test the instrument definitions.
     """
+
     def test_avhrr(self):
         """Test the definition of the avhrr instrument
         """
         avh = avhrr(1, np.array([0, 1023.5, 2047]))
-        self.assertTrue(np.allclose(np.rad2deg(avh.fovs[:, 0]),
+        self.assertTrue(np.allclose(np.rad2deg(avh.fovs[0]),
                                     np.array([55.37, 0, -55.37])))
 
         avh = avhrr(1, np.array([0, 1023.5, 2047]), 10)
-        self.assertTrue(np.allclose(np.rad2deg(avh.fovs[:, 0]),
+        self.assertTrue(np.allclose(np.rad2deg(avh.fovs[0]),
                                     np.array([10, 0, -10])))
 
+        # This is perhaps a bit odd, to require avhrr to accept floats for
+        # the number of scans? FIXME!
+        avh = avhrr(1.1, np.array([0, 1023.5, 2047]), 10)
+        self.assertTrue(np.allclose(np.rad2deg(avh.fovs[0]),
+                                    np.array([10, 0, -10])))
+
+
 def suite():
     """The suite for test_geoloc
     """

pyorbital/tests/test_orbital.py

@@ -32,6 +32,7 @@ from pyorbital import orbital
 
 eps_deg = 10e-3
 
+
 class Test(unittest.TestCase):
 
     def test_get_orbit_number(self):
@@ -56,7 +57,6 @@ class Test(unittest.TestCase):
         self.assertTrue(np.abs(lat - expected_lat) < eps_deg, 'Calculation of sublat failed')
         self.assertTrue(np.abs(alt - expected_alt) < eps_deg, 'Calculation of altitude failed')
         
-
     def test_observer_look(self):
         sat = orbital.Orbital("ISS (ZARYA)",
                               line1="1 25544U 98067A   03097.78853147  .00021906  00000-0  28403-3 0  8652",
@@ -79,7 +79,7 @@ class Test(unittest.TestCase):
         sat = orbital.Orbital("METOP-A",
                               line1="1 29499U 06044A   13060.48822809  .00000017  00000-0  27793-4 0  9819",
                               line2="2 29499  98.6639 121.6164 0001449  71.9056  43.3132 14.21510544330271")
-        dt = timedelta(minutes=98)
+        dt = np.timedelta64(98, 'm')
         self.assertEqual(sat.get_orbit_number(sat.tle.epoch + dt), 33028)
 
     def test_orbit_num_equator(self):
@@ -99,7 +99,6 @@ class Test(unittest.TestCase):
         self.assertTrue(pos2[2] > 0)
 
 
-
 def suite():
     """The suite for test_orbital
     """

pyorbital/tlefile.py

@@ -32,9 +32,14 @@ except ImportError:
     from urllib.request import urlopen
 import os
 import glob
+import numpy as np
 
 TLE_URLS = ('http://celestrak.com/NORAD/elements/weather.txt',
-            'http://celestrak.com/NORAD/elements/resource.txt')
+            'http://celestrak.com/NORAD/elements/resource.txt',
+            'https://www.celestrak.com/NORAD/elements/cubesat.txt',
+            'http://celestrak.com/NORAD/elements/stations.txt',
+            'https://www.celestrak.com/NORAD/elements/sarsat.txt',
+            'https://www.celestrak.com/NORAD/elements/noaa.txt')
 
 LOGGER = logging.getLogger(__name__)
 
@@ -56,6 +61,8 @@ def read_platform_numbers(in_upper=False, num_as_int=False):
             # skip comment lines
             if not row.startswith('#'):
                 parts = row.split()
+                if len(parts) < 2:
+                    continue
                 if in_upper:
                     parts[0] = parts[0].upper()
                 if num_as_int:
@@ -142,19 +149,21 @@ def read(platform, tle_file=None, line1=None, line2=None):
 def fetch(destination):
     """fetch TLE from internet and save it to *destination*.
    """
-    with open(destination, "w") as dest:
+    with io.open(destination, mode="w", encoding="utf-8") as dest:
         for url in TLE_URLS:
             response = urlopen(url)
-            dest.write(response.read())
+            dest.write(response.read().decode("utf-8"))
 
 
 class ChecksumError(Exception):
+
     '''ChecksumError.
     '''
     pass
 
 
 class Tle(object):
+
     """Class holding TLE objects.
    """
 
@@ -294,8 +303,9 @@ class Tle(object):
         self.id_launch_piece = self._line1[14:17]
         self.epoch_year = self._line1[18:20]
         self.epoch_day = float(self._line1[20:32])
-        self.epoch = (datetime.datetime.strptime(self.epoch_year, "%y") +
-                      datetime.timedelta(days=self.epoch_day - 1))
+        self.epoch = \
+            np.datetime64(datetime.datetime.strptime(self.epoch_year, "%y") +
+                          datetime.timedelta(days=self.epoch_day - 1), 'us')
         self.mean_motion_derivative = float(self._line1[33:43])
         self.mean_motion_sec_derivative = _read_tle_decimal(self._line1[44:52])
         self.bstar = _read_tle_decimal(self._line1[53:61])
@@ -326,6 +336,7 @@ class Tle(object):
         pprint.pprint(d_var, s_var)
         return s_var.getvalue()[:-1]
 
+
 def main():
     '''Main for testing TLE reading.
     '''

pyorbital/version.py

@@ -23,4 +23,4 @@
 """Version file.
 """
 
-__version__ = "v1.1.1"
+__version__ = "v1.2.0"