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 d01aa734c12d66c9c78e6b1b1d3d96c608d77a70 Author: Antonio Valentino <[email protected]> Date: Mon Feb 15 22:48:42 2016 +0100 Imported Upstream version 1.0.0 --- changelog.rst | 489 +++++++++++++++++++++++++++++ doc/source/conf.py | 15 +- etc/platforms.txt | 60 ++++ pyorbital/geoloc.py | 111 +++---- pyorbital/geoloc_instrument_definitions.py | 102 +++--- pyorbital/orbital.py | 186 ++++++----- pyorbital/tlefile.py | 267 ++++++++++++---- pyorbital/version.py | 8 +- setup.cfg | 5 + setup.py | 2 +- 10 files changed, 1009 insertions(+), 236 deletions(-) diff --git a/changelog.rst b/changelog.rst new file mode 100644 index 0000000..0783152 --- /dev/null +++ b/changelog.rst @@ -0,0 +1,489 @@ +Changelog +========= + +v1.0.0 (2015-08-25) +------------------- + +- Update changelog. [Martin Raspaud] + +- Bump version: 0.3.2 → 1.0.0. [Martin Raspaud] + +- Cleanup. [Martin Raspaud] + + Signed-off-by: Martin Raspaud <[email protected]> + + +- Fix version number. [Martin Raspaud] + + Signed-off-by: Martin Raspaud <[email protected]> + + +- Cosmetics. [Martin Raspaud] + + Signed-off-by: Martin Raspaud <[email protected]> + + +- Merge pull request #2 from pnuu/feature_tle_lookup. [Martin Raspaud] + + Use NORAD catalog numbers for TLE reading + +- Example file for mapping OSCAR platform names and NORAD catalog + numbers. [Panu Lahtinen] + +- Add setup.cfg for easy rpm generation. [Martin Raspaud] + + Signed-off-by: Martin Raspaud <[email protected]> + + +- Merge branch 'develop' of github.com:mraspaud/pyorbital into develop. + [Martin Raspaud] + +- Merge pull request #1 from spareeth/develop. [Martin Raspaud] + + changes to avhrr_gacfunction and read_tle_decimal + +- Added '+' as a condition in the read_tle function. [Sajid Pareeth] + +- Renaming the variable scans_nb to scan_times in offset in avhrr_gac + function. [Sajid Pareeth] + +- Bugfix: eccentricity too low message formatting. [Martin Raspaud] + + Signed-off-by: Martin Raspaud <[email protected]> + + +- Allow reading TLE from the most recent file described by the TLES env. + [Martin Raspaud] + + Signed-off-by: Martin Raspaud <[email protected]> + + +- Change decimate to frequency in avhrr instruments. [Martin Raspaud] + + Signed-off-by: Martin Raspaud <[email protected]> + + +- Add the avhrr instrument, gac version. [Martin Raspaud] + + Signed-off-by: Martin Raspaud <[email protected]> + + +- Accept missing zeros in TLE (old noaa compatibility). [Martin Raspaud] + +- Add the horizon parameter to get_next_passes to get the + risetime/falltime at given angle. [Martin Raspaud] + + Signed-off-by: Martin Raspaud <[email protected]> + + +- Merge branch 'master' into develop. [Martin Raspaud] + +- Fix backwards numpy compatibility. [Martin Raspaud] + + Signed-off-by: Martin Raspaud <[email protected]> + + +v0.3.2 (2014-04-10) +------------------- + +- Merge branch 'develop' [Martin Raspaud] + +- Bump up version number. [Martin Raspaud] + + Signed-off-by: Martin Raspaud <[email protected]> + + +- Merge branch 'feature-no-scipy' into develop. [Martin Raspaud] + +- Remove scipy dependencies. [Martin Raspaud] + + Was depending on scipy.optimize, brent and brentq function. + Replaced by secant method root finding and successive parabolic + interpolation local minimum finding. + + Signed-off-by: Martin Raspaud <[email protected]> + + +- Correcting the travis file. [Martin Raspaud] + + Signed-off-by: Martin Raspaud <[email protected]> + + +v0.3.1 (2014-02-24) +------------------- + +- Bugfix in travis file. [Martin Raspaud] + +- Bump up version number. [Martin Raspaud] + + Signed-off-by: Martin Raspaud <[email protected]> + + +- Fixed documentation. [Martin Raspaud] + + Signed-off-by: Martin Raspaud <[email protected]> + + +- Cleanup. [Martin Raspaud] + + Signed-off-by: Martin Raspaud <[email protected]> + + +- New nadir computations for geoloc. [Martin Raspaud] + + Signed-off-by: Martin Raspaud <[email protected]> + + +- More unit tests. [Martin Raspaud] + + Signed-off-by: Martin Raspaud <[email protected]> + + +v0.3.0 (2014-01-07) +------------------- + +- Auto update version number in documentation. [Martin Raspaud] + + Signed-off-by: Martin Raspaud <[email protected]> + + +- Change to version file and bump up to v0.3.0. [Martin Raspaud] + + Signed-off-by: Martin Raspaud <[email protected]> + + +- Cleanup the testfiles. [Martin Raspaud] + + Signed-off-by: Martin Raspaud <[email protected]> + + +- Add a test to read tle from file. [Martin Raspaud] + + Signed-off-by: Martin Raspaud <[email protected]> + + +- Fix doc path in MANIFEST.in. [Martin Raspaud] + + Signed-off-by: Martin Raspaud <[email protected]> + + +v0.2.4 (2014-01-07) +------------------- + +- Merge branch 'feature-travis' into pre-master. [Martin Raspaud] + +- Add test for tle reading, cleanup and make ready for travis. [Martin + Raspaud] + + Signed-off-by: Martin Raspaud <[email protected]> + + +- Cleanup. [Martin Raspaud] + + Signed-off-by: Martin Raspaud <[email protected]> + + +- Add function to fetch the tle files from internet manually. [Martin + Raspaud] + + Signed-off-by: Martin Raspaud <[email protected]> + + +- Adding the viirs instrument. [Martin Raspaud] + + Signed-off-by: Martin Raspaud <[email protected]> + + +- Change sphinx theme. [Martin Raspaud] + + Signed-off-by: Martin Raspaud <[email protected]> + + +- Fix doc for readthedocs. [Martin Raspaud] + + Signed-off-by: Martin Raspaud <[email protected]> + + +- Remove unused old file. [Martin Raspaud] + + Signed-off-by: Martin Raspaud <[email protected]> + + +- Merge branch 'geoloc' into pre-master. [Martin Raspaud] + +- Work on geolocation. [Martin Raspaud] + + Signed-off-by: Martin Raspaud <[email protected]> + + +- Numpyze the orbital computation. [Martin Raspaud] + + Signed-off-by: Martin Raspaud <[email protected]> + + +- Add some logging in tle file fetching. [Martin Raspaud] + + Signed-off-by: Martin Raspaud <[email protected]> + + +- Fix syntax error in doc/conf.py. [Martin Raspaud] + +- Make the scan angle of avhrr an argument. [Martin Raspaud] + +- Factorize avhrr code (geoloc definition) [Martin Raspaud] + + Signed-off-by: Martin Raspaud <[email protected]> + + +- Add Mikhail's definition of AMSU-A. [Martin Raspaud] + +- Add instrument examples for geoloc. [Martin Raspaud] + + Signed-off-by: Martin Raspaud <[email protected]> + + +- Merge branch 'geoloc' of github.com:mraspaud/pyorbital into geoloc. + [Martin Raspaud] + +- Try fixing nadir. [Martin Raspaud] + +- Fix attitude. [Martin Raspaud] + +- Updated doc and copyright. [Martin Raspaud] + +- Add geoloc example. [Martin Raspaud] + + Signed-off-by: Martin Raspaud <[email protected]> + + +- Merge branch 'feature-vectorize' into geoloc. [Martin Raspaud] + +- Vectorize the days function. [Martin Raspaud] + +- Merge branch 'master' into geoloc. [Martin Raspaud] + +- Merge branch 'pre-master' into geoloc. [Martin Raspaud] + +- Cosmetics. [Martin Raspaud] + +- Computations for true nadir. [Martin Raspaud] + +- Bugfix in the example and added attitude correction (roll and pitch + for now). [Martin Raspaud] + + Signed-off-by: Martin Raspaud <[email protected]> + + +- Cosmetic, be consistent in name og time argument as 'utc_time' [Lars + Orum Rasmussen] + +- Get_zenith_overpass replaced by Martin's get_next_passes. [Lars Orum + Rasmussen] + +- Add sun_earth_distance_correction function. [Martin Raspaud] + +v0.2.3 (2013-03-07) +------------------- + +- Merge branch 'release-0.2.3' [Martin Raspaud] + +- Merge branch 'pre-master' into release-0.2.3. [Martin Raspaud] + +- Bumped up version number. [Martin Raspaud] + +- Corrected search for previous an_time with a substracted 10 min. dt. + [Esben S. Nielsen] + +- Merge branch 'release-0.2.2' [Martin Raspaud] + +- Import with_statement in test_aiaa.py for python 2.5 compliance. + [Esben S. Nielsen] + +- Made unit tests python 2.5 and 2.6 compliant. [Esben S. Nielsen] + +- Removed download URL from setup.py. [Esben S. Nielsen] + +- Bumped version number and marked as stable. [Esben S. Nielsen] + +- Better handling of time deltas in test_aiaa.py. [Esben S. Nielsen] + +- Updated equator test with position check. [Esben S. Nielsen] + +- Now uses nodal period for orbit number calculation instead of revs/day + for mean motion. [Esben S. Nielsen] + +- Orbit number now handles epoch AN mis-match. Made AIAA unit test path + agnostic. [Esben S. Nielsen] + +- Better __main__ [Lars Orum Rasmussen] + +- Adding risetime and falltime functions, and improving the + get_zenith_overpass function. [Adam Dybbroe] + +- Editorial. [Adam Dybbroe] + +- Cleanup. [Martin Raspaud] + + Signed-off-by: Martin Raspaud <[email protected]> + + +- Feature: Correcting/adding test cases from the aiaa. [Martin Raspaud] + +- Style: raises NotImplementedErrors instead of just Exceptions. [Martin + Raspaud] + +- Merge branch 'pre-master' of github.com:mraspaud/pyorbital into pre- + master. [Martin Raspaud] + +- Adding new function get_zenith_overpass to get the time when the + satellite passes over zenith relative to an observer on ground. [Adam + Dybbroe] + +- Feature: Added checksum for tle lines. [Martin Raspaud] + +v0.2.1 (2012-06-01) +------------------- + +- Updated version number. [Martin Raspaud] + +- Added pyorbital path to doc/source/conf.py. [Esben S. Nielsen] + +- Updated docs and added license and manifest. [Esben S. Nielsen] + +- Merge branch 'pre-master' of https://github.com/mraspaud/pyorbital + into pre-master. [Adam Dybbroe] + +- Merge branch 'pre-master' of https://github.com/mraspaud/pyorbital + into pre-master. [Lars Orum Rasmussen] + +- Added access to line1 and line2 in a Tle instance. [Lars Orum + Rasmussen] + + Change satellite to platform + + +- Spelling error. [Adam Dybbroe] + +v0.2.0 (2012-05-14) +------------------- + +- Prepared for pypi. [Martin Raspaud] + +- Merge branch 'geoloc' into pre-master. [Martin Raspaud] + +- Added now compute pixels on the ellipsoid, not on the sphere anymore. + [Martin Raspaud] + +- Merge branch 'master' into geoloc. [Martin Raspaud] + +- Updated the geoloc todo list. [Martin Raspaud] + +- Added the geoloc module. [Martin Raspaud] + +- Merge branch 'master' into pre-master. [Martin Raspaud] + + Conflicts: + pyorbital/tlefile.py + + +- Corrected handling of mean motion and orbitnumber fields in + tlefiles.py. [Esben S. Nielsen] + +- Testing getting the orbit number from the TLEs. [Adam.Dybbroe] + +- Fixing bug in tle file reading, so that also NPP and other satellites + with orbit numbers less than 9999 can be handled. [Adam.Dybbroe] + +- Typo. [Adam.Dybbroe] + +- Merge branch 'master' into pre-master. [Martin Raspaud] + +- Removed html submodule. [Martin Raspaud] + +- Fixing bug in function sun_zenith_angle. Changing interfaces so that + all public functions expects lon,lat in degrees. All internal + functions us radians. Made the lsmt and local_hour_angle functions + private. [Adam.Dybbroe] + +- Adding main. [Adam.Dybbroe] + +- Gathering unit tests to the tests-directory. [Adam.Dybbroe] + +- Added separate test-script for astronomy.py. [Adam.Dybbroe] + +- Collected all unit test scripts under the tests directory. + [Adam.Dybbroe] + +- Merge branch 'release-0.2.0' [Martin Raspaud] + + Conflicts: + doc/build + setup.py + + +- Bumped version number to 0.2.0. [Martin Raspaud] + +- Added html documentation. [Martin Raspaud] + +- Corrected sgp4's propagate in the case of array as input, and cleaned + up. [Martin Raspaud] + +- Fixed calling test_aiaa from another directory. [Martin Raspaud] + +- Vectorize merge. [Martin Raspaud] + +- Merging master branch. [Martin Raspaud] + +- Remove html submodule. [Martin Raspaud] + +- Remove html submodule. [Martin Raspaud] + +- Added Esben in the author field. [Martin Raspaud] + +- Removed unneded .pyc file. [Martin Raspaud] + +- Added unittests. [Esben S. Nielsen] + +- Corrected observer_look function and added first unittest. [Esben S. + Nielsen] + +- Corrected observer_pos in astronomy. [Esben S. Nielsen] + +- Setting up documentation. [Martin Raspaud] + +v0.1.0 (2011-10-03) +------------------- + +- Merge branch 'release-0.1.0' [Martin Raspaud] + +- Bumped version number to 0.1.0. [Martin Raspaud] + +- Merge branch 'dundee_port' into pre-master. [Martin Raspaud] + +- Cleanup and documentation. [Martin Raspaud] + +- Now using unittest module for aiaa test cases. [Martin Raspaud] + +- Added licences, and removed prints. [Martin Raspaud] + +- Added basic tests to pyorbital. [Martin Raspaud] + +- Ported SGP4 code to Dundee implementation. [Esben S. Nielsen] + +- Ported sgp4 init. [Esben S. Nielsen] + +- Added the first unit test :) [Martin Raspaud] + +- New gmst function (from AIAA paper). Cleaning. [Martin Raspaud] + +- Merged DMI and SMHI versions. [Esben S. Nielsen] + +- Made the package more package-like. [Martin Raspaud] + +- Cleanup of astronomy file. [Martin Raspaud] + +- Added a readme file. [Martin Raspaud] + +- Added astronomy.py file. [Martin Raspaud] + + diff --git a/doc/source/conf.py b/doc/source/conf.py index 5a15784..b753217 100644 --- a/doc/source/conf.py +++ b/doc/source/conf.py @@ -11,7 +11,8 @@ # All configuration values have a default; values that are commented out # serve to show the default. -import sys, os +import sys +import os # If extensions (or modules to document with autodoc) are in another directory, # add these directories to sys.path here. If the directory is relative to the @@ -19,7 +20,7 @@ import sys, os #sys.path.insert(0, os.path.abspath('.')) sys.path.insert(0, os.path.abspath('../../')) sys.path.insert(0, os.path.abspath('../../pyorbital')) -from pyorbital.version import __version__, __major__, __minor__ +from pyorbital.version import __version__ # -- General configuration ----------------------------------------------------- @@ -44,16 +45,16 @@ master_doc = 'index' # General information about the project. project = u'pyorbital' -copyright = u'2012-2014, The Pytroll crew' +copyright = u'2012-2015, The Pytroll crew' # The version info for the project you're documenting, acts as replacement for # |version| and |release|, also used in various other places throughout the # built documents. # -# The short X.Y version. -version = "v" + ".".join([__major__, __minor__]) # The full version, including alpha/beta/rc tags. release = __version__ +# The short X.Y version. +version = ".".join(release.split(".")[:2]) # The language for content autogenerated by Sphinx. Refer to documentation # for a list of supported languages. @@ -181,8 +182,8 @@ htmlhelp_basename = 'pyorbitaldoc' # Grouping the document tree into LaTeX files. List of tuples # (source start file, target name, title, author, documentclass [howto/manual]). latex_documents = [ - ('index', 'pyorbital.tex', u'pyorbital Documentation', - u'The Pytroll crew', 'manual'), + ('index', 'pyorbital.tex', u'pyorbital Documentation', + u'The Pytroll crew', 'manual'), ] # The name of an image file (relative to this directory) to place at the top of diff --git a/etc/platforms.txt b/etc/platforms.txt new file mode 100644 index 0000000..6fc5250 --- /dev/null +++ b/etc/platforms.txt @@ -0,0 +1,60 @@ +# The platform numbers are given in a file $PPP_CONFIG_DIR/platforms.txt +# in the following format. Copy this file to $PPP_CONFIG_DIR +# +# Mappings between satellite catalogue numbers and corresponding +# platform names from OSCAR. +ALOS-2 39766 +CloudSat 29107 +CryoSat-2 36508 +CSK-1 31598 +CSK-2 32376 +CSK-3 33412 +CSK-4 37216 +DMSP-F15 25991 +DMSP-F16 28054 +DMSP-F17 29522 +DMSP-F18 35951 +DMSP-F19 39630 +EOS-Aqua 27424 +EOS-Aura 28376 +EOS-Terra 25994 +FY-2D 29640 +FY-2E 33463 +FY-2F 38049 +FY-2G 40367 +FY-3A 32958 +FY-3B 37214 +FY-3C 39260 +GOES-13 29155 +GOES-14 35491 +GOES-15 36411 +Himawari-6 28622 +Himawari-7 28937 +Himawari-8 40267 +INSAT-3A 27714 +INSAT-3C 27298 +INSAT-3D 39216 +JASON-2 33105 +Kalpana-1 27525 +Landsat-7 25682 +Landsat-8 39084 +Meteosat-7 24932 +Meteosat-8 27509 +Meteosat-9 28912 +Meteosat-10 38552 +Metop-A 29499 +Metop-B 38771 +NOAA-15 25338 +NOAA-16 26536 +NOAA-17 27453 +NOAA-18 28654 +NOAA-19 33591 +RadarSat-2 32382 +Sentinel-1A 39634 +SMOS 36036 +SPOT-5 27421 +SPOT-6 38755 +SPOT-7 40053 +Suomi-NPP 37849 +TanDEM-X 36605 +TerraSAR-X 31698 diff --git a/pyorbital/geoloc.py b/pyorbital/geoloc.py index c543e6a..7ec4932 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. +# Copyright (c) 2011, 2012, 2013, 2014, 2015. # Author(s): @@ -34,44 +34,46 @@ from numpy import cos, sin, sqrt from datetime import timedelta from pyorbital.orbital import Orbital -a = 6378.137 # km -b = 6356.75231414 # km, GRS80 -#b = 6356.752314245 # km, WGS84 +a = 6378.137 # km +b = 6356.75231414 # km, GRS80 +# b = 6356.752314245 # km, WGS84 + def geodetic_lat(point, a=a, b=b): x, y, z = point - r = np.sqrt(x*x + y*y) + r = np.sqrt(x * x + y * y) geoc_lat = np.arctan2(z, r) geod_lat = geoc_lat - e2 = (a*a - b*b) / (a*a) + 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 / sqrt(1 - e2 * sin(phi)**2) + geod_lat = np.arctan2(z + a * C * e2 * 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*. """ x, y, z = query_point - r = sqrt(x*x + y*y) + r = sqrt(x * x + y * y) lat = geodetic_lat(query_point) lon = np.arctan2(y, x) - e2_ = (a*a - b*b) / (a*a) + 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) + nz_ = (1 - e2_) * n__ * sin(lat) return np.vstack([nx_, ny_, nz_]) - + class ScanGeometry(object): + """Description of the geometry of an instrument. *fovs* is the x and y viewing angles of the instrument. y is zero if the we @@ -95,8 +97,7 @@ class ScanGeometry(object): vectors. Returns vectors as stacked rows. """ # TODO: yaw steering mode ! - - + # Fake nadir: This is the intersection point between the satellite # 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 @@ -105,7 +106,7 @@ class ScanGeometry(object): nadir = subpoint(-pos) nadir /= vnorm(nadir) - + # x is along track (roll) x = vel / vnorm(vel) @@ -124,6 +125,7 @@ class ScanGeometry(object): tds = [timedelta(seconds=i) for i in self._times] return np.array(tds) + start_of_scan + class Quaternion(object): def __init__(self, scalar, vector): @@ -135,19 +137,20 @@ class Quaternion(object): zero = np.zeros_like(x) return np.array( ((w**2 + x**2 - y**2 - z**2, - 2*x*y + 2*z*w, - 2*x*z - 2*y*w, + 2 * x * y + 2 * z * w, + 2 * x * z - 2 * y * w, zero), - (2*x*y - 2*z*w, + (2 * x * y - 2 * z * w, w**2 - x**2 + y**2 - z**2, - 2*y*z + 2*x*w, + 2 * y * z + 2 * x * w, zero), - (2*x*z + 2*y*w, - 2*y*z - 2*x*w, + (2 * x * z + 2 * y * w, + 2 * y * z - 2 * x * w, w**2 - x**2 - y**2 + z**2, zero), (zero, zero, zero, w**2 + x**2 + y**2 + z**2))) + def qrotate(vector, axis, angle): """Rotate *vector* around *axis* by *angle* (in radians). @@ -155,25 +158,25 @@ def qrotate(vector, axis, angle): This function uses quaternion rotation. """ n_axis = axis / vnorm(axis) - sin_angle = np.expand_dims(sin(angle/2), 0) - if np.rank(n_axis)==1: + sin_angle = np.expand_dims(sin(angle / 2), 0) + if np.rank(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(cos(angle / 2), p__) 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. +# 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 * + def get_lonlatalt(pos, utc_time): """Calculate sublon, sublat and altitude of satellite, considering the earth an ellipsoid. @@ -185,7 +188,7 @@ def get_lonlatalt(pos, 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) + lon = np.where(lon <= -np.pi, lon + np.pi * 2, lon) r = np.sqrt(pos_x ** 2 + pos_y ** 2) lat = np.arctan2(pos_z, r) @@ -193,15 +196,17 @@ def get_lonlatalt(pos, utc_time): while True: lat2 = lat - c = 1/(np.sqrt(1 - e2 * (np.sin(lat2) ** 2))) - lat = np.arctan2(pos_z + c * e2 *np.sin(lat2), r) + c = 1 / (np.sqrt(1 - e2 * (np.sin(lat2) ** 2))) + lat = np.arctan2(pos_z + c * e2 * np.sin(lat2), r) if np.all(abs(lat - lat2) < 1e-10): break - alt = r / np.cos(lat)- c + alt = r / np.cos(lat) - c alt *= A return np.rad2deg(lon), np.rad2deg(lat), alt -### END OF DIRTY STUFF +# END OF DIRTY STUFF + + def compute_pixels((tle1, tle2), sgeom, times, rpy=(0.0, 0.0, 0.0)): """Compute cartesian coordinates of the pixels in instrument scan. """ @@ -213,18 +218,17 @@ def compute_pixels((tle1, tle2), sgeom, times, rpy=(0.0, 0.0, 0.0)): # now, get the vectors pointing to each pixel vectors = sgeom.vectors(pos, vel, *rpy) - ## compute intersection of lines (directed by vectors and passing through - ## (0, 0, 0)) and ellipsoid. Derived from: - ## http://en.wikipedia.org/wiki/Line%E2%80%93sphere_intersection - + # compute intersection of lines (directed by vectors and passing through + # (0, 0, 0)) and ellipsoid. Derived from: + # http://en.wikipedia.org/wiki/Line%E2%80%93sphere_intersection # do the computation between line and ellipsoid (WGS 84) # NB: AAPP uses GRS 80... centre = -pos - a__ = 6378.137 # km - #b__ = 6356.75231414 # km, GRS80 - b__ = 6356.752314245 # km, WGS84 - radius = np.array([[1/a__, 1/a__, 1/b__]]).T + a__ = 6378.137 # km + # 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 ldotc = np.einsum("ij,ij->j", xr_, cr_) @@ -233,14 +237,14 @@ def compute_pixels((tle1, tle2), 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 - + def norm(v): return np.sqrt(np.dot(v, v.conj())) + def mnorm(m, axis=None): """norm of a matrix of vectors stacked along the *axis* dimension. """ @@ -248,35 +252,36 @@ def mnorm(m, axis=None): axis = np.rank(m) - 1 return np.sqrt((m**2).sum(axis)) + def vnorm(m): """norms of a matrix of column vectors. """ return np.sqrt((m**2).sum(0)) + + def hnorm(m): """norms of a matrix of row vectors. """ 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 - #2 28654 99.0096 235.8581 0014859 135.4286 224.8087 14.11526826329313 - - + # NOAA 18 (from the 2011-10-12, 16:55 utc) + # 1 28654U 05018A 11284.35271227 .00000478 00000-0 28778-3 0 9246 + # 2 28654 99.0096 235.8581 0014859 135.4286 224.8087 14.11526826329313 + noaa18_tle1 = "1 28654U 05018A 11284.35271227 .00000478 00000-0 28778-3 0 9246" noaa18_tle2 = "2 28654 99.0096 235.8581 0014859 135.4286 224.8087 14.11526826329313" from datetime import datetime t = datetime(2011, 10, 12, 13, 45) - ## edge and centre of an avhrr scanline - #sgeom = ScanGeometry([(-0.9664123687741623, 0), + # edge and centre of an avhrr scanline + # sgeom = ScanGeometry([(-0.9664123687741623, 0), # (0, 0)], # [0, 0.0, ]) - #print compute_pixels((noaa18_tle1, noaa18_tle2), sgeom, t) - + # print compute_pixels((noaa18_tle1, noaa18_tle2), sgeom, t) - ## avhrr swath + # avhrr swath scanline_nb = 1 # building the avhrr angles, 2048 pixels from +55.37 to -55.37 degrees diff --git a/pyorbital/geoloc_instrument_definitions.py b/pyorbital/geoloc_instrument_definitions.py index c7070c4..60bac0d 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 Martin Raspaud +# Copyright (c) 2013, 2014, 2015 Martin Raspaud # Author(s): @@ -39,16 +39,15 @@ Both scan angles and scan times are then combined into a ScanGeometry object. import numpy as np from pyorbital.geoloc import ScanGeometry -# number of instrument scans to use. -scans_nb = 10 - ################################################################ # # AVHRR # ################################################################ -def avhrr(scans_nb, scan_points, scan_angle=55.37, decimate=1): + +def avhrr(scans_nb, scan_points, + scan_angle=55.37, frequency=1 / 6.0): """Definition of the avhrr instrument. Source: NOAA KLM User's Guide, Appendix J @@ -61,41 +60,71 @@ def avhrr(scans_nb, scan_points, scan_angle=55.37, decimate=1): avhrr_inst = np.tile(avhrr_inst, [scans_nb, 1]) # building the corresponding times array - offset = np.arange(scans_nb) * decimate / 6.0 - #times = (np.tile(scan_points * 0.000025 + 0.0025415, [scans_nb, 1]) + offset = np.arange(scans_nb) * frequency + # 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]) + np.expand_dims(offset, 1)) return ScanGeometry(avhrr_inst, times.ravel()) -################################################################ -#### avhrr, all pixels -# we take all pixels -scan_points = np.arange(2048) +def avhrr_gac(scan_times, scan_points, + scan_angle=55.37, frequency=0.5): + """Definition of the avhrr instrument, gac version -# build the scan geometry object -avhrr_all_geom = avhrr(scans_nb, scan_points) + Source: NOAA KLM User's Guide, Appendix J + http://www.ncdc.noaa.gov/oa/pod-guide/ncdc/docs/klm/html/j/app-j.htm + """ + try: + offset = np.array([(t - scan_times[0]).seconds + + (t - scan_times[0]).microseconds / 1000000.0 for t in scan_times]) + 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]) + + # 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()) ################################################################ -#### avhrr, edge pixels +# avhrr, all pixels -# we take only edge pixels -scan_points = np.array([0, 2047]) +# build the scan geometry object -# build the scan geometry object -avhrr_edge_geom = avhrr(scans_nb, scan_points) +def avhrr_all_geom(scans_nb): + # we take all pixels + scan_points = np.arange(2048) + return avhrr(scans_nb, scan_points) ################################################################ -#### avhrr, every 40th pixel from the 24th (aapp style) +# avhrr, edge pixels + +# build the scan geometry object + -# we take only every 40th pixel -scan_points = np.arange(24, 2048, 40) +def avhrr_edge_geom(scans_nb): + # we take only edge pixels + scan_points = np.array([0, 2047]) + return avhrr(scans_nb, scan_points) + +################################################################ +# avhrr, every 40th pixel from the 24th (aapp style) # build the scan geometry object -avhrr_40_geom = avhrr(scans_nb, scan_points) + + +def avhrr_40_geom(scans_nb): + # we take only every 40th pixel + scan_points = np.arange(24, 2048, 40) + return avhrr(scans_nb, scan_points) ################################################################ # @@ -103,6 +132,7 @@ avhrr_40_geom = avhrr(scans_nb, scan_points) # ################################################################ + def viirs(scans_nb, scan_indices=slice(0, None)): """Describe VIIRS instrument geometry, I-band. @@ -110,18 +140,18 @@ def viirs(scans_nb, scan_indices=slice(0, None)): entire_width = np.arange(6400) scan_points = entire_width[scan_indices] - + across_track = (scan_points / 3199.5 - 1) * np.deg2rad(-55.84) - y_max_angle = np.arctan2(11.87/2, 824.0) + y_max_angle = np.arctan2(11.87 / 2, 824.0) along_track = np.array([-y_max_angle, 0, y_max_angle]) scan_pixels = len(scan_points) scan = np.vstack((np.tile(across_track, scan_pixels), np.repeat(along_track, 6400))).T - + npp = np.tile(scan, [scans_nb, 1]) - + # from the timestamp in the filenames, a granule takes 1:25.400 to record # (85.4 seconds) so 1.779166667 would be the duration of 1 scanline # dividing the duration of a single scan by a width of 6400 pixels results @@ -149,23 +179,23 @@ def viirs(scans_nb, scan_indices=slice(0, None)): def amsua(scans_nb, edges_only=False): """ Describe AMSU-A instrument geometry - + Parameters: scans_nb | int - number of scan lines - + Keywords: * edges_only - use only edge pixels Returns: pyorbital.geoloc.ScanGeometry object - + """ - scan_len = 30 # 30 samples per scan - scan_rate = 8 # single scan, seconds - scan_angle = -48.3 # swath, degrees - sampling_interval = 0.2 # single view, seconds - sync_time = 0.00355 # delay before the actual scan starts + scan_len = 30 # 30 samples per scan + scan_rate = 8 # single scan, seconds + scan_angle = -48.3 # swath, degrees + sampling_interval = 0.2 # single view, seconds + sync_time = 0.00355 # delay before the actual scan starts if edges_only: scan_points = np.array([0, scan_len - 1]) @@ -173,7 +203,7 @@ def amsua(scans_nb, edges_only=False): scan_points = np.arange(0, scan_len) # build the instrument (scan angles) - samples = np.vstack(((scan_points / (scan_len*0.5-0.5) - 1) + samples = np.vstack(((scan_points / (scan_len * 0.5 - 0.5) - 1) * np.deg2rad(scan_angle), np.zeros((len(scan_points),)))).transpose() samples = np.tile(samples, [scans_nb, 1]) @@ -181,7 +211,7 @@ def amsua(scans_nb, edges_only=False): # building the corresponding times array offset = np.arange(scans_nb) * scan_rate times = (np.tile(scan_points * sampling_interval + sync_time, [scans_nb, 1]) - + np.expand_dims(offset, 1)) + + np.expand_dims(offset, 1)) # build the scan geometry object return ScanGeometry(samples, times.ravel()) diff --git a/pyorbital/orbital.py b/pyorbital/orbital.py index 5409579..41b0201 100644 --- a/pyorbital/orbital.py +++ b/pyorbital/orbital.py @@ -31,9 +31,9 @@ from pyorbital import tlefile from pyorbital import astronomy import warnings -ECC_EPS = 1.0e-6 # Too low for computing further drops. +ECC_EPS = 1.0e-6 # Too low for computing further drops. ECC_LIMIT_LOW = -1.0e-3 -ECC_LIMIT_HIGH = 1.0 - ECC_EPS # Too close to 1 +ECC_LIMIT_HIGH = 1.0 - ECC_EPS # Too close to 1 ECC_ALL = 1.0e-4 EPS_COS = 1.5e-12 @@ -54,8 +54,8 @@ XMNPDA = 1440.0 AE = 1.0 SECDAY = 8.6400E4 -F = 1 / 298.257223563 # Earth flattening WGS-84 -A = 6378.137 # WGS84 Equatorial radius +F = 1 / 298.257223563 # Earth flattening WGS-84 +A = 6378.137 # WGS84 Equatorial radius SGDP4_ZERO_ECC = 0 @@ -66,11 +66,13 @@ SGDP4_NEAR_NORM = 3 KS = AE * (1.0 + S0 / XKMPER) A3OVK2 = (-XJ3 / CK2) * AE**3 + class OrbitalError(Exception): pass class Orbital(object): + """Class for orbital computations. The *satellite* parameter is the name of the satellite to work on and is @@ -138,29 +140,29 @@ class Orbital(object): return pos, vel - def get_lonlatalt(self, utc_time): """Calculate sublon, sublat and altitude of satellite. http://celestrak.com/columns/v02n03/ """ - (pos_x, pos_y, pos_z), (vel_x, vel_y, vel_z) = self.get_position(utc_time, normalize=True) + (pos_x, pos_y, pos_z), (vel_x, vel_y, vel_z) = self.get_position( + utc_time, normalize=True) 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) + lon = np.where(lon <= -np.pi, lon + np.pi * 2, lon) r = np.sqrt(pos_x ** 2 + pos_y ** 2) lat = np.arctan2(pos_z, r) e2 = F * (2 - F) while True: lat2 = lat - c = 1/(np.sqrt(1 - e2 * (np.sin(lat2) ** 2))) - lat = np.arctan2(pos_z + c * e2 *np.sin(lat2), r) + c = 1 / (np.sqrt(1 - e2 * (np.sin(lat2) ** 2))) + lat = np.arctan2(pos_z + c * e2 * np.sin(lat2), r) if np.all(abs(lat - lat2) < 1e-10): break - alt = r / np.cos(lat)- c; + alt = r / np.cos(lat) - c alt *= A return np.rad2deg(lon), np.rad2deg(lat), alt @@ -181,9 +183,10 @@ class Orbital(object): Return: (Azimuth, Elevation) """ - (pos_x, pos_y, pos_z), (vel_x, vel_y, vel_z) = self.get_position(utc_time, normalize=False) + (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) = \ - astronomy.observer_position(utc_time, lon, lat, alt) + astronomy.observer_position(utc_time, lon, lat, alt) lon = np.deg2rad(lon) lat = np.deg2rad(lat) @@ -199,9 +202,11 @@ class Orbital(object): sin_theta = np.sin(theta) cos_theta = np.cos(theta) - top_s = sin_lat * cos_theta * rx + sin_lat * sin_theta * ry - cos_lat * rz + top_s = sin_lat * cos_theta * rx + \ + sin_lat * sin_theta * ry - cos_lat * rz top_e = -sin_theta * rx + cos_theta * ry - top_z = cos_lat * cos_theta * rx + cos_lat * sin_theta * ry + sin_lat * rz + top_z = cos_lat * cos_theta * rx + \ + cos_lat * sin_theta * ry + sin_lat * rz az_ = np.arctan(-top_e / top_s) @@ -222,31 +227,31 @@ class Orbital(object): orbit_period = astronomy._days(self.orbit_elements.an_period) except AttributeError: pos_epoch, vel_epoch = self.get_position(self.tle.epoch, - normalize=False) + normalize=False) if np.abs(pos_epoch[2]) > 1 or not vel_epoch[2] > 0: # Epoch not at ascending node - self.orbit_elements.an_time = self.get_last_an_time(self.tle.epoch) + self.orbit_elements.an_time = self.get_last_an_time( + self.tle.epoch) else: # Epoch at ascending node (z < 1 km) and positive v_z self.orbit_elements.an_time = self.tle.epoch self.orbit_elements.an_period = self.orbit_elements.an_time - \ - self.get_last_an_time(self.orbit_elements.an_time - - timedelta(minutes=10)) + self.get_last_an_time(self.orbit_elements.an_time + - timedelta(minutes=10)) dt = astronomy._days(utc_time - self.orbit_elements.an_time) orbit_period = astronomy._days(self.orbit_elements.an_period) - orbit = int(self.tle.orbit + dt / orbit_period + - self.tle.mean_motion_derivative * dt**2 + - self.tle.mean_motion_sec_derivative * dt**3) + self.tle.mean_motion_derivative * dt**2 + + self.tle.mean_motion_sec_derivative * dt**3) if tbus_style: orbit += 1 return orbit - def get_next_passes(self, utc_time, length, lon, lat, alt, tol=0.001): + def get_next_passes(self, utc_time, length, lon, lat, alt, tol=0.001, horizon=0): """Calculate passes for the next hours for a given start time and a given observer. @@ -258,6 +263,7 @@ class Orbital(object): lat: Latitude of observer position on ground (float) alt: Altitude above sea-level (geoid) of observer position on ground (float) tol: precision of the result in seconds + horizon: the elevation of horizon to compute risetime and falltime. Return: [(rise-time, fall-time, max-elevation-time), ...] """ @@ -266,8 +272,10 @@ class Orbital(object): """elevation """ return self.get_observer_look(utc_time + - timedelta(minutes=minutes), - lon, lat, alt)[1] + timedelta( + minutes=np.float64(minutes)), + lon, lat, alt)[1] - horizon + def elevation_inv(minutes): """inverse of elevation """ @@ -305,7 +313,7 @@ class Orbital(object): while abs(c - a) > tol: x = b - 0.5 * (((b - a) ** 2 * (f_b - f_c) - (b - c) ** 2 * (f_b - f_a)) / - ((b - a) * (f_b - f_c) - (b - c) * (f_b - f_a))) + ((b - a) * (f_b - f_c) - (b - c) * (f_b - f_a))) f_x = fun(x) if x > b: a, b, c = b, x, c @@ -317,15 +325,16 @@ class Orbital(object): return x times = utc_time + np.array([timedelta(minutes=minutes) - for minutes in range(length * 60)]) - elev = self.get_observer_look(times, lon, lat, alt)[1] + for minutes in range(length * 60)]) + elev = self.get_observer_look(times, lon, lat, alt)[1] - horizon zcs = np.where(np.diff(np.sign(elev)))[0] res = [] risetime = None falltime = None for guess in zcs: - horizon_mins = get_root_secant(elevation, guess, guess + 1.0, tol=tol/60.0) + horizon_mins = get_root_secant( + elevation, guess, guess + 1.0, tol=tol / 60.0) horizon_time = utc_time + timedelta(minutes=horizon_mins) if elev[guess] < 0: risetime = horizon_time @@ -338,9 +347,9 @@ class Orbital(object): middle = (risemins + fallmins) / 2.0 highest = utc_time + \ timedelta(minutes=get_max_parab( - elevation_inv, - middle - 0.1, middle + 0.1, - tol=tol/60.0 + elevation_inv, + middle - 0.1, middle + 0.1, + tol=tol / 60.0 )) res += [(risetime, falltime, highest)] risetime = None @@ -366,7 +375,7 @@ class Orbital(object): nmax_iter = 100 sec_step = 0.5 - t_step = timedelta(seconds=sec_step/2.0) + t_step = timedelta(seconds=sec_step / 2.0) # Local derivative: def fprime(timex): @@ -389,9 +398,9 @@ class Orbital(object): #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])) + tx1 = tx0 - timedelta(seconds=(eps * var_scale * fpr[1])) idx = idx + 1 - #print idx, tx0, tx1, var_scale, fpr + # print idx, tx0, tx1, var_scale, fpr if abs(tx1 - utc_time) < precision and idx < 2: tx1 = tx1 + timedelta(seconds=1.0) @@ -400,7 +409,9 @@ class Orbital(object): else: return None + class OrbitElements(object): + """Class holding the orbital elements. """ @@ -413,8 +424,10 @@ class OrbitElements(object): self.mean_anomaly = np.deg2rad(tle.mean_anomaly) self.mean_motion = tle.mean_motion * (np.pi * 2 / XMNPDA) - self.mean_motion_derivative = tle.mean_motion_derivative * np.pi * 2 / XMNPDA ** 2 - self.mean_motion_sec_derivative = tle.mean_motion_sec_derivative * np.pi * 2 / XMNPDA ** 3 + self.mean_motion_derivative = tle.mean_motion_derivative * \ + np.pi * 2 / XMNPDA ** 2 + self.mean_motion_sec_derivative = tle.mean_motion_sec_derivative * \ + np.pi * 2 / XMNPDA ** 3 self.bstar = tle.bstar * AE n_0 = self.mean_motion @@ -423,14 +436,14 @@ class OrbitElements(object): i_0 = self.inclination e_0 = self.excentricity - a_1 = (k_e / n_0) ** (2.0/3) - delta_1 = ((3/2.0) * (k_2 / a_1**2) * ((3 * np.cos(i_0)**2 - 1) / - (1 - e_0**2)**(2.0/3))) + a_1 = (k_e / n_0) ** (2.0 / 3) + delta_1 = ((3 / 2.0) * (k_2 / a_1**2) * ((3 * np.cos(i_0)**2 - 1) / + (1 - e_0**2)**(2.0 / 3))) - a_0 = a_1 * (1 - delta_1/3 - delta_1**2 - (134.0/81) * delta_1**3) + a_0 = a_1 * (1 - delta_1 / 3 - delta_1**2 - (134.0 / 81) * delta_1**3) - delta_0 = ((3/2.0) * (k_2 / a_0**2) * ((3 * np.cos(i_0)**2 - 1) / - (1 - e_0**2)**(2.0/3))) + delta_0 = ((3 / 2.0) * (k_2 / a_0**2) * ((3 * np.cos(i_0)**2 - 1) / + (1 - e_0**2)**(2.0 / 3))) # original mean motion n_0pp = n_0 / (1 + delta_0) @@ -445,13 +458,14 @@ class OrbitElements(object): self.perigee = (a_0pp * (1 - e_0) / AE - AE) * XKMPER self.right_ascension_lon = (self.right_ascension - - astronomy.gmst(self.epoch)) + - astronomy.gmst(self.epoch)) if self.right_ascension_lon > np.pi: self.right_ascension_lon -= 2 * np.pi class _SGDP4(object): + """Class for the SGDP4 computations. """ @@ -497,7 +511,8 @@ class _SGDP4(object): betao = np.sqrt(betao2) temp0 = 1.5 * CK2 * self.x3thm1 / (betao * betao2) del1 = temp0 / (a1**2) - a0 = a1 * (1.0 - del1 * (1.0 / 3.0 + del1 * (1.0 + del1 * 134.0 / 81.0))) + a0 = a1 * \ + (1.0 - del1 * (1.0 / 3.0 + del1 * (1.0 + del1 * 134.0 / 81.0))) del0 = temp0 / (a0**2) self.xnodp = self.xn_0 / (1.0 + del0) self.aodp = (a0 / (1.0 - del0)) @@ -536,27 +551,27 @@ class _SGDP4(object): coef_1 = coef / psisq**3.5 self.c2 = (coef_1 * self.xnodp * (self.aodp * - (1.0 + 1.5 * etasq + eeta * (4.0 + etasq)) + - (0.75 * CK2) * tsi / psisq * self.x3thm1 * - (8.0 + 3.0 * etasq * (8.0 + etasq)))) + (1.0 + 1.5 * etasq + eeta * (4.0 + etasq)) + + (0.75 * CK2) * tsi / psisq * self.x3thm1 * + (8.0 + 3.0 * etasq * (8.0 + etasq)))) 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)))) + (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 - if self.mode == SGDP4_NEAR_NORM: self.c5 = (2.0 * coef_1 * self.aodp * betao2 * - (1.0 + 2.75 * (etasq + eeta) + eeta * etasq)) + (1.0 + 2.75 * (etasq + eeta) + eeta * etasq)) if self.eo > ECC_ALL: - self.c3 = coef * tsi * A3OVK2 * self.xnodp * AE * self.sinIO / self.eo + self.c3 = coef * tsi * A3OVK2 * \ + self.xnodp * AE * self.sinIO / self.eo self.omgcof = self.bstar * self.c3 * np.cos(self.omegao) temp1 = 3.0 * CK2 * pinvsq * self.xnodp @@ -564,18 +579,18 @@ class _SGDP4(object): temp3 = 1.25 * CK4 * pinvsq**2 * self.xnodp self.xmdot = (self.xnodp + (0.5 * temp1 * betao * self.x3thm1 + 0.0625 * - temp2 * betao * (13.0 - 78.0 * theta2 + - 137.0 * theta4))) + temp2 * betao * (13.0 - 78.0 * theta2 + + 137.0 * theta4))) x1m5th = 1.0 - 5.0 * theta2 self.omgdot = (-0.5 * temp1 * x1m5th + 0.0625 * temp2 * - (7.0 - 114.0 * theta2 + 395.0 * theta4) + - temp3 * (3.0 - 36.0 * theta2 + 49.0 * theta4)) + (7.0 - 114.0 * theta2 + 395.0 * theta4) + + temp3 * (3.0 - 36.0 * theta2 + 49.0 * theta4)) xhdot1 = -temp1 * self.cosIO self.xnodot = (xhdot1 + (0.5 * temp2 * (4.0 - 19.0 * theta2) + - 2.0 * temp3 * (3.0 - 7.0 * theta2)) * self.cosIO) + 2.0 * temp3 * (3.0 - 7.0 * theta2)) * self.cosIO) if self.eo > ECC_ALL: self.xmcof = (-(2. / 3) * AE) * coef * self.bstar / eeta @@ -585,12 +600,14 @@ class _SGDP4(object): self.xnodcf = 3.5 * betao2 * xhdot1 * self.c1 self.t2cof = 1.5 * self.c1 - # Check for possible divide-by-zero for X/(1+cos(xincl)) when calculating xlcof */ - temp0 = 1.0 + self.cosIO - if np.abs(temp0) < EPS_COS: - temp0 = np.sign(temp0) * EPS_COS + # Check for possible divide-by-zero for X/(1+cos(xincl)) when + # calculating xlcof */ + temp0 = 1.0 + self.cosIO + if np.abs(temp0) < EPS_COS: + temp0 = np.sign(temp0) * EPS_COS - self.xlcof = 0.125 * A3OVK2 * self.sinIO * (3.0 + 5.0 * self.cosIO) / temp0 + self.xlcof = 0.125 * A3OVK2 * self.sinIO * \ + (3.0 + 5.0 * self.cosIO) / temp0 self.aycof = 0.25 * A3OVK2 * self.sinIO @@ -603,11 +620,13 @@ class _SGDP4(object): self.d2 = 4.0 * self.aodp * tsi * c1sq temp0 = self.d2 * tsi * self.c1 / 3.0 self.d3 = (17.0 * self.aodp + s4) * temp0 - self.d4 = 0.5 * temp0 * self.aodp * tsi * (221.0 * self.aodp + 31.0 * s4) * self.c1 + self.d4 = 0.5 * temp0 * self.aodp * tsi * \ + (221.0 * self.aodp + 31.0 * s4) * self.c1 self.t3cof = self.d2 + 2.0 * c1sq - self.t4cof = 0.25 * (3.0 * self.d3 + self.c1 * (12.0 * self.d2 + 10.0 * c1sq)) + self.t4cof = 0.25 * \ + (3.0 * self.d3 + self.c1 * (12.0 * self.d2 + 10.0 * c1sq)) self.t5cof = (0.2 * (3.0 * self.d4 + 12.0 * self.c1 * self.d3 + 6.0 * self.d2**2 + - 15.0 * c1sq * (2.0 * self.d2 + c1sq))) + 15.0 * c1sq * (2.0 * self.d2 + c1sq))) elif self.mode == SGDP4_DEEP_NORM: raise NotImplementedError('Deep space calculations not supported') @@ -620,7 +639,7 @@ class _SGDP4(object): em = self.eo xinc = self.xincl - xmp = self.xmo + self.xmdot * ts + xmp = self.xmo + self.xmdot * ts xnode = self.xnodeo + ts * (self.xnodot + ts * self.xnodcf) omega = self.omegao + self.omgdot * ts @@ -630,25 +649,31 @@ class _SGDP4(object): raise NotImplementedError('Mode "Near-space, simplified equations"' ' not implemented') elif self.mode == SGDP4_NEAR_NORM: - delm = self.xmcof * ((1.0 + self.eta * np.cos(xmp))**3 - self.delmo) + delm = self.xmcof * \ + ((1.0 + self.eta * np.cos(xmp))**3 - self.delmo) temp0 = ts * self.omgcof + delm xmp += temp0 omega -= temp0 - tempa = 1.0 - (ts * (self.c1 + ts * (self.d2 + ts * (self.d3 + ts * self.d4)))) - tempe = self.bstar * (self.c4 * ts + self.c5 * (np.sin(xmp) - self.sinXMO)) - templ = ts * ts * (self.t2cof + ts * (self.t3cof + ts * (self.t4cof + ts * self.t5cof))) + tempa = 1.0 - \ + (ts * + (self.c1 + ts * (self.d2 + ts * (self.d3 + ts * self.d4)))) + tempe = self.bstar * \ + (self.c4 * ts + self.c5 * (np.sin(xmp) - self.sinXMO)) + templ = ts * ts * \ + (self.t2cof + ts * + (self.t3cof + ts * (self.t4cof + ts * self.t5cof))) a = self.aodp * tempa**2 e = em - tempe xl = xmp + omega + xnode + self.xnodp * templ else: - raise NotImplementedError('Deep space calculations not supported') + raise NotImplementedError('Deep space calculations not supported') if np.any(a < 1): raise Exception('Satellite crased at time %s', utc_time) elif np.any(e < ECC_LIMIT_LOW): - raise ValueError('Satellite modified eccentricity to low: %e < %e' - % (e, ECC_LIMIT_LOW)) + raise ValueError('Satellite modified eccentricity to low: %s < %e' + % (str(e[e < ECC_LIMIT_LOW]), ECC_LIMIT_LOW)) e = np.where(e < ECC_EPS, ECC_EPS, e) e = np.where(e > ECC_LIMIT_HIGH, ECC_LIMIT_HIGH, e) @@ -693,7 +718,7 @@ class _SGDP4(object): # 2nd order Newton-Raphson correction. nr = np.where(np.logical_and(i == 0, np.abs(nr) > 1.25 * maxnr), np.sign(nr) * maxnr, - f / (df + 0.5*esinE*nr)) + f / (df + 0.5 * esinE * nr)) epw += nr # Short period preliminary quantities @@ -715,7 +740,8 @@ class _SGDP4(object): # Update for short term periodics to position terms. - rk = r * (1.0 - 1.5 * temp2 * betal * self.x3thm1) + 0.5 * temp1 * self.x1mth2 * cos2u + rk = r * (1.0 - 1.5 * temp2 * betal * self.x3thm1) + \ + 0.5 * temp1 * self.x1mth2 * cos2u uk = u - 0.25 * temp2 * self.x7thm1 * sin2u xnodek = xnode + 1.5 * temp2 * self.cosIO * sin2u xinck = xinc + 1.5 * temp2 * self.cosIO * self.sinIO * cos2u @@ -725,10 +751,10 @@ class _SGDP4(object): temp0 = np.sqrt(a) temp2 = XKE / (a * temp0) - rdotk = ((XKE * temp0 * esinE * invR -temp2 * temp1 * self.x1mth2 * sin2u) * + rdotk = ((XKE * temp0 * esinE * invR - temp2 * temp1 * self.x1mth2 * sin2u) * (XKMPER / AE * XMNPDA / 86400.0)) rfdotk = ((XKE * np.sqrt(pl) * invR + temp2 * temp1 * - (self.x1mth2 * cos2u + 1.5 * self.x3thm1)) * + (self.x1mth2 * cos2u + 1.5 * self.x3thm1)) * (XKMPER / AE * XMNPDA / 86400.0)) kep['radius'] = rk * XKMPER / AE diff --git a/pyorbital/tlefile.py b/pyorbital/tlefile.py index a939c62..1d57775 100644 --- a/pyorbital/tlefile.py +++ b/pyorbital/tlefile.py @@ -1,12 +1,13 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- -# Copyright (c) 2011, 2012, 2013, 2014. +# Copyright (c) 2011, 2012, 2013, 2014, 2015. # Author(s): # Esben S. Nielsen <[email protected]> # Martin Raspaud <[email protected]> +# Panu Lahtinen <[email protected]> # 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 @@ -25,84 +26,184 @@ import logging import datetime import urllib2 +import os +import glob -tle_urls = ('http://celestrak.com/NORAD/elements/weather.txt', +TLE_URLS = ('http://celestrak.com/NORAD/elements/weather.txt', 'http://celestrak.com/NORAD/elements/resource.txt') -logger = logging.getLogger(__name__) +LOGGER = logging.getLogger(__name__) + + +def read_platform_numbers(in_upper=False, num_as_int=False): + '''Read platform numbers from $PPP_CONFIG_DIR/platforms.txt if available. + ''' + + out_dict = {} + if "PPP_CONFIG_DIR" in os.environ: + platform_file = os.path.join(os.environ["PPP_CONFIG_DIR"], + "platforms.txt") + try: + fid = open(platform_file, 'r') + except IOError: + LOGGER.error("Platform file %s not found.", platform_file) + return out_dict + for row in fid: + # skip comment lines + if not row.startswith('#'): + parts = row.split() + if in_upper: + parts[0] = parts[0].upper() + if num_as_int: + parts[1] = int(parts[1]) + out_dict[parts[0]] = parts[1] + fid.close() + + return out_dict + + +SATELLITES = read_platform_numbers(in_upper=True, num_as_int=False) +''' +The platform numbers are given in a file $PPP_CONFIG/platforms.txt +in the following format: + +# Mappings between satellite catalogue numbers and corresponding +# platform names from OSCAR. +ALOS-2 39766 +CloudSat 29107 +CryoSat-2 36508 +CSK-1 31598 +CSK-2 32376 +CSK-3 33412 +CSK-4 37216 +DMSP-F15 25991 +DMSP-F16 28054 +DMSP-F17 29522 +DMSP-F18 35951 +DMSP-F19 39630 +EOS-Aqua 27424 +EOS-Aura 28376 +EOS-Terra 25994 +FY-2D 29640 +FY-2E 33463 +FY-2F 38049 +FY-2G 40367 +FY-3A 32958 +FY-3B 37214 +FY-3C 39260 +GOES-13 29155 +GOES-14 35491 +GOES-15 36411 +Himawari-6 28622 +Himawari-7 28937 +Himawari-8 40267 +INSAT-3A 27714 +INSAT-3C 27298 +INSAT-3D 39216 +JASON-2 33105 +Kalpana-1 27525 +Landsat-7 25682 +Landsat-8 39084 +Meteosat-7 24932 +Meteosat-8 27509 +Meteosat-9 28912 +Meteosat-10 38552 +Metop-A 29499 +Metop-B 38771 +NOAA-15 25338 +NOAA-16 26536 +NOAA-17 27453 +NOAA-18 28654 +NOAA-19 33591 +RadarSat-2 32382 +Sentinel-1A 39634 +SMOS 36036 +SPOT-5 27421 +SPOT-6 38755 +SPOT-7 40053 +Suomi-NPP 37849 +TanDEM-X 36605 +TerraSAR-X 31698 +''' + def read(platform, tle_file=None, line1=None, line2=None): - """Read TLE for *satellite* from *tle_file*, from *line1* and *line2*, or - from internet if none is provided. - """ + """Read TLE for *satellite* from *tle_file*, from *line1* and *line2*, from + the newest file provided in the TLES pattern, or from internet if none is + provided. + """ return Tle(platform, tle_file=tle_file, line1=line1, line2=line2) + def fetch(destination): """fetch TLE from internet and save it to *destination*. - """ + """ with open(destination, "w") as dest: - for url in tle_urls: + for url in TLE_URLS: response = urllib2.urlopen(url) dest.write(response.read()) + class ChecksumError(Exception): + '''ChecksumError. + ''' pass class Tle(object): """Class holding TLE objects. - """ + """ def __init__(self, platform, tle_file=None, line1=None, line2=None): - platform = platform.strip().upper() + self._platform = platform.strip().upper() + self._tle_file = tle_file + self._line1 = line1 + self._line2 = line2 - if line1 is not None and line2 is not None: - tle = line1.strip() + "\n" + line2.strip() - else: - if tle_file: - urls = (tle_file,) - open_func = open - else: - logger.debug("Fetch tle from the internet.") - urls = tle_urls - open_func = urllib2.urlopen - - tle = "" - for url in urls: - fp = open_func(url) - for l0 in fp: - l1, l2 = fp.next(), fp.next() - if l0.strip() == platform: - tle = l1.strip() + "\n" + l2.strip() - break - fp.close() - if tle: - break - - if not tle: - raise AttributeError, "Found no TLE entry for '%s'" % platform + self.satnumber = None + self.classification = None + self.id_launch_year = None + self.id_launch_number = None + self.id_launch_piece = None + self.epoch_year = None + self.epoch_day = None + self.epoch = None + self.mean_motion_derivative = None + self.mean_motion_sec_derivative = None + self.bstar = None + self.ephemeris_type = None + self.element_number = None + self.inclination = None + self.right_ascension = None + self.excentricity = None + self.arg_perigee = None + self.mean_anomaly = None + self.mean_motion = None + self.orbit = None - self._platform = platform - self._line1, self._line2 = tle.split('\n') - self._checksum() self._read_tle() + self._checksum() + self._parse_tle() @property def line1(self): + '''Return first TLE line.''' return self._line1 @property def line2(self): + '''Return second TLE line.''' return self._line2 @property def platform(self): + '''Return satellite platform name.''' return self._platform def _checksum(self): """Performs the checksum for the current TLE. """ - for line in [self.line1, self.line2]: + for line in [self._line1, self._line2]: check = 0 for char in line[:-1]: if char.isdigit(): @@ -114,15 +215,69 @@ class Tle(object): raise ChecksumError(self._platform + " " + line) def _read_tle(self): + '''Read TLE data. + ''' + if self._line1 is not None and self._line2 is not None: + tle = self._line1.strip() + "\n" + self._line2.strip() + else: + if self._tle_file: + urls = (self._tle_file,) + open_func = open + elif "TLES" in os.environ: + # TODO: get the TLE file closest in time to the actual satellite + # overpass, NOT the latest! + urls = (max(glob.glob(os.environ["TLES"]), + key=os.path.getctime), ) + LOGGER.debug("Reading TLE from %s", urls[0]) + open_func = open + else: + LOGGER.debug("Fetch TLE from the internet.") + urls = TLE_URLS + open_func = urllib2.urlopen + + tle = "" + designator = "1 " + SATELLITES.get(self._platform, '') + for url in urls: + fid = open_func(url) + for l_0 in fid: + if l_0.strip() == self._platform: + l_1, l_2 = fid.next(), fid.next() + tle = l_1.strip() + "\n" + l_2.strip() + break + if(self._platform in SATELLITES and + l_0.strip().startswith(designator)): + l_1 = l_0 + l_2 = fid.next() + tle = l_1.strip() + "\n" + l_2.strip() + LOGGER.debug("Found platform %s, ID: %s", + self._platform, + SATELLITES[self._platform]) + break + fid.close() + if tle: + break + + if not tle: + raise KeyError("Found no TLE entry for '%s'" % self._platform) + + self._line1, self._line2 = tle.split('\n') + + def _parse_tle(self): + '''Parse values from TLE data. + ''' def _read_tle_decimal(rep): - if rep[0] in ["-", " "]: - val = rep[0] + "." + rep[1:-2] + "e" + rep[-2:] + '''Convert *rep* to decimal value. + ''' + if rep[0] in ["-", " ", "+"]: + digits = rep[1:-2].strip() + val = rep[0] + "." + digits + "e" + rep[-2:] else: - val = "." + rep[:-2] + "e" + rep[-2:] + digits = rep[:-2].strip() + val = "." + digits + "e" + rep[-2:] return float(val) - + self.satnumber = self._line1[2:7] self.classification = self._line1[7] self.id_launch_year = self._line1[9:11] @@ -131,11 +286,10 @@ class Tle(object): 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)) + datetime.timedelta(days=self.epoch_day - 1)) self.mean_motion_derivative = float(self._line1[33:43]) self.mean_motion_sec_derivative = _read_tle_decimal(self._line1[44:52]) - self.bstar = float(self._line1[53] + "." + self._line1[54:59] + "e" + self._line1[59:61]) - _read_tle_decimal(self._line1[53:61]) + self.bstar = _read_tle_decimal(self._line1[53:61]) try: self.ephemeris_type = int(self._line1[62]) except ValueError: @@ -151,12 +305,19 @@ class Tle(object): self.orbit = int(self._line2[63:68]) def __str__(self): - import pprint, StringIO - s = StringIO.StringIO() - d = dict(([(k, v) for k, v in self.__dict__.items() if k[0] != '_'])) - pprint.pprint(d, s) - return s.getvalue()[:-1] + import pprint + import StringIO + s_var = StringIO.StringIO() + d_var = dict(([(k, v) for k, v in + self.__dict__.items() if k[0] != '_'])) + pprint.pprint(d_var, s_var) + return s_var.getvalue()[:-1] + +def main(): + '''Main for testing TLE reading. + ''' + tle_data = read('Noaa-19') + print tle_data if __name__ == '__main__': - tle = read('noaa 19') - print tle + main() diff --git a/pyorbital/version.py b/pyorbital/version.py index db58857..f6bb660 100644 --- a/pyorbital/version.py +++ b/pyorbital/version.py @@ -1,7 +1,7 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- -# Copyright (c) 2014 Martin Raspaud +# Copyright (c) 2014, 2015 Martin Raspaud # Author(s): @@ -23,8 +23,4 @@ """Version file. """ -__major__ = "0" -__minor__ = "3" -__patch__ = "2" - -__version__ = "v" + ".".join([__major__, __minor__, __patch__]) +__version__ = "v1.0.0" diff --git a/setup.cfg b/setup.cfg new file mode 100644 index 0000000..1ce39b9 --- /dev/null +++ b/setup.cfg @@ -0,0 +1,5 @@ +[bdist_rpm] +requires=numpy +release=1 +doc_files = doc/Makefile doc/source/*.rst + diff --git a/setup.py b/setup.py index 3ce9449..6a8232d 100644 --- a/setup.py +++ b/setup.py @@ -42,7 +42,7 @@ setup(name='pyorbital', test_suite='pyorbital.tests.suite', package_dir = {'pyorbital': 'pyorbital'}, packages = ['pyorbital'], - install_requires=['numpy'], + install_requires=['numpy>=1.6.0'], zip_safe=False, ) -- Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/pkg-grass/pyorbital.git _______________________________________________ Pkg-grass-devel mailing list [email protected] http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/pkg-grass-devel
