Merge pull request #18 from SWxTREC/msis2.1

Add MSIS 2.1 and update code

.flake8

@@ -1,5 +1,7 @@
 [flake8]
-exclude = .git,__pycache__,build,dist
+exclude = .git,__pycache__,build,dist,examples,tests
+max-line-length = 88
 per-file-ignores =
     # imports not at top of file
     setup.py: E402
+ignore = D100,D104

.github/workflows/build-wheels.yml

@@ -53,7 +53,7 @@ jobs:
         run: |
           python -m cibuildwheel --output-dir wheelhouse
         env:
-          CIBW_BUILD: "cp37-* cp38-* cp39-* cp310-*"
+          CIBW_BUILD: "cp38-* cp39-* cp310-* cp311-*"
 
       - uses: actions/upload-artifact@v2
         with:

.github/workflows/ci.yml

@@ -17,7 +17,7 @@ jobs:
     strategy:
       matrix:
         os: [windows-latest, ubuntu-latest, macos-latest]
-        python-version: [3.7, 3.8, 3.9, '3.10']
+        python-version: [3.8, 3.9, '3.10', '3.11-dev']
     env:
       FC: gfortran-9
       F77: gfortran-9
@@ -56,6 +56,8 @@ jobs:
         pip install -r requirements-test.txt
 
     - name: Install pymsis
+      # Enable the legacy editable install capability for now
+      # numpy distutils fails compiling the wrappers without it
       run: pip install -ve .
 
     - name: Lint with flake8

.github/workflows/download_mirror.py

@@ -1,23 +1,38 @@
 """This is only for CI downloading/testing."""
-import os
+from pathlib import Path
 import tarfile
 import urllib.request
 
+MSIS21_FILE = (
+    "https://gist.github.com/greglucas/"
+    "6a545a4ccbfdb93290a96ad13a0b6f81/"
+    "raw/3f82de72d2ea5df7b1cfe492a6b2efee55349f95/"
+    "nrlmsis2.1.tar.gz"
+)
+
+with urllib.request.urlopen(MSIS21_FILE) as stream:
+    tf = tarfile.open(fileobj=stream, mode="r|gz")
+    tf.extractall(path=Path("src/msis2.1/"))
+
 # MSIS2 Tar file
-SOURCE_FILE = ("https://gist.github.com/greglucas/"
-               "6a545a4ccbfdb93290a96ad13a0b6f81/"
-               "raw/8a853f3e1c2b324f1c4ea9e0fa1e6d073258a456/"
-               "NRLMSIS2.0.tar.gz")
+MSIS20_FILE = (
+    "https://gist.github.com/greglucas/"
+    "6a545a4ccbfdb93290a96ad13a0b6f81/"
+    "raw/8a853f3e1c2b324f1c4ea9e0fa1e6d073258a456/"
+    "NRLMSIS2.0.tar.gz"
+)
 
-with urllib.request.urlopen(SOURCE_FILE) as stream:
+with urllib.request.urlopen(MSIS20_FILE) as stream:
     tf = tarfile.open(fileobj=stream, mode="r|gz")
-    tf.extractall(path=os.path.join('src', 'msis2', ''))
+    tf.extractall(path=Path("src/msis2.0/"))
 
 # MSIS-00 fortran file
-MSIS00_FILE = ("https://gist.githubusercontent.com/greglucas/"
-               "6a545a4ccbfdb93290a96ad13a0b6f81/"
-               "raw/2c1f5d899d7b42392a6b19a041d2cc213589a5f1/"
-               "NRLMSISE-00.FOR")
+MSIS00_FILE = (
+    "https://gist.githubusercontent.com/greglucas/"
+    "6a545a4ccbfdb93290a96ad13a0b6f81/"
+    "raw/2c1f5d899d7b42392a6b19a041d2cc213589a5f1/"
+    "NRLMSISE-00.FOR"
+)
 with urllib.request.urlopen(MSIS00_FILE) as response:
-    with open(os.path.join('src', 'msis00', 'NRLMSISE-00.FOR'), 'wb') as f:
+    with open(Path("src/msis00/NRLMSISE-00.FOR"), "wb") as f:
         f.write(response.read())

.gitignore

@@ -1,12 +1,12 @@
-src/msis2/*
-!src/msis2/msis2.F90
-!src/msis2/msis2.pyf
-src/msis00/*
-!src/msis00/msis00.F90
-!src/msis00/msis00.pyf
-
-pymsis/msis2.0.parm
-pymsis/msis2f*
+src/msis**/*
+!src/msis*/.gitkeep
+src/wrappers/*
+!src/wrappers/*.pyf
+!src/wrappers/*.F90
+
+pymsis/msis*.parm
+pymsis/msis20f*
+pymsis/msis21f*
 pymsis/msis00f*
 
 *.tar.gz

MANIFEST.in

@@ -1,4 +1,8 @@
-include README.rst
-include pymsis/msis2.0.parm
+include MSIS2_LICENSE
+include tests/*
 include tools/*
-include src/**/*.pyf
+# Exclude external fortran files
+exclude src/msis**/*
+exclude src/wrappers/*.[fc]
+include src/wrappers/*.pyf
+exclude pymsis/*.parm

docs/source/conf.py

@@ -8,22 +8,24 @@
 
 # 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
-# documentation root, use os.path.abspath to make it absolute, like shown here.
+# documentation root, use resolve() to make it absolute, like shown here.
 #
-import os
+from pathlib import Path
 import sys
 from sphinx_gallery.sorting import ExampleTitleSortKey
 
-sys.path.insert(0, os.path.abspath('../../pymsis'))
+import pymsis
+
+sys.path.insert(0, Path("../../pymsis").resolve())
 
 # -- Project information -----------------------------------------------------
 
-project = 'pymsis'
-copyright = '2020, Regents of the University of Colorado'
-author = 'Greg Lucas'
+project = "pymsis"
+copyright = "2020, Regents of the University of Colorado"
+author = "Greg Lucas"
 
 # The full version, including alpha/beta/rc tags
-release = '0.1.0'
+release = pymsis.__version__
 
 
 # -- General configuration ---------------------------------------------------
@@ -32,16 +34,16 @@
 # extensions coming with Sphinx (named 'sphinx.ext.*') or your custom
 # ones.
 extensions = [
-    'sphinx.ext.autodoc',
-    'sphinx.ext.autosummary',
-    'sphinx.ext.githubpages',  # Helpful for publishing to gh-pages
-    'sphinx.ext.napoleon',
-    'matplotlib.sphinxext.plot_directive',
-    'sphinx_gallery.gen_gallery',
+    "sphinx.ext.autodoc",
+    "sphinx.ext.autosummary",
+    "sphinx.ext.githubpages",  # Helpful for publishing to gh-pages
+    "sphinx.ext.napoleon",
+    "matplotlib.sphinxext.plot_directive",
+    "sphinx_gallery.gen_gallery",
 ]
 
 # Add any paths that contain templates here, relative to this directory.
-templates_path = ['_templates']
+templates_path = ["_templates"]
 
 # List of patterns, relative to source directory, that match files and
 # directories to ignore when looking for source files.
@@ -54,26 +56,26 @@
 # The theme to use for HTML and HTML Help pages.  See the documentation for
 # a list of builtin themes.
 #
-html_theme = 'pydata_sphinx_theme'
+html_theme = "pydata_sphinx_theme"
 
 html_logo = "_static/pymsis-logo.png"
 
 html_theme_options = {
-  "github_url": "https://github.com/SWxTREC/pymsis",
+    "github_url": "https://github.com/SWxTREC/pymsis",
 }
 
 # Add any paths that contain custom static files (such as style sheets) here,
 # relative to this directory. They are copied after the builtin static files,
 # so a file named "default.css" will overwrite the builtin "default.css".
-html_static_path = ['_static']
+html_static_path = ["_static"]
 
 # Autosummary
 autosummary_generate = True
 
 # Sphinx gallery
 sphinx_gallery_conf = {
-     'examples_dirs': '../../examples',   # path to example scripts
-     'gallery_dirs': 'examples',  # path to where to save generated output
-     'matplotlib_animations': True,
-     'within_subsection_order': ExampleTitleSortKey,
+    "examples_dirs": "../../examples",  # path to example scripts
+    "gallery_dirs": "examples",  # path to where to save generated output
+    "matplotlib_animations": True,
+    "within_subsection_order": ExampleTitleSortKey,
 }

examples/plot_altitude_profiles.py

@@ -20,36 +20,48 @@
 f107 = 150
 f107a = 150
 ap = 7
-aps = [[ap]*7]
+aps = [[ap] * 7]
 
-date = np.datetime64('2003-01-01T00:00')
+date = np.datetime64("2003-01-01T00:00")
 output_midnight = msis.run(date, lon, lat, alts, f107, f107a, aps)
-date = np.datetime64('2003-01-01T12:00')
+date = np.datetime64("2003-01-01T12:00")
 output_noon = msis.run(date, lon, lat, alts, f107, f107a, aps)
 
 #  output is now of the shape (1, 1, 1, 1000, 11)
 # Get rid of the single dimensions
 output_midnight = np.squeeze(output_midnight)
 output_noon = np.squeeze(output_noon)
 
-variables = ['Total mass density', 'N2', 'O2', 'O', 'He',
-             'H', 'Ar', 'N', 'Anomalous O', 'NO', 'Temperature']
+variables = [
+    "Total mass density",
+    "N2",
+    "O2",
+    "O",
+    "He",
+    "H",
+    "Ar",
+    "N",
+    "Anomalous O",
+    "NO",
+    "Temperature",
+]
 
 _, ax = plt.subplots()
 for i, label in enumerate(variables):
-    if label in ('NO', 'Total mass density', 'Temperature'):
+    if label in ("NO", "Total mass density", "Temperature"):
         # There is currently no NO data, also ignore non-number densities
         continue
-    line, = ax.plot(output_midnight[:, i], alts, linestyle='--')
+    (line,) = ax.plot(output_midnight[:, i], alts, linestyle="--")
     ax.plot(output_noon[:, i], alts, c=line.get_color(), label=label)
 
-ax.legend(loc='upper center', bbox_to_anchor=(0.5, 1.15),
-          fancybox=True, shadow=True, ncol=4)
+ax.legend(
+    loc="upper center", bbox_to_anchor=(0.5, 1.15), fancybox=True, shadow=True, ncol=4
+)
 ax.set_title(f"Longitude: {lon}, Latitude: {lat}")
-ax.set_xscale('log')
+ax.set_xscale("log")
 ax.set_xlim(1e8, 1e18)
 ax.set_ylim(0, 1000)
-ax.set_xlabel('Number density (/m$^3$)')
-ax.set_ylabel('Altitude (km)')
+ax.set_xlabel("Number density (/m$^3$)")
+ax.set_ylabel("Altitude (km)")
 
 plt.show()

examples/plot_annual_variation.py

@@ -19,35 +19,46 @@
 f107a = 150
 ap = 7
 # One years worth of data at the 12th hour every day
-dates = np.arange('2003-01', '2004-01',
-                  dtype='datetime64[D]') + np.timedelta64(12, 'h')
+dates = np.arange("2003-01", "2004-01", dtype="datetime64[D]") + np.timedelta64(12, "h")
 ndates = len(dates)
 # (F107, F107a, ap) all need to be specified at the same length as dates
-f107s = [f107]*ndates
-f107as = [f107a]*ndates
-aps = [[ap]*7]*ndates
+f107s = [f107] * ndates
+f107as = [f107a] * ndates
+aps = [[ap] * 7] * ndates
 
 output = msis.run(dates, lon, lat, alt, f107s, f107as, aps)
 #  output is now of the shape (ndates, 1, 1, 1, 11)
 # Get rid of the single dimensions
 output = np.squeeze(output)
 
 # Lets get the percent variation from the annual mean for each variable
-variation = 100*(output/output.mean(axis=0) - 1)
+variation = 100 * (output / output.mean(axis=0) - 1)
 
-variables = ['Total mass density', 'N2', 'O2', 'O', 'He',
-             'H', 'Ar', 'N', 'Anomalous O', 'NO', 'Temperature']
+variables = [
+    "Total mass density",
+    "N2",
+    "O2",
+    "O",
+    "He",
+    "H",
+    "Ar",
+    "N",
+    "Anomalous O",
+    "NO",
+    "Temperature",
+]
 
 _, ax = plt.subplots()
 for i, label in enumerate(variables):
-    if label == 'NO':
+    if label == "NO":
         # There is currently no NO data
         continue
     ax.plot(dates, variation[:, i], label=label)
 
-ax.legend(loc='upper center', bbox_to_anchor=(0.5, 1.15),
-          fancybox=True, shadow=True, ncol=5)
+ax.legend(
+    loc="upper center", bbox_to_anchor=(0.5, 1.15), fancybox=True, shadow=True, ncol=5
+)
 ax.set_xlabel(f"Longitude: {lon}, Latitude: {lat}, Altitude: {alt} km")
-ax.set_ylabel('Difference from annual mean (%)')
+ax.set_ylabel("Difference from annual mean (%)")
 
 plt.show()

examples/plot_diurnal_variation.py

@@ -20,37 +20,49 @@
 f107a = 150
 ap = 7
 # One years worth of data at the 12th hour every day
-dates = np.arange('2003-01-01', '2003-01-02', dtype='datetime64[m]')
+dates = np.arange("2003-01-01", "2003-01-02", dtype="datetime64[m]")
 ndates = len(dates)
 # (F107, F107a, ap) all need to be specified at the same length as dates
-f107s = [f107]*ndates
-f107as = [f107a]*ndates
-aps = [[ap]*7]*ndates
+f107s = [f107] * ndates
+f107as = [f107a] * ndates
+aps = [[ap] * 7] * ndates
 
 output = msis.run(dates, lon, lat, alt, f107s, f107as, aps)
 #  output is now of the shape (ndates, 1, 1, 1, 11)
 # Get rid of the single dimensions
 output = np.squeeze(output)
 
 # Lets get the percent variation from the annual mean for each variable
-variation = 100*(output/output.mean(axis=0) - 1)
+variation = 100 * (output / output.mean(axis=0) - 1)
 
-variables = ['Total mass density', 'N2', 'O2', 'O', 'He',
-             'H', 'Ar', 'N', 'Anomalous O', 'NO', 'Temperature']
+variables = [
+    "Total mass density",
+    "N2",
+    "O2",
+    "O",
+    "He",
+    "H",
+    "Ar",
+    "N",
+    "Anomalous O",
+    "NO",
+    "Temperature",
+]
 
 _, ax = plt.subplots()
 for i, label in enumerate(variables):
-    if label == 'NO':
+    if label == "NO":
         # There is currently no NO data
         continue
     ax.plot(dates, variation[:, i], label=label)
 
-ax.legend(loc='upper center', bbox_to_anchor=(0.5, 1.15),
-          fancybox=True, shadow=True, ncol=5)
+ax.legend(
+    loc="upper center", bbox_to_anchor=(0.5, 1.15), fancybox=True, shadow=True, ncol=5
+)
 ax.set_xlabel(f"Longitude: {lon}, Latitude: {lat}, Altitude: {alt} km")
-ax.set_ylabel('Difference from the daily mean (%)')
+ax.set_ylabel("Difference from the daily mean (%)")
 ax.set_xlim(dates[0], dates[-1])
 ax.xaxis.set_major_locator(mdates.HourLocator(interval=3))
-ax.xaxis.set_major_formatter(mdates.DateFormatter('%H:%M'))
+ax.xaxis.set_major_formatter(mdates.DateFormatter("%H:%M"))
 
 plt.show()