makedeclinationtable.py

Go to the documentation of this file.

#!/usr/bin/env python 
'''Create magnetic declination table from UTM coordinates in degrees
Requires `pip install docopt numpy utm`
 
Usage:
    makedeclinationtable.py <easting> <northing> <zone> <letter>
 
Arguments:
    easting   UTM Easting of the site
    northing  UTM Northing of the site
    zone      UTM zone
    letter    UTM letter
'''
from docopt import docopt
import requests
import json
import utm
 
#convert date to timestamp (float). 
#Taken from https://stackoverflow.com/questions/8777753/converting-datetime-date-to-utc-timestamp-in-python
 
#from datetime import datetime, date, timezone
import datetime
import calendar
import sys
def gettimestamp(year, month, day):
    d = datetime.date(year, month, day)
    ts = calendar.timegm(d.timetuple())
    if sys.version_info < (3,12,0):
      dt = datetime.datetime.utcfromtimestamp(ts)
      fts = dt.replace(tzinfo=datetime.timezone.utc).timestamp()
    else:
      dt = datetime.datetime.fromtimestamp(ts, datetime.UTC)
      fts = dt.timestamp()
    return int(fts)
 
args = docopt(__doc__, version='0.1')
easting =  args['<easting>'] 
northing = args['<northing>']  
zone = args['<zone>']
letter = args['<letter>']
latitude, longitude = utm.to_latlon(float(easting), float(northing), int(zone), letter)
print("#lat long",latitude,longitude)
    
model = 'IGRF' #data for 1590-2024
day = 1
key = 'zNEw7' #access key
 
for year in range(2000,2030,1):
    for month in range(1,13,1):
        request_text = 'https://www.ngdc.noaa.gov/geomag-web/calculators/calculateDeclination?lat1='+str(latitude) + \
                '&lon1='+str(longitude) + \
                '&model='+model +  \
                '&startYear='+str(year) + \
                '&startMonth='+str(month) + \
                '&startDay='+str(day) + \
                '&resultFormat=json' + \
                '&key='+key
#        print(request_text)
        x = requests.get(request_text)
#        print(x)
        y = json.loads(x.text)
        result = y['result'][0]
        declination = result['declination']
        print("(*this)[", gettimestamp(year,month, day), "] = ", declination, ";")