This is an automated email from the git hooks/post-receive script. a_valentino-guest pushed a commit to branch master in repository pyorbital.
commit 2fd39a2114e81e089b3ef486471b2af6b6bad0d4 Author: Antonio Valentino <antonio.valent...@tiscali.it> Date: Sun Feb 18 09:47:50 2018 +0100 New upstream version 1.2.0 --- .bumpversion.cfg | 2 +- .github/ISSUE_TEMPLATE.md | 18 +++++ .github/PULL_REQUEST_TEMPLATE.md | 9 +++ changelog.rst | 123 +++++++++++++++++++++++++++++ doc/source/index.rst | 36 ++++++++- etc/platforms.txt | 5 ++ pyorbital/__init__.py | 29 +++++++ pyorbital/astronomy.py | 72 ++++++++--------- pyorbital/geoloc.py | 58 ++++++++------ pyorbital/geoloc_instrument_definitions.py | 23 +++--- pyorbital/orbital.py | 30 ++++--- pyorbital/tests/test_aiaa.py | 45 ++++++----- pyorbital/tests/test_geoloc.py | 107 +++++++++++++++---------- pyorbital/tests/test_orbital.py | 45 ++++++----- pyorbital/tlefile.py | 21 +++-- pyorbital/version.py | 2 +- 16 files changed, 447 insertions(+), 178 deletions(-) diff --git a/.bumpversion.cfg b/.bumpversion.cfg index 01c9521..234e98e 100644 --- a/.bumpversion.cfg +++ b/.bumpversion.cfg @@ -1,5 +1,5 @@ [bumpversion] -current_version = 1.1.1 +current_version = 1.2.0 commit = True tag = True diff --git a/.github/ISSUE_TEMPLATE.md b/.github/ISSUE_TEMPLATE.md new file mode 100644 index 0000000..18c258a --- /dev/null +++ b/.github/ISSUE_TEMPLATE.md @@ -0,0 +1,18 @@ +#### 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 ! diff --git a/.github/PULL_REQUEST_TEMPLATE.md b/.github/PULL_REQUEST_TEMPLATE.md new file mode 100644 index 0000000..346c742 --- /dev/null +++ b/.github/PULL_REQUEST_TEMPLATE.md @@ -0,0 +1,9 @@ +<!-- 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 --> diff --git a/changelog.rst b/changelog.rst index be5c2a0..c2f96a7 100644 --- a/changelog.rst +++ b/changelog.rst @@ -1,6 +1,129 @@ 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.dybb...@smhi.se> + + +- Add github issue and PR templates. [Adam.Dybbroe] + + Signed-off-by: Adam.Dybbroe <adam.dybb...@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.dybb...@smhi.se> + + +- Add GOES-16 and Himawari-9. [Adam.Dybbroe] + + Signed-off-by: Adam.Dybbroe <adam.dybb...@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 <a000...@c20671.ad.smhi.se> + + +- Numpy 1.12 compatible - accepting number of scans to be a float. + [Adam.Dybbroe] + + Signed-off-by: Adam.Dybbroe <a000...@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.rasp...@smhi.se> + + +- Convert input times to datetime64. [Martin Raspaud] + + Signed-off-by: Martin Raspaud <martin.rasp...@smhi.se> + + +- Allow offset application to be turned off (avhrr) [Martin Raspaud] + + Signed-off-by: Martin Raspaud <martin.rasp...@smhi.se> + + +- Adapt to datetime64. [Martin Raspaud] + + Signed-off-by: Martin Raspaud <martin.rasp...@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) ------------------- diff --git a/doc/source/index.rst b/doc/source/index.rst index a330d7d..a14b425 100644 --- a/doc/source/index.rst +++ b/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 --------------------------------- diff --git a/etc/platforms.txt b/etc/platforms.txt index 6fc5250..14a716c 100644 --- a/etc/platforms.txt +++ b/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 diff --git a/pyorbital/__init__.py b/pyorbital/__init__.py index e69de29..2ce9000 100644 --- a/pyorbital/__init__.py +++ b/pyorbital/__init__.py @@ -0,0 +1,29 @@ +# -*- coding: utf-8 -*- + +# Copyright (c) 2017 + +# Author(s): + +# Martin Raspaud <martin.rasp...@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]') diff --git a/pyorbital/astronomy.py b/pyorbital/astronomy.py index 3f76efb..cd49f89 100644 --- a/pyorbital/astronomy.py +++ b/pyorbital/astronomy.py @@ -24,46 +24,36 @@ 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 -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): """Get the julian day of *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): """Greenwich mean sidereal utc_time, in radians. - + As defined in the AIAA 2006 implementation: http://www.celestrak.com/publications/AIAA/2006-6753/ """ @@ -86,24 +76,25 @@ def sun_ecliptic_longitude(utc_time): jdate = jdays2000(utc_time) / 36525.0 # mean anomaly, rad m_a = np.deg2rad(357.52910 + - 35999.05030*jdate - - 0.0001559*jdate*jdate - - 0.00000048*jdate*jdate*jdate) + 35999.05030 * jdate - + 0.0001559 * jdate * jdate - + 0.00000048 * jdate * jdate * jdate) # mean longitude, deg - l_0 = 280.46645 + 36000.76983*jdate + 0.0003032*jdate*jdate - d_l = ((1.914600 - 0.004817*jdate - 0.000014*jdate*jdate)*np.sin(m_a) + - (0.019993 - 0.000101*jdate)*np.sin(2*m_a) + 0.000290*np.sin(3*m_a)) + l_0 = 280.46645 + 36000.76983 * jdate + 0.0003032 * jdate * jdate + d_l = ((1.914600 - 0.004817 * jdate - 0.000014 * jdate * jdate) * np.sin(m_a) + + (0.019993 - 0.000101 * jdate) * np.sin(2 * m_a) + 0.000290 * np.sin(3 * m_a)) # true longitude, deg 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*. """ jdate = jdays2000(utc_time) / 36525.0 - eps = np.deg2rad(23.0 + 26.0/60.0 + 21.448/3600.0 - - (46.8150*jdate + 0.00059*jdate*jdate - - 0.001813*jdate*jdate*jdate) / 3600) + eps = np.deg2rad(23.0 + 26.0 / 60.0 + 21.448 / 3600.0 - + (46.8150 * jdate + 0.00059 * jdate * jdate - + 0.001813 * jdate * jdate * jdate) / 3600) eclon = sun_ecliptic_longitude(utc_time) x__ = np.cos(eclon) y__ = np.cos(eps) * np.sin(eclon) @@ -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 @@ -132,10 +125,11 @@ def get_alt_az(utc_time, lon, lat): ra_, dec = sun_ra_dec(utc_time) h__ = _local_hour_angle(utc_time, lon, ra_) - return (np.arcsin(np.sin(lat)*np.sin(dec) + + return (np.arcsin(np.sin(lat) * np.sin(dec) + np.cos(lat) * np.cos(dec) * np.cos(h__)), - np.arctan2(-np.sin(h__), (np.cos(lat)*np.tan(dec) - - np.sin(lat)*np.cos(h__)))) + 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*. @@ -147,7 +141,8 @@ def cos_zen(utc_time, lon, lat): r_a, dec = sun_ra_dec(utc_time) 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__)) + 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*. @@ -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. """ @@ -170,19 +166,20 @@ def sun_earth_distance_correction(utc_time): # corr_me = (r / AU) ** 2 # from known software. - corr = 1 - 0.0334 * np.cos(2 * np.pi * (jdays2000(utc_time) - 2) / year) + corr = 1 - 0.0334 * np.cos(2 * np.pi * (jdays2000(utc_time) - 2) / year) return corr + def observer_position(time, lon, lat, alt): """Calculate observer ECI position. http://celestrak.com/columns/v02n03/ """ - + lon = np.deg2rad(lon) lat = np.deg2rad(lat) - + theta = (gmst(time) + lon) % (2 * np.pi) c = 1 / np.sqrt(1 + F * (F - 2) * np.sin(lat)**2) sq = c * (1 - F)**2 @@ -192,9 +189,8 @@ def observer_position(time, lon, lat, alt): y = achcp * np.sin(theta) z = (A * sq + alt) * np.sin(lat) - vx = -MFACTOR*y # kilometers/second - vy = MFACTOR*x + vx = -MFACTOR * y # kilometers/second + vy = MFACTOR * x vz = 0 return (x, y, z), (vx, vy, vz) - diff --git a/pyorbital/geoloc.py b/pyorbital/geoloc.py index 2589705..46b7519 100644 --- a/pyorbital/geoloc.py +++ b/pyorbital/geoloc.py @@ -1,7 +1,7 @@ #!/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,11 +212,12 @@ 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 - orb = Orbital("mysatellite", line1=tle1, line2=tle2) + 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 pos, vel = orb.get_position(times, normalize=False) @@ -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)) diff --git a/pyorbital/geoloc_instrument_definitions.py b/pyorbital/geoloc_instrument_definitions.py index b8a3c2c..87967b0 100644 --- a/pyorbital/geoloc_instrument_definitions.py +++ b/pyorbital/geoloc_instrument_definitions.py @@ -1,7 +1,7 @@ #!/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 diff --git a/pyorbital/orbital.py b/pyorbital/orbital.py index 431578a..d7e99a6 100644 --- a/pyorbital/orbital.py +++ b/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) diff --git a/pyorbital/tests/test_aiaa.py b/pyorbital/tests/test_aiaa.py index 5c1fd07..a820e68 100644 --- a/pyorbital/tests/test_aiaa.py +++ b/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 @@ -55,7 +58,7 @@ def get_results(satnumber, delay): while(line): if line.endswith(" xx\n") and int(line[:-3]) == satnumber: line = f_2.readline() - while(not line.startswith("%.8f"%delay)): + while(not line.startswith("%.8f" % delay)): line = f_2.readline() sline = line.split() if delay == 0: @@ -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,10 +78,11 @@ def get_results(satnumber, delay): utc_time) line = f_2.readline() + class AIAAIntegrationTest(unittest.TestCase): """Test against the AIAA test cases. """ - + def test_aiaa(self): """Do the tests against 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 @@ -119,10 +127,10 @@ class AIAAIntegrationTest(unittest.TestCase): # from warnings import warn # warn(test_name + ' ' + str(e)) # break - - delta_pos = 5e-6 # km = 5 mm - delta_vel = 5e-9 # km/s = 5 um/s - delta_time = 1e-3 # 1 milisecond + + delta_pos = 5e-6 # km = 5 mm + delta_vel = 5e-9 # km/s = 5 um/s + delta_time = 1e-3 # 1 milisecond self.assertTrue(abs(res[0] - pos[0]) < delta_pos) self.assertTrue(abs(res[1] - pos[1]) < delta_pos) self.assertTrue(abs(res[2] - pos[2]) < delta_pos) @@ -132,9 +140,9 @@ class AIAAIntegrationTest(unittest.TestCase): if res[6] is not None: dt = astronomy._days(res[6] - test_time) * 24 * 60 self.assertTrue(abs(dt) < delta_time) - + test_line = f__.readline() - + def suite(): """The suite for test_aiaa @@ -142,6 +150,5 @@ def suite(): loader = unittest.TestLoader() mysuite = unittest.TestSuite() mysuite.addTest(loader.loadTestsFromTestCase(AIAAIntegrationTest)) - - return mysuite + return mysuite diff --git a/pyorbital/tests/test_geoloc.py b/pyorbital/tests/test_geoloc.py index 79c1fbb..c105e80 100644 --- a/pyorbital/tests/test_geoloc.py +++ b/pyorbital/tests/test_geoloc.py @@ -1,7 +1,7 @@ #!/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. """ @@ -59,55 +62,69 @@ class TestQuaternion(unittest.TestCase): self.assertTrue(np.allclose(qrotate(vector, axis, angle), np.array([[0, 0, 1], [-1, 0, 0]]).T)) - + 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. """ - a = 6378.137 # km - b = 6356.75231414 # km, GRS80 + a = 6378.137 # km + b = 6356.75231414 # km, GRS80 point = np.array([7000, 0, 7000]) self.assertEqual(geodetic_lat(point), 0.78755832699854733) @@ -116,48 +133,58 @@ class TestGeoloc(unittest.TestCase): self.assertTrue(np.allclose(geodetic_lat(points), np.array([0.78755832699854733, 0.78755832699854733]))) - - + def test_subpoint(self): """Test nadir determination. """ - a = 6378.137 # km - b = 6356.75231414 # km, GRS80 + a = 6378.137 # km + 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]), + avh = avhrr(1, np.array([0, 1023.5, 2047])) + 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]), + avh = avhrr(1, np.array([0, 1023.5, 2047]), 10) + 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 """ @@ -166,5 +193,5 @@ def suite(): mysuite.addTest(loader.loadTestsFromTestCase(TestQuaternion)) mysuite.addTest(loader.loadTestsFromTestCase(TestGeoloc)) mysuite.addTest(loader.loadTestsFromTestCase(TestGeolocDefs)) - + return mysuite diff --git a/pyorbital/tests/test_orbital.py b/pyorbital/tests/test_orbital.py index 5238136..f67735d 100644 --- a/pyorbital/tests/test_orbital.py +++ b/pyorbital/tests/test_orbital.py @@ -32,21 +32,22 @@ from pyorbital import orbital eps_deg = 10e-3 + class Test(unittest.TestCase): def test_get_orbit_number(self): """Testing getting the orbitnumber from the tle""" sat = orbital.Orbital("NPP", - line1="1 37849U 11061A 12017.90990040 -.00000112 00000-0 -32693-4 0 772", - line2="2 37849 98.7026 317.8811 0001845 92.4533 267.6830 14.19582686 11574") + line1="1 37849U 11061A 12017.90990040 -.00000112 00000-0 -32693-4 0 772", + line2="2 37849 98.7026 317.8811 0001845 92.4533 267.6830 14.19582686 11574") dobj = datetime(2012, 1, 18, 8, 4, 19) orbnum = sat.get_orbit_number(dobj) self.assertEqual(orbnum, 1163) def test_sublonlat(self): - sat = orbital.Orbital("ISS (ZARYA)", - line1="1 25544U 98067A 03097.78853147 .00021906 00000-0 28403-3 0 8652", - line2="2 25544 51.6361 13.7980 0004256 35.6671 59.2566 15.58778559250029") + sat = orbital.Orbital("ISS (ZARYA)", + line1="1 25544U 98067A 03097.78853147 .00021906 00000-0 28403-3 0 8652", + line2="2 25544 51.6361 13.7980 0004256 35.6671 59.2566 15.58778559250029") d = datetime(2003, 3, 23, 0, 3, 22) lon, lat, alt = sat.get_lonlatalt(d) expected_lon = -68.199894472013213 @@ -56,11 +57,10 @@ 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", - line2="2 25544 51.6361 13.7980 0004256 35.6671 59.2566 15.58778559250029") + sat = orbital.Orbital("ISS (ZARYA)", + line1="1 25544U 98067A 03097.78853147 .00021906 00000-0 28403-3 0 8652", + line2="2 25544 51.6361 13.7980 0004256 35.6671 59.2566 15.58778559250029") d = datetime(2003, 3, 23, 0, 3, 22) az, el = sat.get_observer_look(d, -84.39733, 33.775867, 0) expected_az = 122.45169655331965 @@ -69,23 +69,23 @@ class Test(unittest.TestCase): self.assertTrue(np.abs(el - expected_el) < eps_deg, 'Calculation of elevation failed') def test_orbit_num_an(self): - sat = orbital.Orbital("METOP-A", - line1="1 29499U 06044A 11254.96536486 .00000092 00000-0 62081-4 0 5221", - line2="2 29499 98.6804 312.6735 0001758 111.9178 248.2152 14.21501774254058") + sat = orbital.Orbital("METOP-A", + line1="1 29499U 06044A 11254.96536486 .00000092 00000-0 62081-4 0 5221", + line2="2 29499 98.6804 312.6735 0001758 111.9178 248.2152 14.21501774254058") d = datetime(2011, 9, 14, 5, 30) self.assertEqual(sat.get_orbit_number(d), 25437) - + def test_orbit_num_non_an(self): - 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) + 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 = np.timedelta64(98, 'm') self.assertEqual(sat.get_orbit_number(sat.tle.epoch + dt), 33028) - + def test_orbit_num_equator(self): - sat = orbital.Orbital("SUOMI NPP", - line1="1 37849U 11061A 13061.24611272 .00000048 00000-0 43679-4 0 4334", - line2="2 37849 98.7444 1.0588 0001264 63.8791 102.8546 14.19528338 69643") + sat = orbital.Orbital("SUOMI NPP", + line1="1 37849U 11061A 13061.24611272 .00000048 00000-0 43679-4 0 4334", + line2="2 37849 98.7444 1.0588 0001264 63.8791 102.8546 14.19528338 69643") t1 = datetime(2013, 3, 2, 22, 2, 25) t2 = datetime(2013, 3, 2, 22, 2, 26) on1 = sat.get_orbit_number(t1) @@ -99,12 +99,11 @@ class Test(unittest.TestCase): self.assertTrue(pos2[2] > 0) - def suite(): """The suite for test_orbital """ loader = unittest.TestLoader() mysuite = unittest.TestSuite() mysuite.addTest(loader.loadTestsFromTestCase(Test)) - + return mysuite diff --git a/pyorbital/tlefile.py b/pyorbital/tlefile.py index d2d016d..114f8d8 100644 --- a/pyorbital/tlefile.py +++ b/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. ''' diff --git a/pyorbital/version.py b/pyorbital/version.py index 46e3a50..ac8c7c2 100644 --- a/pyorbital/version.py +++ b/pyorbital/version.py @@ -23,4 +23,4 @@ """Version file. """ -__version__ = "v1.1.1" +__version__ = "v1.2.0" -- Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/pkg-grass/pyorbital.git _______________________________________________ Pkg-grass-devel mailing list Pkg-grass-devel@lists.alioth.debian.org http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/pkg-grass-devel