From 0f2f4a686909c08ed30611afdb64b6e37097fc73 Mon Sep 17 00:00:00 2001 From: Stephane Poss Date: Mon, 3 Oct 2016 17:05:05 +0200 Subject: [PATCH] Refactor to make it more pythonic. Some more changes are possible to have easier build/deploy (using pbr). --- .gitignore | 92 +++++++++ geomag/MANIFEST.in => MANIFEST.in | 4 +- geomag/README.txt => README.txt | 10 +- geomag/{geomag => }/WMM.COF | 0 geomag/{geomag => }/WMM2010.COF | 186 +++++++++--------- geomag/{geomag => }/__init__.py | 74 +++---- geomag/geomag.py | 296 ++++++++++++++++++++++++++++ geomag/geomag/geomag.py | 309 ------------------------------ geomag/setup.py | 32 ---- geomag/test.py | 39 ++++ setup.py | 33 ++++ 11 files changed, 598 insertions(+), 477 deletions(-) create mode 100644 .gitignore rename geomag/MANIFEST.in => MANIFEST.in (95%) rename geomag/README.txt => README.txt (97%) rename geomag/{geomag => }/WMM.COF (100%) rename geomag/{geomag => }/WMM2010.COF (98%) rename geomag/{geomag => }/__init__.py (91%) create mode 100644 geomag/geomag.py delete mode 100644 geomag/geomag/geomag.py delete mode 100644 geomag/setup.py create mode 100644 geomag/test.py create mode 100644 setup.py diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..3b9c8da --- /dev/null +++ b/.gitignore @@ -0,0 +1,92 @@ +# Created by .ignore support plugin (hsz.mobi) +### Python template +# Byte-compiled / optimized / DLL files +__pycache__/ +*.py[cod] +*$py.class + +# C extensions +*.so + +# Distribution / packaging +.Python +env/ +build/ +develop-eggs/ +dist/ +downloads/ +eggs/ +.eggs/ +lib/ +lib64/ +parts/ +sdist/ +var/ +*.egg-info/ +.installed.cfg +*.egg + +# PyInstaller +# Usually these files are written by a python script from a template +# before PyInstaller builds the exe, so as to inject date/other infos into it. +*.manifest +*.spec + +# Installer logs +pip-log.txt +pip-delete-this-directory.txt + +# Unit test / coverage reports +htmlcov/ +.tox/ +.coverage +.coverage.* +.cache +nosetests.xml +coverage.xml +*,cover +.hypothesis/ + +# Translations +*.mo +*.pot + +# Django stuff: +*.log +local_settings.py + +# Flask stuff: +instance/ +.webassets-cache + +# Scrapy stuff: +.scrapy + +# Sphinx documentation +docs/_build/ + +# PyBuilder +target/ + +# IPython Notebook +.ipynb_checkpoints + +# pyenv +.python-version + +# celery beat schedule file +celerybeat-schedule + +# dotenv +.env + +# virtualenv +venv/ +ENV/ + +# Spyder project settings +.spyderproject + +# Rope project settings +.ropeproject + diff --git a/geomag/MANIFEST.in b/MANIFEST.in similarity index 95% rename from geomag/MANIFEST.in rename to MANIFEST.in index dadfdd2..1cffb43 100644 --- a/geomag/MANIFEST.in +++ b/MANIFEST.in @@ -1,2 +1,2 @@ -include geomag/*.COF -exclude geomag/.svn +include geomag/*.COF +exclude geomag/.svn diff --git a/geomag/README.txt b/README.txt similarity index 97% rename from geomag/README.txt rename to README.txt index 7276621..b591b81 100644 --- a/geomag/README.txt +++ b/README.txt @@ -1,5 +1,5 @@ -Magnetic variation/declination ------------------------------- - -Calculates magnetic variation/declination for any latitude/longitude/altitude, -for any date. Uses the NOAA National Geophysical Data Center, epoch 2015 data. +Magnetic variation/declination +------------------------------ + +Calculates magnetic variation/declination for any latitude/longitude/altitude, +for any date. Uses the NOAA National Geophysical Data Center, epoch 2015 data. diff --git a/geomag/geomag/WMM.COF b/geomag/WMM.COF similarity index 100% rename from geomag/geomag/WMM.COF rename to geomag/WMM.COF diff --git a/geomag/geomag/WMM2010.COF b/geomag/WMM2010.COF similarity index 98% rename from geomag/geomag/WMM2010.COF rename to geomag/WMM2010.COF index eb78503..6bfacd6 100644 --- a/geomag/geomag/WMM2010.COF +++ b/geomag/WMM2010.COF @@ -1,93 +1,93 @@ - 2010.0 WMM-2010 11/20/2009 - 1 0 -29496.6 0.0 11.6 0.0 - 1 1 -1586.3 4944.4 16.5 -25.9 - 2 0 -2396.6 0.0 -12.1 0.0 - 2 1 3026.1 -2707.7 -4.4 -22.5 - 2 2 1668.6 -576.1 1.9 -11.8 - 3 0 1340.1 0.0 0.4 0.0 - 3 1 -2326.2 -160.2 -4.1 7.3 - 3 2 1231.9 251.9 -2.9 -3.9 - 3 3 634.0 -536.6 -7.7 -2.6 - 4 0 912.6 0.0 -1.8 0.0 - 4 1 808.9 286.4 2.3 1.1 - 4 2 166.7 -211.2 -8.7 2.7 - 4 3 -357.1 164.3 4.6 3.9 - 4 4 89.4 -309.1 -2.1 -0.8 - 5 0 -230.9 0.0 -1.0 0.0 - 5 1 357.2 44.6 0.6 0.4 - 5 2 200.3 188.9 -1.8 1.8 - 5 3 -141.1 -118.2 -1.0 1.2 - 5 4 -163.0 0.0 0.9 4.0 - 5 5 -7.8 100.9 1.0 -0.6 - 6 0 72.8 0.0 -0.2 0.0 - 6 1 68.6 -20.8 -0.2 -0.2 - 6 2 76.0 44.1 -0.1 -2.1 - 6 3 -141.4 61.5 2.0 -0.4 - 6 4 -22.8 -66.3 -1.7 -0.6 - 6 5 13.2 3.1 -0.3 0.5 - 6 6 -77.9 55.0 1.7 0.9 - 7 0 80.5 0.0 0.1 0.0 - 7 1 -75.1 -57.9 -0.1 0.7 - 7 2 -4.7 -21.1 -0.6 0.3 - 7 3 45.3 6.5 1.3 -0.1 - 7 4 13.9 24.9 0.4 -0.1 - 7 5 10.4 7.0 0.3 -0.8 - 7 6 1.7 -27.7 -0.7 -0.3 - 7 7 4.9 -3.3 0.6 0.3 - 8 0 24.4 0.0 -0.1 0.0 - 8 1 8.1 11.0 0.1 -0.1 - 8 2 -14.5 -20.0 -0.6 0.2 - 8 3 -5.6 11.9 0.2 0.4 - 8 4 -19.3 -17.4 -0.2 0.4 - 8 5 11.5 16.7 0.3 0.1 - 8 6 10.9 7.0 0.3 -0.1 - 8 7 -14.1 -10.8 -0.6 0.4 - 8 8 -3.7 1.7 0.2 0.3 - 9 0 5.4 0.0 -0.0 0.0 - 9 1 9.4 -20.5 -0.1 -0.0 - 9 2 3.4 11.5 0.0 -0.2 - 9 3 -5.2 12.8 0.3 0.0 - 9 4 3.1 -7.2 -0.4 -0.1 - 9 5 -12.4 -7.4 -0.3 0.1 - 9 6 -0.7 8.0 0.1 -0.0 - 9 7 8.4 2.1 -0.1 -0.2 - 9 8 -8.5 -6.1 -0.4 0.3 - 9 9 -10.1 7.0 -0.2 0.2 - 10 0 -2.0 0.0 0.0 0.0 - 10 1 -6.3 2.8 -0.0 0.1 - 10 2 0.9 -0.1 -0.1 -0.1 - 10 3 -1.1 4.7 0.2 0.0 - 10 4 -0.2 4.4 -0.0 -0.1 - 10 5 2.5 -7.2 -0.1 -0.1 - 10 6 -0.3 -1.0 -0.2 -0.0 - 10 7 2.2 -3.9 0.0 -0.1 - 10 8 3.1 -2.0 -0.1 -0.2 - 10 9 -1.0 -2.0 -0.2 0.0 - 10 10 -2.8 -8.3 -0.2 -0.1 - 11 0 3.0 0.0 0.0 0.0 - 11 1 -1.5 0.2 0.0 -0.0 - 11 2 -2.1 1.7 -0.0 0.1 - 11 3 1.7 -0.6 0.1 0.0 - 11 4 -0.5 -1.8 -0.0 0.1 - 11 5 0.5 0.9 0.0 0.0 - 11 6 -0.8 -0.4 -0.0 0.1 - 11 7 0.4 -2.5 -0.0 0.0 - 11 8 1.8 -1.3 -0.0 -0.1 - 11 9 0.1 -2.1 0.0 -0.1 - 11 10 0.7 -1.9 -0.1 -0.0 - 11 11 3.8 -1.8 -0.0 -0.1 - 12 0 -2.2 0.0 -0.0 0.0 - 12 1 -0.2 -0.9 0.0 -0.0 - 12 2 0.3 0.3 0.1 0.0 - 12 3 1.0 2.1 0.1 -0.0 - 12 4 -0.6 -2.5 -0.1 0.0 - 12 5 0.9 0.5 -0.0 -0.0 - 12 6 -0.1 0.6 0.0 0.1 - 12 7 0.5 -0.0 0.0 0.0 - 12 8 -0.4 0.1 -0.0 0.0 - 12 9 -0.4 0.3 0.0 -0.0 - 12 10 0.2 -0.9 0.0 -0.0 - 12 11 -0.8 -0.2 -0.1 0.0 - 12 12 0.0 0.9 0.1 0.0 -999999999999999999999999999999999999999999999999 -999999999999999999999999999999999999999999999999 + 2010.0 WMM-2010 11/20/2009 + 1 0 -29496.6 0.0 11.6 0.0 + 1 1 -1586.3 4944.4 16.5 -25.9 + 2 0 -2396.6 0.0 -12.1 0.0 + 2 1 3026.1 -2707.7 -4.4 -22.5 + 2 2 1668.6 -576.1 1.9 -11.8 + 3 0 1340.1 0.0 0.4 0.0 + 3 1 -2326.2 -160.2 -4.1 7.3 + 3 2 1231.9 251.9 -2.9 -3.9 + 3 3 634.0 -536.6 -7.7 -2.6 + 4 0 912.6 0.0 -1.8 0.0 + 4 1 808.9 286.4 2.3 1.1 + 4 2 166.7 -211.2 -8.7 2.7 + 4 3 -357.1 164.3 4.6 3.9 + 4 4 89.4 -309.1 -2.1 -0.8 + 5 0 -230.9 0.0 -1.0 0.0 + 5 1 357.2 44.6 0.6 0.4 + 5 2 200.3 188.9 -1.8 1.8 + 5 3 -141.1 -118.2 -1.0 1.2 + 5 4 -163.0 0.0 0.9 4.0 + 5 5 -7.8 100.9 1.0 -0.6 + 6 0 72.8 0.0 -0.2 0.0 + 6 1 68.6 -20.8 -0.2 -0.2 + 6 2 76.0 44.1 -0.1 -2.1 + 6 3 -141.4 61.5 2.0 -0.4 + 6 4 -22.8 -66.3 -1.7 -0.6 + 6 5 13.2 3.1 -0.3 0.5 + 6 6 -77.9 55.0 1.7 0.9 + 7 0 80.5 0.0 0.1 0.0 + 7 1 -75.1 -57.9 -0.1 0.7 + 7 2 -4.7 -21.1 -0.6 0.3 + 7 3 45.3 6.5 1.3 -0.1 + 7 4 13.9 24.9 0.4 -0.1 + 7 5 10.4 7.0 0.3 -0.8 + 7 6 1.7 -27.7 -0.7 -0.3 + 7 7 4.9 -3.3 0.6 0.3 + 8 0 24.4 0.0 -0.1 0.0 + 8 1 8.1 11.0 0.1 -0.1 + 8 2 -14.5 -20.0 -0.6 0.2 + 8 3 -5.6 11.9 0.2 0.4 + 8 4 -19.3 -17.4 -0.2 0.4 + 8 5 11.5 16.7 0.3 0.1 + 8 6 10.9 7.0 0.3 -0.1 + 8 7 -14.1 -10.8 -0.6 0.4 + 8 8 -3.7 1.7 0.2 0.3 + 9 0 5.4 0.0 -0.0 0.0 + 9 1 9.4 -20.5 -0.1 -0.0 + 9 2 3.4 11.5 0.0 -0.2 + 9 3 -5.2 12.8 0.3 0.0 + 9 4 3.1 -7.2 -0.4 -0.1 + 9 5 -12.4 -7.4 -0.3 0.1 + 9 6 -0.7 8.0 0.1 -0.0 + 9 7 8.4 2.1 -0.1 -0.2 + 9 8 -8.5 -6.1 -0.4 0.3 + 9 9 -10.1 7.0 -0.2 0.2 + 10 0 -2.0 0.0 0.0 0.0 + 10 1 -6.3 2.8 -0.0 0.1 + 10 2 0.9 -0.1 -0.1 -0.1 + 10 3 -1.1 4.7 0.2 0.0 + 10 4 -0.2 4.4 -0.0 -0.1 + 10 5 2.5 -7.2 -0.1 -0.1 + 10 6 -0.3 -1.0 -0.2 -0.0 + 10 7 2.2 -3.9 0.0 -0.1 + 10 8 3.1 -2.0 -0.1 -0.2 + 10 9 -1.0 -2.0 -0.2 0.0 + 10 10 -2.8 -8.3 -0.2 -0.1 + 11 0 3.0 0.0 0.0 0.0 + 11 1 -1.5 0.2 0.0 -0.0 + 11 2 -2.1 1.7 -0.0 0.1 + 11 3 1.7 -0.6 0.1 0.0 + 11 4 -0.5 -1.8 -0.0 0.1 + 11 5 0.5 0.9 0.0 0.0 + 11 6 -0.8 -0.4 -0.0 0.1 + 11 7 0.4 -2.5 -0.0 0.0 + 11 8 1.8 -1.3 -0.0 -0.1 + 11 9 0.1 -2.1 0.0 -0.1 + 11 10 0.7 -1.9 -0.1 -0.0 + 11 11 3.8 -1.8 -0.0 -0.1 + 12 0 -2.2 0.0 -0.0 0.0 + 12 1 -0.2 -0.9 0.0 -0.0 + 12 2 0.3 0.3 0.1 0.0 + 12 3 1.0 2.1 0.1 -0.0 + 12 4 -0.6 -2.5 -0.1 0.0 + 12 5 0.9 0.5 -0.0 -0.0 + 12 6 -0.1 0.6 0.0 0.1 + 12 7 0.5 -0.0 0.0 0.0 + 12 8 -0.4 0.1 -0.0 0.0 + 12 9 -0.4 0.3 0.0 -0.0 + 12 10 0.2 -0.9 0.0 -0.0 + 12 11 -0.8 -0.2 -0.1 0.0 + 12 12 0.0 0.9 0.1 0.0 +999999999999999999999999999999999999999999999999 +999999999999999999999999999999999999999999999999 diff --git a/geomag/geomag/__init__.py b/geomag/__init__.py similarity index 91% rename from geomag/geomag/__init__.py rename to geomag/__init__.py index c210ed9..4011b71 100644 --- a/geomag/geomag/__init__.py +++ b/geomag/__init__.py @@ -1,36 +1,38 @@ -"""geomag package -by Christopher Weiss cmweiss@gmail.com - -Adapted from the geomagc software and World Magnetic Model of the NOAA -Satellite and Information Service, National Geophysical Data Center -http://www.ngdc.noaa.gov/geomag/WMM/DoDWMM.shtml - -Suggestions for improvements are appreciated. - -USAGE: ->>> import geomag ->>> geomag.declination(80,0) --3.382344140520556 -""" - -from . import geomag - -__singleton__ = geomag.GeoMag() - -def declination(*args, **kargs): - """Calculate magnetic declination in degrees - dlat = latitude in degrees - dlon = longitude in degrees - h = altitude in feet, default=0 - time = date for computing declination, default=today - """ - mag = __singleton__.GeoMag(*args, **kargs) - return mag.dec - -def mag_heading(hdg, *args, **kargs): - """Calculates the magnetic heading from a true heading. - hdg = true heading in degrees - All other parameters are the same as declination. - """ - dec = declination(*args, **kargs) - return (hdg - dec + 360.0) % 360 +"""geomag package +by Christopher Weiss cmweiss@gmail.com + +Adapted from the geomagc software and World Magnetic Model of the NOAA +Satellite and Information Service, National Geophysical Data Center +http://www.ngdc.noaa.gov/geomag/WMM/DoDWMM.shtml + +Suggestions for improvements are appreciated. + +USAGE: +>>> import geomag +>>> geomag.declination(80,0) +-3.382344140520556 +""" + +from geomag import * + +__singleton__ = GeoMag() + + +def declination(*args, **kargs): + """Calculate magnetic declination in degrees + dlat = latitude in degrees + dlon = longitude in degrees + h = altitude in feet, default=0 + time = date for computing declination, default=today + """ + mag = __singleton__.GeoMag(*args, **kargs) + return mag.dec + + +def mag_heading(hdg, *args, **kargs): + """Calculates the magnetic heading from a true heading. + hdg = true heading in degrees + All other parameters are the same as declination. + """ + dec = declination(*args, **kargs) + return (hdg - dec + 360.0) % 360 diff --git a/geomag/geomag.py b/geomag/geomag.py new file mode 100644 index 0000000..0249d32 --- /dev/null +++ b/geomag/geomag.py @@ -0,0 +1,296 @@ +# -*- coding: utf-8 -*- +""" + +geomag.py +by Christopher Weiss cmweiss@gmail.com +Improved by Stephane Poss stephane.poss@airnavigation.aero + +Adapted from the geomagc software and World Magnetic Model of the NOAA +Satellite and Information Service, National Geophysical Data Center +http://www.ngdc.noaa.gov/geomag/WMM/DoDWMM.shtml + +Suggestions for improvements are appreciated. + +USAGE: +>>> gm = GeoMag("WMM2010.COF") +>>> mag = gm.GeoMag(80,0) +>>> mag.dec +-3.6269555737620465 + +""" + +import math +import os +from datetime import date + +WMM_FILE = os.path.join(os.path.dirname(__file__), 'WMM.COF') + + +class GeoMag(object): + def GeoMag(self, dlat, dlon, h=0, + time=date.today()): # latitude (decimal degrees), longitude (decimal degrees), altitude (feet), date + # time = date('Y') + date('z')/365 + time = time.year + ((time - date(time.year, 1, 1)).days / 365.0) + alt = h / 3280.8399 + + otime = oalt = olat = olon = -1000.0 + + dt = time - self.epoch + glat = dlat + glon = dlon + rlat = math.radians(glat) + rlon = math.radians(glon) + srlon = math.sin(rlon) + srlat = math.sin(rlat) + crlon = math.cos(rlon) + crlat = math.cos(rlat) + srlat2 = srlat * srlat + crlat2 = crlat * crlat + self.sp[1] = srlon + self.cp[1] = crlon + + # /* CONVERT FROM GEODETIC COORDS. TO SPHERICAL COORDS. */ + if alt != oalt or glat != olat: + q = math.sqrt(self.a2 - self.c2 * srlat2) + q1 = alt * q + q2 = ((q1 + self.a2) / (q1 + self.b2)) * ((q1 + self.a2) / (q1 + self.b2)) + ct = srlat / math.sqrt(q2 * crlat2 + srlat2) + st = math.sqrt(1.0 - (ct * ct)) + r2 = (alt * alt) + 2.0 * q1 + (self.a4 - self.c4 * srlat2) / (q * q) + r = math.sqrt(r2) + d = math.sqrt(self.a2 * crlat2 + self.b2 * srlat2) + ca = (alt + d) / r + sa = self.c2 * crlat * srlat / (r * d) + + if glon != olon: + for m in range(2, self.maxord + 1): + self.sp[m] = self.sp[1] * self.cp[m - 1] + self.cp[1] * self.sp[m - 1] + self.cp[m] = self.cp[1] * self.cp[m - 1] - self.sp[1] * self.sp[m - 1] + + aor = self.re / r + ar = aor * aor + br = bt = bp = bpp = 0.0 + for n in range(1, self.maxord + 1): + ar *= aor + + # for (m=0,D3=1,D4=(n+m+D3)/D3;D4>0;D4--,m+=D3): + m = 0 + D3 = 1 + # D4=(n+m+D3)/D3 + D4 = (n + m + 1) + while D4 > 0: + + # /* + # COMPUTE UNNORMALIZED ASSOCIATED LEGENDRE POLYNOMIALS + # AND DERIVATIVES VIA RECURSION RELATIONS + # */ + if alt != oalt or glat != olat: + if n == m: + self.p[m][n] = st * self.p[m - 1][n - 1] + self.dp[m][n] = st * self.dp[m - 1][n - 1] + ct * self.p[m - 1][n - 1] + + elif n == 1 and m == 0: + self.p[m][n] = ct * self.p[m][n - 1] + self.dp[m][n] = ct * self.dp[m][n - 1] - st * self.p[m][n - 1] + + elif n > 1 and n != m: + if m > n - 2: + self.p[m][n - 2] = 0 + if m > n - 2: + self.dp[m][n - 2] = 0.0 + self.p[m][n] = ct * self.p[m][n - 1] - self.k[m][n] * self.p[m][n - 2] + self.dp[m][n] = ct * self.dp[m][n - 1] - st * self.p[m][n - 1] - self.k[m][n] * self.dp[m][ + n - 2] + + # /* + # TIME ADJUST THE GAUSS COEFFICIENTS + # */ + if time != otime: + self.tc[m][n] = self.c[m][n] + dt * self.cd[m][n] + if m != 0: + self.tc[n][m - 1] = self.c[n][m - 1] + dt * self.cd[n][m - 1] + + # /* + # ACCUMULATE TERMS OF THE SPHERICAL HARMONIC EXPANSIONS + # */ + par = ar * self.p[m][n] + + if m == 0: + temp1 = self.tc[m][n] * self.cp[m] + temp2 = self.tc[m][n] * self.sp[m] + else: + temp1 = self.tc[m][n] * self.cp[m] + self.tc[n][m - 1] * self.sp[m] + temp2 = self.tc[m][n] * self.sp[m] - self.tc[n][m - 1] * self.cp[m] + + bt -= ar * temp1 * self.dp[m][n] + bp += self.fm[m] * temp2 * par + br += self.fn[n] * temp1 * par + # /* + # SPECIAL CASE: NORTH/SOUTH GEOGRAPHIC POLES + # */ + if st == 0.0 and m == 1: + if n == 1: + self.pp[n] = self.pp[n - 1] + else: + self.pp[n] = ct * self.pp[n - 1] - self.k[m][n] * self.pp[n - 2] + parp = ar * self.pp[n] + bpp += self.fm[m] * temp2 * parp + + D4 -= 1 + m += 1 + + if st == 0.0: + bp = bpp + else: + bp /= st + # /* + # ROTATE MAGNETIC VECTOR COMPONENTS FROM SPHERICAL TO + # GEODETIC COORDINATES + # */ + bx = -bt * ca - br * sa + by = bp + bz = bt * sa - br * ca + # /* + # COMPUTE DECLINATION (DEC), INCLINATION (DIP) AND + # TOTAL INTENSITY (TI) + # */ + bh = math.sqrt((bx * bx) + (by * by)) + ti = math.sqrt((bh * bh) + (bz * bz)) + dec = math.degrees(math.atan2(by, bx)) + dip = math.degrees(math.atan2(bz, bh)) + # /* + # COMPUTE MAGNETIC GRID VARIATION IF THE CURRENT + # GEODETIC POSITION IS IN THE ARCTIC OR ANTARCTIC + # (I.E. GLAT > +55 DEGREES OR GLAT < -55 DEGREES) + + # OTHERWISE, SET MAGNETIC GRID VARIATION TO -999.0 + # */ + gv = -999.0 + if math.fabs(glat) >= 55.: + if glat > 0.0 and glon >= 0.0: + gv = dec - glon + if glat > 0.0 > glon: + gv = dec + math.fabs(glon) + if glat < 0.0 <= glon: + gv = dec + glon + if glat < 0.0 and glon < 0.0: + gv = dec - math.fabs(glon) + if gv > +180.0: + gv -= 360.0 + if gv < -180.0: + gv += 360.0 + + otime = time + oalt = alt + olat = glat + olon = glon + + class RetObj: + pass + + retobj = RetObj() + retobj.dec = dec + retobj.dip = dip + retobj.ti = ti + retobj.bh = bh + retobj.bx = bx + retobj.by = by + retobj.bz = bz + retobj.lat = dlat + retobj.lon = dlon + retobj.alt = h + retobj.time = time + + return retobj + + def __init__(self, wmm_filename=None): + if not wmm_filename: + wmm_filename = WMM_FILE + wmm = [] + with open(wmm_filename) as wmm_file: + for line in wmm_file: + linevals = line.strip().split() + if len(linevals) == 3: + self.epoch = float(linevals[0]) + self.model = linevals[1] + self.modeldate = linevals[2] + elif len(linevals) == 6: + linedict = {'n': int(float(linevals[0])), + 'm': int(float(linevals[1])), + 'gnm': float(linevals[2]), + 'hnm': float(linevals[3]), + 'dgnm': float(linevals[4]), + 'dhnm': float(linevals[5])} + wmm.append(linedict) + + z = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] + self.maxord = self.maxdeg = 12 + self.tc = [z[0:13], z[0:13], z[0:13], z[0:13], z[0:13], z[0:13], z[0:13], z[0:13], z[0:13], z[0:13], z[0:13], + z[0:13], z[0:13], z[0:13]] + self.sp = z[0:14] + self.cp = z[0:14] + self.cp[0] = 1.0 + self.pp = z[0:13] + self.pp[0] = 1.0 + self.p = [z[0:14], z[0:14], z[0:14], z[0:14], z[0:14], z[0:14], z[0:14], z[0:14], z[0:14], z[0:14], z[0:14], + z[0:14], z[0:14], z[0:14]] + self.p[0][0] = 1.0 + self.dp = [z[0:13], z[0:13], z[0:13], z[0:13], z[0:13], z[0:13], z[0:13], z[0:13], z[0:13], z[0:13], z[0:13], + z[0:13], z[0:13], z[0:13]] + self.a = 6378.137 + self.b = 6356.7523142 + self.re = 6371.2 + self.a2 = self.a * self.a + self.b2 = self.b * self.b + self.c2 = self.a2 - self.b2 + self.a4 = self.a2 * self.a2 + self.b4 = self.b2 * self.b2 + self.c4 = self.a4 - self.b4 + + self.c = [z[0:14], z[0:14], z[0:14], z[0:14], z[0:14], z[0:14], z[0:14], z[0:14], z[0:14], z[0:14], z[0:14], + z[0:14], z[0:14], z[0:14]] + self.cd = [z[0:14], z[0:14], z[0:14], z[0:14], z[0:14], z[0:14], z[0:14], z[0:14], z[0:14], z[0:14], z[0:14], + z[0:14], z[0:14], z[0:14]] + + for wmmnm in wmm: + m = wmmnm['m'] + n = wmmnm['n'] + gnm = wmmnm['gnm'] + hnm = wmmnm['hnm'] + dgnm = wmmnm['dgnm'] + dhnm = wmmnm['dhnm'] + if m <= n: + self.c[m][n] = gnm + self.cd[m][n] = dgnm + if m != 0: + self.c[n][m - 1] = hnm + self.cd[n][m - 1] = dhnm + + # /* CONVERT SCHMIDT NORMALIZED GAUSS COEFFICIENTS TO UNNORMALIZED */ + self.snorm = [z[0:13], z[0:13], z[0:13], z[0:13], z[0:13], z[0:13], z[0:13], z[0:13], z[0:13], z[0:13], z[0:13], + z[0:13], z[0:13]] + self.snorm[0][0] = 1.0 + self.k = [z[0:13], z[0:13], z[0:13], z[0:13], z[0:13], z[0:13], z[0:13], z[0:13], z[0:13], z[0:13], z[0:13], + z[0:13], z[0:13]] + self.k[1][1] = 0.0 + self.fn = [0.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0] + self.fm = [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0] + for n in range(1, self.maxord + 1): + self.snorm[0][n] = self.snorm[0][n - 1] * (2.0 * n - 1) / n + j = 2.0 + # for (m=0,D1=1,D2=(n-m+D1)/D1;D2>0;D2--,m+=D1): + m = 0 + D1 = 1 + D2 = (n - m + D1) / D1 + while D2 > 0: + self.k[m][n] = (((n - 1) * (n - 1)) - (m * m)) / ((2.0 * n - 1) * (2.0 * n - 3.0)) + if m > 0: + flnmj = ((n - m + 1.0) * j) / (n + m) + self.snorm[m][n] = self.snorm[m - 1][n] * math.sqrt(flnmj) + j = 1.0 + self.c[n][m - 1] = self.snorm[m][n] * self.c[n][m - 1] + self.cd[n][m - 1] = self.snorm[m][n] * self.cd[n][m - 1] + self.c[m][n] = self.snorm[m][n] * self.c[m][n] + self.cd[m][n] = self.snorm[m][n] * self.cd[m][n] + D2 -= 1 + m += D1 diff --git a/geomag/geomag/geomag.py b/geomag/geomag/geomag.py deleted file mode 100644 index d546d55..0000000 --- a/geomag/geomag/geomag.py +++ /dev/null @@ -1,309 +0,0 @@ -# geomag.py -# by Christopher Weiss cmweiss@gmail.com - -# Adapted from the geomagc software and World Magnetic Model of the NOAA -# Satellite and Information Service, National Geophysical Data Center -# http://www.ngdc.noaa.gov/geomag/WMM/DoDWMM.shtml -# -# Suggestions for improvements are appreciated. - -# USAGE: -# -# >>> gm = geomag.GeoMag("WMM.COF") -# >>> mag = gm.GeoMag(80,0) -# >>> mag.dec -# -6.1335150785195536 -# >>> - -import math, os, unittest -from datetime import date - -class GeoMag: - - def GeoMag(self, dlat, dlon, h=0, time=date.today()): # latitude (decimal degrees), longitude (decimal degrees), altitude (feet), date - #time = date('Y') + date('z')/365 - time = time.year+((time - date(time.year,1,1)).days/365.0) - alt = h/3280.8399 - - otime = oalt = olat = olon = -1000.0 - - dt = time - self.epoch - glat = dlat - glon = dlon - rlat = math.radians(glat) - rlon = math.radians(glon) - srlon = math.sin(rlon) - srlat = math.sin(rlat) - crlon = math.cos(rlon) - crlat = math.cos(rlat) - srlat2 = srlat*srlat - crlat2 = crlat*crlat - self.sp[1] = srlon - self.cp[1] = crlon - - #/* CONVERT FROM GEODETIC COORDS. TO SPHERICAL COORDS. */ - if (alt != oalt or glat != olat): - q = math.sqrt(self.a2-self.c2*srlat2) - q1 = alt*q - q2 = ((q1+self.a2)/(q1+self.b2))*((q1+self.a2)/(q1+self.b2)) - ct = srlat/math.sqrt(q2*crlat2+srlat2) - st = math.sqrt(1.0-(ct*ct)) - r2 = (alt*alt)+2.0*q1+(self.a4-self.c4*srlat2)/(q*q) - r = math.sqrt(r2) - d = math.sqrt(self.a2*crlat2+self.b2*srlat2) - ca = (alt+d)/r - sa = self.c2*crlat*srlat/(r*d) - - if (glon != olon): - for m in range(2,self.maxord+1): - self.sp[m] = self.sp[1]*self.cp[m-1]+self.cp[1]*self.sp[m-1] - self.cp[m] = self.cp[1]*self.cp[m-1]-self.sp[1]*self.sp[m-1] - - aor = self.re/r - ar = aor*aor - br = bt = bp = bpp = 0.0 - for n in range(1,self.maxord+1): - ar = ar*aor - - #for (m=0,D3=1,D4=(n+m+D3)/D3;D4>0;D4--,m+=D3): - m=0 - D3=1 - #D4=(n+m+D3)/D3 - D4=(n+m+1) - while D4>0: - - # /* - # COMPUTE UNNORMALIZED ASSOCIATED LEGENDRE POLYNOMIALS - # AND DERIVATIVES VIA RECURSION RELATIONS - # */ - if (alt != oalt or glat != olat): - if (n == m): - self.p[m][n] = st * self.p[m-1][n-1] - self.dp[m][n] = st*self.dp[m-1][n-1]+ct*self.p[m-1][n-1] - - elif (n == 1 and m == 0): - self.p[m][n] = ct*self.p[m][n-1] - self.dp[m][n] = ct*self.dp[m][n-1]-st*self.p[m][n-1] - - elif (n > 1 and n != m): - if (m > n-2): - self.p[m][n-2] = 0 - if (m > n-2): - self.dp[m][n-2] = 0.0 - self.p[m][n] = ct*self.p[m][n-1]-self.k[m][n]*self.p[m][n-2] - self.dp[m][n] = ct*self.dp[m][n-1] - st*self.p[m][n-1]-self.k[m][n]*self.dp[m][n-2] - - # /* - # TIME ADJUST THE GAUSS COEFFICIENTS - # */ - if (time != otime): - self.tc[m][n] = self.c[m][n]+dt*self.cd[m][n] - if (m != 0): - self.tc[n][m-1] = self.c[n][m-1]+dt*self.cd[n][m-1] - - # /* - # ACCUMULATE TERMS OF THE SPHERICAL HARMONIC EXPANSIONS - # */ - par = ar*self.p[m][n] - - if (m == 0): - temp1 = self.tc[m][n]*self.cp[m] - temp2 = self.tc[m][n]*self.sp[m] - else: - temp1 = self.tc[m][n]*self.cp[m]+self.tc[n][m-1]*self.sp[m] - temp2 = self.tc[m][n]*self.sp[m]-self.tc[n][m-1]*self.cp[m] - - bt = bt-ar*temp1*self.dp[m][n] - bp = bp + (self.fm[m] * temp2 * par) - br = br + (self.fn[n] * temp1 * par) - # /* - # SPECIAL CASE: NORTH/SOUTH GEOGRAPHIC POLES - # */ - if (st == 0.0 and m == 1): - if (n == 1): - self.pp[n] = self.pp[n-1] - else: - self.pp[n] = ct*self.pp[n-1]-self.k[m][n]*self.pp[n-2] - parp = ar*self.pp[n] - bpp = bpp + (self.fm[m]*temp2*parp) - - D4=D4-1 - m=m+1 - - if (st == 0.0): - bp = bpp - else: - bp = bp/st - # /* - # ROTATE MAGNETIC VECTOR COMPONENTS FROM SPHERICAL TO - # GEODETIC COORDINATES - # */ - bx = -bt*ca-br*sa - by = bp - bz = bt*sa-br*ca - # /* - # COMPUTE DECLINATION (DEC), INCLINATION (DIP) AND - # TOTAL INTENSITY (TI) - # */ - bh = math.sqrt((bx*bx)+(by*by)) - ti = math.sqrt((bh*bh)+(bz*bz)) - dec = math.degrees(math.atan2(by,bx)) - dip = math.degrees(math.atan2(bz,bh)) - # /* - # COMPUTE MAGNETIC GRID VARIATION IF THE CURRENT - # GEODETIC POSITION IS IN THE ARCTIC OR ANTARCTIC - # (I.E. GLAT > +55 DEGREES OR GLAT < -55 DEGREES) - - # OTHERWISE, SET MAGNETIC GRID VARIATION TO -999.0 - # */ - gv = -999.0 - if (math.fabs(glat) >= 55.): - if (glat > 0.0 and glon >= 0.0): - gv = dec-glon - if (glat > 0.0 and glon < 0.0): - gv = dec+math.fabs(glon); - if (glat < 0.0 and glon >= 0.0): - gv = dec+glon - if (glat < 0.0 and glon < 0.0): - gv = dec-math.fabs(glon) - if (gv > +180.0): - gv = gv - 360.0 - if (gv < -180.0): - gv = gv + 360.0 - - otime = time - oalt = alt - olat = glat - olon = glon - - class RetObj: - pass - retobj = RetObj() - retobj.dec = dec - retobj.dip = dip - retobj.ti = ti - retobj.bh = bh - retobj.bx = bx - retobj.by = by - retobj.bz = bz - retobj.lat = dlat - retobj.lon = dlon - retobj.alt = h - retobj.time = time - - return retobj - - def __init__(self, wmm_filename=None): - if not wmm_filename: - wmm_filename = os.path.join(os.path.dirname(__file__), 'WMM.COF') - wmm=[] - with open(wmm_filename) as wmm_file: - for line in wmm_file: - linevals = line.strip().split() - if len(linevals) == 3: - self.epoch = float(linevals[0]) - self.model = linevals[1] - self.modeldate = linevals[2] - elif len(linevals) == 6: - linedict = {'n': int(float(linevals[0])), - 'm': int(float(linevals[1])), - 'gnm': float(linevals[2]), - 'hnm': float(linevals[3]), - 'dgnm': float(linevals[4]), - 'dhnm': float(linevals[5])} - wmm.append(linedict) - - z = [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0] - self.maxord = self.maxdeg = 12 - self.tc = [z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13]] - self.sp = z[0:14] - self.cp = z[0:14] - self.cp[0] = 1.0 - self.pp = z[0:13] - self.pp[0] = 1.0 - self.p = [z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14]] - self.p[0][0] = 1.0 - self.dp = [z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13]] - self.a = 6378.137 - self.b = 6356.7523142 - self.re = 6371.2 - self.a2 = self.a*self.a - self.b2 = self.b*self.b - self.c2 = self.a2-self.b2 - self.a4 = self.a2*self.a2 - self.b4 = self.b2*self.b2 - self.c4 = self.a4 - self.b4 - - self.c = [z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14]] - self.cd = [z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14]] - - for wmmnm in wmm: - m = wmmnm['m'] - n = wmmnm['n'] - gnm = wmmnm['gnm'] - hnm = wmmnm['hnm'] - dgnm = wmmnm['dgnm'] - dhnm = wmmnm['dhnm'] - if (m <= n): - self.c[m][n] = gnm - self.cd[m][n] = dgnm - if (m != 0): - self.c[n][m-1] = hnm - self.cd[n][m-1] = dhnm - - #/* CONVERT SCHMIDT NORMALIZED GAUSS COEFFICIENTS TO UNNORMALIZED */ - self.snorm = [z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13]] - self.snorm[0][0] = 1.0 - self.k = [z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13]] - self.k[1][1] = 0.0 - self.fn = [0.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0,10.0,11.0,12.0,13.0] - self.fm = [0.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0,10.0,11.0,12.0] - for n in range(1,self.maxord+1): - self.snorm[0][n] = self.snorm[0][n-1]*(2.0*n-1)/n - j=2.0 - #for (m=0,D1=1,D2=(n-m+D1)/D1;D2>0;D2--,m+=D1): - m=0 - D1=1 - D2=(n-m+D1)/D1 - while (D2 > 0): - self.k[m][n] = (((n-1)*(n-1))-(m*m))/((2.0*n-1)*(2.0*n-3.0)) - if (m > 0): - flnmj = ((n-m+1.0)*j)/(n+m) - self.snorm[m][n] = self.snorm[m-1][n]*math.sqrt(flnmj) - j = 1.0 - self.c[n][m-1] = self.snorm[m][n]*self.c[n][m-1] - self.cd[n][m-1] = self.snorm[m][n]*self.cd[n][m-1] - self.c[m][n] = self.snorm[m][n]*self.c[m][n] - self.cd[m][n] = self.snorm[m][n]*self.cd[m][n] - D2=D2-1 - m=m+D1 - -class GeoMagTest(unittest.TestCase): - - d1=date(2015,1,1) - d2=date(2017,7,2) - - test_values = ( - # date, alt, lat, lon, var - (d1, 0, 80, 0, -3.85), - (d1, 0, 0, 120, 0.57), - (d1, 0, -80, 240, 69.81), - (d1, 328083.99, 80, 0, -4.27), - (d1, 328083.99, 0, 120, 0.56), - (d1, 328083.99, -80, 240, 69.22), - (d2, 0, 80, 0, -2.75), - (d2, 0, 0, 120, 0.32), - (d2, 0, -80, 240, 69.58), - (d2, 328083.99, 80, 0, -3.17), - (d2, 328083.99, 0, 120, 0.32), - (d2, 328083.99, -80, 240, 69.00), - ) - - def test_declination(self): - gm = GeoMag() - for values in self.test_values: - calcval=gm.GeoMag(values[2], values[3], values[1], values[0]) - self.assertAlmostEqual(values[4], calcval.dec, 2, 'Expected %s, result %s' % (values[4], calcval.dec)) - -if __name__ == '__main__': - unittest.main() diff --git a/geomag/setup.py b/geomag/setup.py deleted file mode 100644 index dec8558..0000000 --- a/geomag/setup.py +++ /dev/null @@ -1,32 +0,0 @@ -from distutils.core import setup -setup( - name = "geomag", - packages = ["geomag"], - #data_files = [('geomag', ('geomag/WMM.COF',))], - package_data = {'geomag': ['WMM.COF','WMM2010.COF']}, - version = "0.9.2015", - description = "Magnetic variation/declination", - author = "Christopher Weiss", - author_email = "cmweiss@gmail.com", - url = "http://geomag.googlecode.com/", - download_url = "//pypi.python.org/packages/source/g/geomag/geomag-0.9.2015.zip", - keywords = ["magnetic", "variation", "declination"], - classifiers = [ - "Programming Language :: Python", - "Development Status :: 4 - Beta", - "Environment :: Other Environment", - "Intended Audience :: Developers", - "Intended Audience :: Science/Research", - "License :: OSI Approved :: GNU Library or Lesser General Public License (LGPL)", - "Operating System :: OS Independent", - "Topic :: Scientific/Engineering :: GIS", - "Topic :: Utilities", - ], - long_description = """\ -Magnetic variation/declination ------------------------------- - -Calculates magnetic variation/declination for any latitude/longitude/altitude, -for any date. Uses the NOAA National Geophysical Data Center, epoch 2015 data. -""" -) diff --git a/geomag/test.py b/geomag/test.py new file mode 100644 index 0000000..3946235 --- /dev/null +++ b/geomag/test.py @@ -0,0 +1,39 @@ +# -*- coding: utf-8 -*- +""" +© 2012 - 2016 Xample Sàrl + +Author: Stephane Poss +Date: 03.10.16 +""" +import unittest +import datetime + +from geomag import GeoMag + + +class GeoMagTest(unittest.TestCase): + def setUp(self): + d1 = datetime.date(2015, 1, 1) + d2 = datetime.date(2017, 7, 2) + + self.test_values = ( + # date, alt, lat, lon, var + (d1, 0, 80, 0, -3.85), + (d1, 0, 0, 120, 0.57), + (d1, 0, -80, 240, 69.81), + (d1, 328083.99, 80, 0, -4.27), + (d1, 328083.99, 0, 120, 0.56), + (d1, 328083.99, -80, 240, 69.22), + (d2, 0, 80, 0, -2.75), + (d2, 0, 0, 120, 0.32), + (d2, 0, -80, 240, 69.58), + (d2, 328083.99, 80, 0, -3.17), + (d2, 328083.99, 0, 120, 0.32), + (d2, 328083.99, -80, 240, 69.00), + ) + + def test_declination(self): + gm = GeoMag() + for values in self.test_values: + calcval = gm.GeoMag(values[2], values[3], values[1], values[0]) + self.assertAlmostEqual(values[4], calcval.dec, 2, 'Expected %s, result %s' % (values[4], calcval.dec)) diff --git a/setup.py b/setup.py new file mode 100644 index 0000000..fe70b57 --- /dev/null +++ b/setup.py @@ -0,0 +1,33 @@ +from distutils.core import setup + +setup( + name="geomag", + packages=["geomag"], + # data_files = [('geomag', ('geomag/WMM.COF',))], + package_data={'geomag': ['WMM.COF', 'WMM2010.COF']}, + version="0.9.2015", + description="Magnetic variation/declination", + author="Christopher Weiss", + author_email="cmweiss@gmail.com", + url="http://geomag.googlecode.com/", + download_url="//pypi.python.org/packages/source/g/geomag/geomag-0.9.2015.zip", + keywords=["magnetic", "variation", "declination"], + classifiers=[ + "Programming Language :: Python", + "Development Status :: 4 - Beta", + "Environment :: Other Environment", + "Intended Audience :: Developers", + "Intended Audience :: Science/Research", + "License :: OSI Approved :: GNU Library or Lesser General Public License (LGPL)", + "Operating System :: OS Independent", + "Topic :: Scientific/Engineering :: GIS", + "Topic :: Utilities", + ], + long_description="""\ +Magnetic variation/declination +------------------------------ + +Calculates magnetic variation/declination for any latitude/longitude/altitude, +for any date. Uses the NOAA National Geophysical Data Center, epoch 2015 data. +""" +)