Squashed commit of the following:

* Nhirs 148 fix thread (#119)
* Add script to run through all permutations of wind parameters
* Add readthedocs config file, update pylintrc
* Move queries to separate file
* Move definition statements to separate file
* Update tcrm-tests.yml to include Python 3.9
* Remove travis CI tests
* Update DOI badge
* Change build status badge to github actions
* Read netcdf-format track file from GA SST: GA's Scenario selection tool allows users to download a selected track file, but due to technical constraints this is in a different format. This change allows TCRM to read the modified format.
* Bugfix: getPoci returns np.nan incorrectly
* track.ncReadTrackData returns true datetime objects: track.ncReadTrackData previously returned cftime.DatetimeGregorian objects, which caused newer versions of matplotlib.dates.num2date to fail. This is because we write the tracks with units of 'hours since 1900-01-01 00:00', but matplotlib.dates uses 1970-01-01 as the epoch, and works in units of days (with no way to specify units in the num2date function).
* getPoci() returns np.nan incorrectly

.github/workflows/tcrm-pylint.yml

@@ -0,0 +1,34 @@
+name: Pylint tests for TCRM 
+
+on: 
+  push:
+    branches:  [ master, develop ]
+
+jobs:
+  build:
+    name: Pylint TCRM
+    runs-on: ubuntu-latest
+
+    steps:
+    - uses: actions/checkout@v2
+    - name: Set up Python env
+      uses: conda-incubator/[email protected]
+      with: 
+        activate-environment: tcrm
+        environment-file: tcrmenv.yml
+        python-version: 3.7
+        auto-activate-base: false
+
+    - name: Install dependencies
+      run: |
+        python -m pip install --upgrade pip
+        pip install pylint
+    - name: Analysing the code with pylint
+      run: |
+        pylint --rcfile pylintrc --fail-under=7 `find -regextype egrep -regex '(.*.py)$'` |
+        tee pylint.txt
+    - name: Upload pylint.txt as artifact
+      uses: actions/upload-artifact@v2
+      with:
+        name: pylint report
+        path: pylint.txt

.github/workflows/tcrm-tests.yml

@@ -10,17 +10,20 @@ on:
     branches: [ master, develop ]
 
 jobs:
-  Hazimp:
-    name: Test HazImp
+  TCRM:
+    name: Test TCRM
     runs-on: ubuntu-latest
+    strategy: 
+      matrix:
+        python-version: [3.7, 3.8, 3.9]
     steps:
     - uses: actions/checkout@v2
     - name: Set up environment
       uses: conda-incubator/[email protected]
       with:
         activate-environment: tcrm
         environment-file: tcrmenv.yml
-        python-version: 3.7
+        python-version: ${{ matrix.python-version }}
         auto-activate-base: false
 
     - name: Test with nose

.readthedocs.yaml

@@ -0,0 +1,22 @@
+# .readthedocs.yaml
+# Read the Docs configuration file
+# See https://docs.readthedocs.io/en/stable/config-file/v2.html for details
+
+# Required
+version: 2
+
+# Build documentation in the docs/ directory with Sphinx
+sphinx:
+  builder: html
+   configuration: conf.py
+
+# Optionally build your docs in additional formats such as PDF
+formats:
+   - pdf
+
+# Optionally set the version of Python and requirements required to build your docs
+python:
+   version: 3.7
+
+conda:
+   environment: tcrmenv.yml

Evaluate/windFieldValidation.py

@@ -0,0 +1,276 @@
+import itertools
+#import matplotlib.pyplot as plt
+#from matplotlib import cm as cmap
+import numpy as np
+
+import xarray as xr
+#import seaborn as sns
+
+from wind import windmodels
+from Utilities import metutils
+from Utilities.maputils import bearing2theta, makeGrid, meshLatLon
+from Utilities.parallel import attemptParallel
+
+
+#sns.set_style('ticks', {'image.cmap':'coolwarm'})
+#sns.set_context('poster')
+#palette = [(1, 1, 1), (0.000, 0.627, 0.235), (0.412, 0.627, 0.235), (0.663, 0.780, 0.282),
+#        (0.957, 0.812, 0.000), (0.925, 0.643, 0.016), (0.835, 0.314, 0.118),
+#        (0.780, 0.086, 0.118)]
+#cmap = sns.blend_palette(palette, as_cmap=True)
+
+def polarGridAroundEye(lon, lat, margin=2, resolution=0.02):
+    R, theta = makeGrid(lon, lat, margin, resolution)
+    return R, theta
+
+def meshGrid(lon, lat, margin=2, resolution=0.02):
+    xgrid, ygrid = meshLatLon(lon, lat, margin, resolution)
+    return xgrid, ygrid
+
+def calculateWindField(lon, lat, pEnv, pCentre, rMax, vFm, thetaFm, beta,
+                       profileType='powell', windFieldType='kepert'):
+
+    pCentre = metutils.convert(pCentre, 'hPa', 'Pa')
+    pEnv = metutils.convert(pEnv, 'hPa', 'Pa')
+    vFm = metutils.convert(vFm, 'kmh', 'mps')
+    thetaFm = bearing2theta(np.pi * thetaFm / 180.)
+    thetaMax = 70.
+    rmax = metutils.convert(rMax, 'km', 'm')
+    cls = windmodels.profile(profileType)
+    if profileType=="holland":
+        profile = cls(lat, lon, pEnv, pCentre, rmax, beta)
+    else:
+        profile = cls(lat, lon, pEnv, pCentre, rmax)
+    R, theta = polarGridAroundEye(lon, lat, 5.)
+    gradV = profile.velocity(R*1000)
+    cls = windmodels.field(windFieldType)
+    windfield = cls(profile)
+    Ux, Vy = windfield.field(R*1000, theta, vFm, thetaFm, thetaMax)
+
+    surfV = np.sqrt(Ux*Ux+Vy*Vy)*1.268 # Gust conversion factor
+    return gradV, surfV
+
+"""
+lat = np.arange(-30, -4, 2, dtype=float)
+pc = np.arange(900, 991, 5, dtype=float)
+pe = np.arange(995, 1016, dtype=float)
+rm = np.arange(10, 91, 5, dtype=float)
+vfm = np.arange(0, 51, 5, dtype=float)
+gwind = np.zeros((len(lat), len(pc), len(pe), len(rm), len(vfm)))
+swind = np.zeros((len(lat), len(pc), len(pe), len(rm), len(vfm)))
+it = np.nditer(gwind, flags=['multi_index'])
+nn = gwind.size
+print(nn)
+
+lon = 120.
+thetaFm = 70
+beta = 1.6
+profileType = "powell"
+blmodel = "kepert"
+i = 0
+
+
+
+for x in it:
+    il, ic, ip, ir, iv = it.multi_index
+    gradV, surfV = calculateWindField(lon, lat[il], pe[ip], pc[ic],
+                            rm[ir], vfm[iv], thetaFm, beta,
+                            profileType=profileType,
+                            windFieldType=blmodel)
+    gwind[it.multi_index] = np.max(gradV)
+    swind[it.multi_index] = np.max(surfV)
+    i += 1
+    print(f"{100*i/nn:0.4f} %")
+
+coords = [
+    ("latitude", lat, dict(long_name="Latitude",
+                           units="degrees_south")),
+    ("pcentre", pc, dict(long_name="Central pressure",
+                         units="hPa")),
+    ("penv", pe, dict(long_name="Environmental pressure",
+                      units="hPa")),
+    ("rmax", rm, dict(long_name="Radius to maximum winds",
+                      units="km")),
+    ("vfm", vfm, dict(long_name="Forward speed",
+                      units="km/h"))
+]
+
+dims = ["latitude", 'pcentre', 'penv', 'rmax', 'vfm']
+gattrs = {
+    "long_name": "Gradient level wind speed",
+    "profile": profileType,
+    "blmodel": blmodel,
+    "description": "maximum gradient level wind speed",
+    "units": "m s-1",
+    }
+sattrs = {
+    "long_name": "Surface wind speed",
+    "profile": profileType,
+    "blmodel": blmodel,
+    "description": "maximum 0.2-s wind gust",
+    "units": "m s-1",
+    }
+
+
+gda = xr.DataArray(gwind, dims=dims, coords=coords, attrs=gattrs)
+sda = xr.DataArray(swind, dims=dims, coords=coords, attrs=sattrs)
+ds = xr.Dataset()
+ds['gradwind'] = gda
+ds['surfwind'] = sda
+ds.to_netcdf("output.nc")
+"""
+
+def balanced(iterable):
+    """
+    Balance an iterator across processors.
+
+    This partitions the work evenly across processors. However, it
+    requires the iterator to have been generated on all processors
+    before hand. This is only some magical slicing of the iterator,
+    i.e., a poor man version of scattering.
+    """
+    P, p = MPI.COMM_WORLD.size, MPI.COMM_WORLD.rank
+    return itertools.islice(iterable, p, None, P)
+
+def run():
+    lat = np.arange(-30, -4, 2, dtype=float)
+    pc = np.arange(900, 991, 5, dtype=float)
+    pe = np.arange(995, 1016, dtype=float)
+    rm = np.arange(10, 91, 5, dtype=float)
+    vfm = np.arange(0, 51, 5, dtype=float)
+    gwind = np.zeros((len(lat), len(pc), len(pe), len(rm), len(vfm)))
+    swind = np.zeros((len(lat), len(pc), len(pe), len(rm), len(vfm)))
+    it = np.nditer(gwind, flags=['multi_index'])
+    nn = gwind.size
+    #print(nn)
+
+    lon = 120.
+    thetaFm = 70
+    beta = 1.6
+    profileType = "powell"
+    blmodel = "kepert"
+    i = 0
+
+    # Attempt to start the track generator in parallel
+    global MPI
+    MPI = attemptParallel()
+    comm = MPI.COMM_WORLD
+
+    status = MPI.Status()
+    worktag = 0
+    resulttag = 1
+    idx = [it.multi_index for x in it]
+
+    if (comm.rank == 0) and (comm.size > 1):
+        w = 0
+        p = comm.size -1
+        for d in range(1, comm.size):
+            print(w)
+            if w < len(idx):
+                comm.send(idx[w], dest=d, tag=worktag)
+                w += 1
+            else:
+                comm.send(None, dest=d, tag=worktag)
+                p = w
+
+        terminated = 0
+
+        while terminated < p:
+            try:
+                result = comm.recv(source=MPI.ANY_SOURCE, tag=MPI.ANY_TAG, status=status)
+            except Exception:
+                pass
+
+            d = status.source
+            if result:
+                gV, sV, workidx = result
+                gwind[workidx] = gV
+                swind[workidx] = sV
+                #gwind[idx[w]], swind[idx[w]] = result
+
+            if w < len(idx):
+                comm.send(idx[w], dest=d, tag=worktag)
+                w += 1
+            else:
+                comm.send(None, dest=d, tag=worktag)
+                terminated += 1
+
+    elif (comm.rank != 0) and (comm.size > 1):
+        while True:
+            workidx = comm.recv(source=0, tag=worktag, status=status)
+            if workidx is None:
+                break
+            il, ic, ip, ir, iv = workidx
+            print(f"Processing {workidx}")
+            gradV, surfV = calculateWindField(lon, lat[il], pe[ip], pc[ic],
+                                rm[ir], vfm[iv], thetaFm, beta,
+                                profileType=profileType,
+                                windFieldType=blmodel)
+            results = (np.max(np.abs(gradV)), np.max(surfV), workidx)
+            comm.send(results, dest=0, tag=resulttag)
+
+    elif (comm.rank == 0) and (comm.size == 1):
+        for x in idx:
+            il, ic, ip, ir, iv = x
+            print(lat[il], pc[ic], pe[ip], rm[ir], vfm[iv])
+            gradV, surfV = calculateWindField(lon, lat[il], pe[ip], pc[ic],
+                                              rm[ir], vfm[iv], thetaFm, beta,
+                                              profileType=profileType,
+                                              windFieldType=blmodel)
+            gwind[x] = np.max(np.abs(gradV))
+            swind[x] = np.max(surfV)
+
+    comm.barrier()
+
+    coords = [
+        ("latitude", lat, dict(long_name="Latitude",
+                            units="degrees_south")),
+        ("pcentre", pc, dict(long_name="Central pressure",
+                            units="hPa")),
+        ("penv", pe, dict(long_name="Environmental pressure",
+                        units="hPa")),
+        ("rmax", rm, dict(long_name="Radius to maximum winds",
+                        units="km")),
+        ("vfm", vfm, dict(long_name="Forward speed",
+                        units="km/h"))
+    ]
+
+    dims = ["latitude", 'pcentre', 'penv', 'rmax', 'vfm']
+    gattrs = {
+        "long_name": "Gradient level wind speed",
+        "profile": profileType,
+        "blmodel": blmodel,
+        "description": "maximum gradient level wind speed",
+        "units": "m s-1",
+        }
+    sattrs = {
+        "long_name": "Surface wind speed",
+        "profile": profileType,
+        "blmodel": blmodel,
+        "description": "maximum 0.2-s wind gust",
+        "units": "m s-1",
+        }
+
+    if comm.rank == 0:
+        gda = xr.DataArray(gwind, dims=dims, coords=coords, attrs=gattrs)
+        sda = xr.DataArray(swind, dims=dims, coords=coords, attrs=sattrs)
+        ds = xr.Dataset()
+        ds['gradwind'] = gda
+        ds['surfwind'] = sda
+        ds.to_netcdf("output.nc")
+
+    MPI.Finalize()
+
+if __name__ == '__main__':
+    print("Starting")
+    global MPI, comm
+    print("Initialiszing MPI")
+    MPI = attemptParallel()
+    #import atexit
+    #atexit.register(MPI.Finalize)
+    comm = MPI.COMM_WORLD
+
+    print("Executing run()")
+    run()
+
+    #MPI.Finalize()