log

LICENSE.txt

@@ -0,0 +1,19 @@
+Copyright (c) 2018 Chunxiao Li
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.

README.md

@@ -0,0 +1,73 @@
+Atmospheric model for Python. Currently, only nrlmsise00 is feasible, which is valid from altitude z = 0..1000 km.
+
+## Install
+
+```sh
+pip install pyatmos
+```
+
+## Examples
+
+```python
+from pyatmos.msise import download_sw,read_sw
+from pyatmos.atmosclasses import Coordinate
+
+# Download or update the space weather file from www.celestrak.com
+swfile = download_sw() 
+# Read the space weather data
+sw_obs_pre = read_sw(swfile) 
+
+# Test 1
+# Set a specific time and location
+t = '2015-10-05 03:00:00' # time(UTC)
+lat,lon = 25,102 # latitude and longitude [degree]
+alt = 70 # altitude [km]
+# Initialize a coordinate instance by a space-time point
+st = Coordinate(t,lat,lon,alt)
+
+para_input,para_output = st.nrlmsise00(sw_obs_pre)
+print(para_input)
+{'doy': 278, 'year': 2015, 'sec': 10800.0, 'alt': 70, 'g_lat': 25, 'g_long': 102, 'lst': 9.8, 'f107A': 150, 'f107': 150, 'ap': 4, 'ap_a': array([4, 4, 4, 4, 4, 4, 4])}
+print(para_output)
+{'d': {'He': 9100292488300570.0, 'O': 0, 'N2': 1.3439413974205876e+21, 'O2': 3.52551376755781e+20, 'AR': 1.6044163757370681e+19, 'RHO': 8.225931818480755e-05, 'H': 0, 'N': 0, 'ANM O': 0}, 't': {'TINF': 1027.3184649, 'TG': 219.9649472491653}}
+
+# Test 2
+t = '2004-07-08 10:30:50' 
+lat,lon,alt = -65,-120,100 
+st = Coordinate(t,lat,lon,alt)
+para_input,para_output = st.nrlmsise00(sw_obs_pre)
+print(para_input)
+{'doy': 190, 'year': 2004, 'sec': 37850.0, 'alt': 100, 'g_lat': -65, 'g_long': -120, 'lst': 2.5138888888888893, 'f107A': 109.0, 'f107': 79.3, 'ap': 2, 'ap_a': array([2.   , 2.   , 2.   , 2.   , 2.   , 3.125, 4.625])}
+print(para_output)
+{'d': {'He': 119477307274636.89, 'O': 4.1658304136233e+17, 'N2': 7.521248904485598e+18, 'O2': 1.7444969074975662e+18, 'AR': 7.739495767665198e+16, 'RHO': 4.584596293339505e-07, 'H': 22215754381448.5, 'N': 152814261016.3964, 'ANM O': 1.8278224834873257e-37}, 't': {'TINF': 1027.3184649, 'TG': 192.5868649143824}}
+
+# Test 3
+t = '2010-02-15 12:18:37' 
+lat,lon,alt = 85,210,500 
+st = Coordinate(t,lat,lon,alt)
+para_input,para_output = st.nrlmsise00(sw_obs_pre,'NoOxygen','Aph')
+print(para_input)
+{'doy': 46, 'year': 2010, 'sec': 44317.0, 'alt': 500, 'g_lat': 85, 'g_long': 210, 'lst': 2.310277777777779, 'f107A': 83.4, 'f107': 89.4, 'ap': 14, 'ap_a': array([14.   ,  5.   ,  7.   ,  6.   , 15.   ,  5.375,  4.   ])}
+print(para_output)
+{'d': {'He': 3314507585382.5425, 'O': 3855595951659.0874, 'N2': 19285497858.028534, 'O2': 395599656.3119481, 'AR': 146073.85956102316, 'RHO': 1.2650700238089615e-13, 'H': 171775437382.8238, 'N': 38359828672.39737, 'ANM O': 5345258193.554493}, 't': {'TINF': 776.3155804924045, 'TG': 776.3139192714452}}
+
+# Test 4
+t = '2019-08-20 23:10:59' 
+lat,lon,alt = 3,5,900 
+st = Coordinate(t,lat,lon,alt)
+para_input,para_output = st.nrlmsise00(sw_obs_pre,aphmode = 'Aph')
+print(para_input)
+{'doy': 232, 'year': 2019, 'sec': 83459.0, 'alt': 900, 'g_lat': 3, 'g_long': 5, 'lst': 23.51638888888889, 'f107A': 67.4, 'f107': 67.7, 'ap': 4, 'ap_a': array([4.   , 4.   , 3.   , 3.   , 5.   , 3.625, 3.5  ])}
+print(para_output)
+{'d': {'He': 74934329990.0412, 'O': 71368139.39199762, 'N2': 104.72048033793158, 'O2': 0.09392848471935447, 'AR': 1.3231114543012155e-07, 'RHO': 8.914971667362366e-16, 'H': 207405192640.34592, 'N': 3785341.821909535, 'ANM O': 1794317839.638502}, 't': {'TINF': 646.8157488121493, 'TG': 646.8157488108872}}
+```
+
+## Reference
+
+* [Original Fortran and C code](https://ccmc.gsfc.nasa.gov/pub/modelweb/atmospheric/msis/)
+* [MSISE-00 in Python and Matlab](https://github.com/space-physics/msise00)
+* [NRLMSISE-00 Atmosphere Model - Matlab](https://ww2.mathworks.cn/matlabcentral/fileexchange/56253-nrlmsise-00-atmosphere-model?requestedDomain=zh)
+* [NRLMSISE-00 Atmosphere Model - Aerospace Blockset](https://www.mathworks.com/help/aeroblks/nrlmsise00atmospheremodel.html?requestedDomain=)
+* [NRLMSISE-00 Atmosphere Model - CCMC](https://ccmc.gsfc.nasa.gov/modelweb/models/nrlmsise00.php)
+* [NRLMSISE-00 empirical model of the atmosphere: Statistical comparisons and scientific issues](http://onlinelibrary.wiley.com/doi/10.1029/2002JA009430/pdf)
+

pyatmos/__init__.py

@@ -0,0 +1,6 @@
+'''
+pyatmos package
+
+This package is an archive of scientific routines that can be used to
+perform the estimation of papameters for various atmosphere models.  
+'''    

pyatmos/atmosclasses/__init__.py

@@ -0,0 +1,10 @@
+'''
+pyatmos atmosclasses subpackage
+
+This subpackage defines a class that facilitates the parameter estimation for various atmosphere model 
+
+Class structure:
+
+    Coordinate
+'''
+from .coordinate import Coordinate

pyatmos/atmosclasses/coordinate.py

@@ -0,0 +1,110 @@
+from astropy.time import Time
+
+from ..msise.nrlmsise00 import nrlmsise00 
+
+class Coordinate(object):
+    '''
+    space-time coordinate class
+
+    The coordinate of this class can be initialized using the following
+    constructor method:
+
+    x = coordinate(t,lat,lon,alt)
+
+    Once initialized, each class instance defines the following class
+    attributes:
+
+    t   : time, default in UTC
+    lat : latitude, default in degrees
+    lon : longitude, default in degrees
+    alt : altitude, default in km
+
+    Each class instance provides the following methods:
+
+    nrlmsise00() : Estimate the atmosphere parameters at a specific coordinate using the nrlmsise00 model 
+    '''
+    def __init__(self,t,lat,lon,alt):
+        self.t = t
+        self.lat = lat
+        self.lon = lon
+        self.alt = alt
+
+    def __repr__(self):
+    
+        return 't = {:s} UTC\nlat = {:f} deg\nlon = {:f} deg\nalt = {:f} km'.format(self.t,self.lat,self.lon,self.alt)
+
+    def nrlmsise00(self,sw_obs_pre,omode='Oxygen',aphmode='NoAph'):
+        '''
+        Estimate the atmosphere parameters at a specific coordinate(time and location) using the nrlmsise00 model.
+
+        Usage: 
+        para_input,para_output = st.nrlmsis00_data(sw_obs_pre,[omode,aphmode])
+
+        Inputs: 
+        st -> [Coordinate class instance] It can be initialized by defining a specific set of time and location
+        sw_obs_pre -> [2d str array] space weather data
+        omode -> [str, optional, default = 'Oxygen'] If 'Oxygen', Anomalous Oxygen density will be included. If 'NoOxygen', Anomalous Oxygen density will be excluded.
+        aphmode -> [str, optional, default = 'NoAph'] If NoAph', 3-hour geomagnetic index will not be used. If Aph', 3h geomagnetic index will be used. 
+
+        Outputs: 
+        para_input -> [dictionary] parameters for time, location, solar radiation index, and geomagnetic index
+        para_output -> [dictionary] parameters for atmospheric density and temperature
+
+        Examples:
+        >>> from pyatmos.msise import download_sw,read_sw
+        >>> from pyatmos.atmosclasses import Coordinate
+
+        >>> # Download or update the space weather file from www.celestrak.com
+        >>> swfile = download_sw() 
+        >>> # Read the space weather data
+        >>> sw_obs_pre = read_sw(swfile) 
+        >>> 
+        >>> # Test 1
+        >>> # Set a specific time and location
+        >>> t = '2015-10-05 03:00:00' # time(UTC)
+        >>> lat,lon = 25,102 # latitude and longitude [degree]
+        >>> alt = 70 # altitude [km]
+        >>> # Initialize a coordinate instance by a space-time point
+        >>> st = Coordinate(t,lat,lon,alt)
+        >>> 
+        >>> para_input,para_output = st.nrlmsise00(sw_obs_pre)
+        >>> print(para_input)
+        {'doy': 278, 'year': 2015, 'sec': 10800.0, 'alt': 70, 'g_lat': 25, 'g_long': 102, 'lst': 9.8, 'f107A': 150, 'f107': 150, 'ap': 4, 'ap_a': array([4, 4, 4, 4, 4, 4, 4])}
+        >>> print(para_output)
+        >>> {'d': {'He': 9100292488300570.0, 'O': 0, 'N2': 1.3439413974205876e+21, 'O2': 3.52551376755781e+20, 'AR': 1.6044163757370681e+19, 'RHO': 8.225931818480755e-05, 'H': 0, 'N': 0, 'ANM O': 0}, 't': {'TINF': 1027.3184649, 'TG': 219.9649472491653}}
+        >>> 
+        >>> # Test 2
+        >>> t = '2004-07-08 10:30:50' 
+        >>> lat,lon,alt = -65,-120,100 
+        >>> st = Coordinate(t,lat,lon,alt)
+        >>> para_input,para_output = st.nrlmsise00(sw_obs_pre)
+        >>> print(para_input)
+        {'doy': 190, 'year': 2004, 'sec': 37850.0, 'alt': 100, 'g_lat': -65, 'g_long': -120, 'lst': 2.5138888888888893, 'f107A': 109.0, 'f107': 79.3, 'ap': 2, 'ap_a': array([2.   , 2.   , 2.   , 2.   , 2.   , 3.125, 4.625])}
+        >>> print(para_output)
+        {'d': {'He': 119477307274636.89, 'O': 4.1658304136233e+17, 'N2': 7.521248904485598e+18, 'O2': 1.7444969074975662e+18, 'AR': 7.739495767665198e+16, 'RHO': 4.584596293339505e-07, 'H': 22215754381448.5, 'N': 152814261016.3964, 'ANM O': 1.8278224834873257e-37}, 't': {'TINF': 1027.3184649, 'TG': 192.5868649143824}}
+        >>>
+        >>> # Test 3
+        >>> t = '2010-02-15 12:18:37' 
+        >>> lat,lon,alt = 85,210,500 
+        >>> st = Coordinate(t,lat,lon,alt)
+        >>> para_input,para_output = st.nrlmsise00(sw_obs_pre,'NoOxygen','Aph')
+        >>> print(para_input)
+        {'doy': 46, 'year': 2010, 'sec': 44317.0, 'alt': 500, 'g_lat': 85, 'g_long': 210, 'lst': 2.310277777777779, 'f107A': 83.4, 'f107': 89.4, 'ap': 14, 'ap_a': array([14.   ,  5.   ,  7.   ,  6.   , 15.   ,  5.375,  4.   ])}
+        >>> print(para_output)
+        {'d': {'He': 3314507585382.5425, 'O': 3855595951659.0874, 'N2': 19285497858.028534, 'O2': 395599656.3119481, 'AR': 146073.85956102316, 'RHO': 1.2650700238089615e-13, 'H': 171775437382.8238, 'N': 38359828672.39737, 'ANM O': 5345258193.554493}, 't': {'TINF': 776.3155804924045, 'TG': 776.3139192714452}}
+        >>> 
+        >>> # Test 4
+        >>> t = '2019-08-20 23:10:59' 
+        >>> lat,lon,alt = 3,5,900 
+        >>> st = Coordinate(t,lat,lon,alt)
+        >>> para_input,para_output = st.nrlmsise00(sw_obs_pre,aphmode = 'Aph')
+        >>> print(para_input)
+        {'doy': 232, 'year': 2019, 'sec': 83459.0, 'alt': 900, 'g_lat': 3, 'g_long': 5, 'lst': 23.51638888888889, 'f107A': 67.4, 'f107': 67.7, 'ap': 4, 'ap_a': array([4.   , 4.   , 3.   , 3.   , 5.   , 3.625, 3.5  ])}
+        >> print(para_output)
+        {'d': {'He': 74934329990.0412, 'O': 71368139.39199762, 'N2': 104.72048033793158, 'O2': 0.09392848471935447, 'AR': 1.3231114543012155e-07, 'RHO': 8.914971667362366e-16, 'H': 207405192640.34592, 'N': 3785341.821909535, 'ANM O': 1794317839.638502}, 't': {'TINF': 646.8157488121493, 'TG': 646.8157488108872}}
+        '''
+        para_input,para_output = nrlmsise00(Time(self.t),self.lat,self.lon,self.alt,sw_obs_pre,omode,aphmode)
+        return para_input,para_output   
+
+
+            

pyatmos/msise/__init__.py

@@ -0,0 +1,23 @@
+'''
+pyatmos msise subpackage
+
+This subpackage defines the following functions:
+
+# =================== spaceweather functions ================= #
+
+download_sw - Download or update the space weather file from www.celestrak.com
+
+read_sw     - Read the space weather file
+
+get_sw      - Extract the space weather data  
+
+# ===================== utility functions ==================== #
+
+wraplon     - Wrap a longitude in range of [0,360] to [-180,180]
+
+hms2s       - Convert hour/minute/second to seconds
+
+hms2h       - Convert hour/minute/second to hours  
+
+'''              
+from .spaceweather import download_sw,read_sw

pyatmos/msise/spaceweather.py

@@ -0,0 +1,150 @@
+# -------------------------------------------------------------------- #
+# ------------------------  space weather  --------------------------- #
+# -------------------------------------------------------------------- #
+
+'''
+pyatmos spaceweather
+
+This submodule defines the following functions:
+
+download_sw - Download or update the space weather file from www.celestrak.com
+
+read_sw     - Read the space weather file
+
+get_sw      - Extract the space weather data    
+'''
+
+import numpy as np
+from datetime import datetime,timedelta
+from os import getenv,path,makedirs,remove
+from urllib.request import urlretrieve
+
+# =================== download and update sw data  =================== #
+
+def download_sw(direc=None):
+    '''
+    Download or update the space weather file from www.celestrak.com
+
+    Usage: 
+    swfile = download_sw([direc])
+
+    Inputs: 
+    direc -> [str, optionanl, default = $HOME+'/src/sw-data/'] Directory for storing sw file
+    
+    Outputs: 
+    swfile -> [str] Path of sw file
+
+    Examples:
+    >>> swfile = download_sw()
+    Downloading the latest space weather data ... Finished
+    >>> print(swfile)
+    /Users/lichunxiao/src/sw-data/SW-All.txt
+    >>> swfile = download_sw('sw-data/')
+    Downloading the latest space weather data ... Finished
+    >>> swfile = download_sw('sw-data/')
+    The existing space weather data is already up to date
+    >>> print(swfile)
+    sw-data/SW-All.txt
+    '''
+    
+    if direc is None:
+        home = getenv('HOME')
+        direc = home + '/src/sw-data/'
+    
+    swfile = direc + 'SW-All.txt'
+    url = 'https://www.celestrak.com/SpaceData/SW-All.txt'
+
+    if not path.exists(direc): makedirs(direc)
+
+    if not path.exists(swfile):
+        print('Downloading the latest space weather data',end=' ... ')
+        urlretrieve(url, swfile)
+        print('Finished')
+
+    else:
+        modified_time = datetime.fromtimestamp(path.getmtime(swfile))
+        if datetime.now() > modified_time + timedelta(days=1):
+            remove(swfile)
+            print('Updating the space weather data',end=' ... ')
+            urlretrieve(url, swfile)    
+            print('Finished')
+        else:
+            print('The existing space weather data is already up to date')    
+    return swfile
+
+# =========================== read sw file ========================== #
+
+def read_sw(swfile):
+    '''
+    Parse and read the space weather file
+
+    Usage: 
+    sw_obs_pre = read_sw(swfile)
+
+    Inputs: 
+    swfile -> [str] Path of sw file
+    
+    Outputs: 
+    sw_obs_pre -> [2d str array] sw data
+
+    Examples:
+    >>> swfile = 'sw-data/SW-All.txt'
+    >>> sw_obs_pre = read_sw(swfile)
+    >>> print(sw_obs_pre)
+    [['2020' '01' '07' ... '72.4' '68.0' '71.0']
+    ['2020' '01' '06' ... '72.4' '68.1' '70.9']
+    ['2020' '01' '05' ... '72.4' '68.2' '70.9']
+    ...
+    ['1957' '10' '03' ... '266.3' '268.1' '232.7']
+    ['1957' '10' '02' ... '253.3' '267.4' '231.7']
+    ['1957' '10' '01' ... '269.3' '266.6' '230.9']]
+    '''
+    sw_data = open(swfile,'r').readlines()
+    SW_OBS,SW_PRE = [],[]
+    flag1 = flag2 = 0
+    for line in sw_data:
+        if line.startswith('BEGIN OBSERVED'): 
+            flag1 = 1
+            continue
+        if line.startswith('END OBSERVED'): flag1 = 0 
+        if flag1 == 1: 
+            sw_p = line.split()
+            if len(sw_p) == 30:
+                del sw_p[24]
+            elif len(sw_p) == 31: 
+                sw_p = np.delete(sw_p,[23,25]) 
+            else: 
+                sw_p = np.delete(sw_p,[23,24,25,27])
+            SW_OBS.append(sw_p)
+            
+        if line.startswith('BEGIN DAILY_PREDICTED'): 
+            flag2 = 1
+            continue    
+        if line.startswith('END DAILY_PREDICTED'): break 
+        if flag2 == 1: SW_PRE.append(line.split())    
+    SW_OBS_PRE = np.vstack((np.array(SW_OBS),np.array(SW_PRE)))   
+    # inverse sort
+    SW_OBS_PRE = np.flip(SW_OBS_PRE,0)
+    return SW_OBS_PRE
+
+# ========================== extract sw data ========================== #   
+
+def get_sw(SW_OBS_PRE,t_ymd,hour):
+    j = 0
+    for ymd in SW_OBS_PRE[:,:3]:
+        if np.array_equal(t_ymd,ymd): break
+        j+=1 
+    f107A,f107,ap = float(SW_OBS_PRE[j,27]),float(SW_OBS_PRE[j+1,26]),int(SW_OBS_PRE[j,22])
+    aph_tmp_b0 = SW_OBS_PRE[j,14:22]   
+    i = int(np.floor_divide(hour,3))
+    ap_c = aph_tmp_b0[i]
+    aph_tmp_b1 = SW_OBS_PRE[j+1,14:22]
+    aph_tmp_b2 = SW_OBS_PRE[j+2,14:22]
+    aph_tmp_b3 = SW_OBS_PRE[j+3,14:22]
+    aph_tmp = np.hstack((aph_tmp_b3,aph_tmp_b2,aph_tmp_b1,aph_tmp_b0))[::-1].astype(np.float)
+    apc_index = 7-i
+    aph_c369 = aph_tmp[apc_index:apc_index+4]
+    aph_1233 = np.average(aph_tmp[apc_index+4:apc_index+12])
+    aph_3657 = np.average(aph_tmp[apc_index+12:apc_index+20])
+    aph = np.hstack((ap,aph_c369,aph_1233,aph_3657))
+    return f107A,f107,ap,aph

pyatmos/msise/utils.py

@@ -0,0 +1,33 @@
+# ------------------------------------------------------------------- #
+# ----------------------------- utilities --------------------------- #
+# ------------------------------------------------------------------- #
+
+'''
+pyatmos utils
+
+This submodule defines the following functions:
+
+wraplon - Wrap a longitude in range of [0,360] to [-180,180]
+
+hms2s   - Convert hour/minute/second to seconds
+
+hms2h   - Convert hour/minute/second to hours        
+'''
+
+# ========================= convert position ======================== #
+
+def wraplon(lon):
+    if lon > 180:
+        lonwrap = lon - 360
+    else:
+        lonwrap = lon
+    return lonwrap 
+
+# =========================== convert time ========================== #
+
+def hms2s(h,m,s):
+    return h*3.6E3 + m*60 + s
+
+def hms2h(h,m,s):
+    return h + m/60 + s/3.6E3
+

setup.py

@@ -0,0 +1,30 @@
+import setuptools
+from setuptools import setup 
+
+setup(
+    name='pyatmos',
+    version='1.0',
+    long_description_content_type='text/markdown',
+    description='A package to estimate the atmosphere parameters',
+    long_description=open('README.md', 'rb').read().decode('utf-8'),
+    license='MIT',
+    author='Chunxiao Li',
+    author_email='lcx366@126.com',
+    url='https://github.com/lcx366/ATMOS',
+    classifiers=[
+        'Intended Audience :: Education',
+        'Intended Audience :: Science/Research',
+        'Programming Language :: Python :: 3',
+        'License :: OSI Approved :: MIT License',
+        ],
+    packages=setuptools.find_packages(),
+    package_data = {
+        'pyatmos.data': ['*.npz'],
+        },
+    install_requires=[
+        'scipy',
+        'numpy',
+        'pyshtools',
+        'astropy'
+        ],
+)