diff --git a/.github/workflows/tcrm-pylint.yml b/.github/workflows/tcrm-pylint.yml new file mode 100644 index 00000000..d4f1bb93 --- /dev/null +++ b/.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/setup-miniconda@v2.0.0 + 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 diff --git a/.github/workflows/tcrm-tests.yml b/.github/workflows/tcrm-tests.yml new file mode 100644 index 00000000..5f6798a3 --- /dev/null +++ b/.github/workflows/tcrm-tests.yml @@ -0,0 +1,34 @@ +# This workflow will install Python dependencies, run tests and lint with a variety of Python versions +# For more information see: https://help.github.com/actions/language-and-framework-guides/using-python-with-github-actions + +name: Unit tests for TCRM + +on: + push: + branches: [ master, develop ] + pull_request: + branches: [ master, develop ] + +jobs: + 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/setup-miniconda@v2.0.0 + with: + activate-environment: tcrm + environment-file: tcrmenv.yml + python-version: ${{ matrix.python-version }} + auto-activate-base: false + + - name: Test with pytest + env: + PYTHONPATH: ~/tcrm;~/tcrm/Utilities + shell: bash -l {0} + run: | + pytest -x --cov=. --cov-report xml diff --git a/.readthedocs.yaml b/.readthedocs.yaml new file mode 100644 index 00000000..c0c9b0d5 --- /dev/null +++ b/.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 \ No newline at end of file diff --git a/.travis.yml b/.travis.yml deleted file mode 100644 index aef70720..00000000 --- a/.travis.yml +++ /dev/null @@ -1,30 +0,0 @@ -language: shell - -env: - - PYTHON_VERSION=3.6 - - PYTHON_VERSION=3.7 - - PYTHON_VERSION=3.8 - -os: - - linux - - windows - -before_install: - - source ./preinstall.sh - -install: - - source ./postinstall.sh - -branches: - except: - - config - - notebooks - -script: - - python installer/setup.py build_ext -i - - nosetests -v --with-coverage --cover-package=. - -after_success: coveralls -notifications: - slack: - secure: Ckmwy59ytS1GPRZ5Tmvzad6+W9AzvfjNJAa4orgdKS/WktoK4b9W2rbTHxi8V3hBLIDUCso8vIQi3rVXpWY3cFMvb/uRbXO4GiIW1iua3CKjxd+dEw4E6/8DEknS1qdGJRDhN9/3ucZNvSGHY3EQQDfxb/R+OGd2jT6+jed8pss= diff --git a/Evaluate/evaluate.py b/Evaluate/evaluate.py index c306c8c1..9b8b92c6 100644 --- a/Evaluate/evaluate.py +++ b/Evaluate/evaluate.py @@ -1600,5 +1600,25 @@ def process_args(argv): log.info("Completed {0}".format(sys.argv[0])) + +def bom2tcrm(df, trackId): + """ + Transforms a dataframe in BoM format into a tcrm track. + + """ + df['Datetime'] = pd.to_datetime(df.validtime) + df['Speed'] = df.translation_speed + df['CentralPressure'] = df.pcentre + df['Longitude'] = df.longitude + df['Latitude'] = df.latitude + df['EnvPressure'] = df.poci + df['rMax'] = df.rmax + df['trackId'] = df.disturbance.values + + track = Track(df) + track.trackId = [trackId] + return track + + if __name__ == '__main__': process_args(sys.argv[1:]) diff --git a/Evaluate/interpolateTracks.py b/Evaluate/interpolateTracks.py index c1b551d6..3762c0b3 100644 --- a/Evaluate/interpolateTracks.py +++ b/Evaluate/interpolateTracks.py @@ -3,12 +3,17 @@ import numpy as np from datetime import datetime, timedelta -from matplotlib.dates import num2date +from matplotlib.dates import num2date, date2num from scipy.interpolate import interp1d, splev, splrep from Utilities.maputils import latLon2Azi from Utilities.loadData import loadTrackFile, maxWindSpeed +from Utilities.parallel import disableOnWorkers from Utilities.track import Track, ncSaveTracks +# from pycxml.pycxml import loadfile +import pandas as pd + + LOG = logging.getLogger(__name__) LOG.addHandler(logging.NullHandler()) @@ -70,7 +75,7 @@ def interpolate(track, delta, interpolation_type=None): track.Minute)] else: day_ = track.Datetime - + timestep = timedelta(delta/24.) try: time_ = np.array([d.toordinal() + (d.hour + d.minute/60.)/24.0 @@ -84,7 +89,7 @@ def interpolate(track, delta, interpolation_type=None): else: raise dt_ = 24.0 * np.diff(time_) - dt = np.empty(len(track.data), dtype=float) + dt = np.zeros(len(track.data), dtype=float) dt[1:] = dt_ # Convert all times to a time after initial observation: @@ -92,16 +97,16 @@ def interpolate(track, delta, interpolation_type=None): newtime = np.arange(timestep[0], timestep[-1] + .01, delta) newtime[-1] = timestep[-1] - _newtime = (newtime / 24.) + time_[0] - + _newtime = (newtime / 24.) + time_[0] + date2num(np.datetime64('0000-12-31')) # Before Matplotlib 3.3, the epoch was 0000-12-31, later it changes to 1970-01-01 UTC newdates = num2date(_newtime) newdates = np.array([n.replace(tzinfo=None) for n in newdates]) if not hasattr(track, 'Speed'): idx = np.zeros(len(track.data)) idx[0] = 1 + # TODO: Possibly could change `np.mean(dt)` to `dt`? track.WindSpeed = maxWindSpeed(idx, np.mean(dt), track.Longitude, - track.Latitude, track.CentralPressure, + track.Latitude, track.CentralPressure, track.EnvPressure) # Find the indices of valid pressure observations: validIdx = np.where(track.CentralPressure < sys.maxsize)[0] @@ -139,8 +144,8 @@ def interpolate(track, delta, interpolation_type=None): if interpolation_type == 'akima': # Use the Akima interpolation method: try: - import akima - except ImportError: + from Utilities import akima + except ModuleNotFoundError: LOG.exception(("Akima interpolation module unavailable " " - default to scipy.interpolate")) nLon = splev(newtime, splrep(timestep, track.Longitude, s=0), @@ -262,6 +267,7 @@ def saveTracks(tracks, outputFile): if len(data) > 0: np.savetxt(fp, data, fmt=OUTPUT_FMTS) +@disableOnWorkers def parseTracks(configFile, trackFile, source, delta, outputFile=None, interpolation_type=None): """ @@ -299,6 +305,9 @@ def parseTracks(configFile, trackFile, source, delta, outputFile=None, if trackFile.endswith("nc"): from Utilities.track import ncReadTrackData tracks = ncReadTrackData(trackFile) + elif trackFile.endswith("xml"): + dfs = loadfile(trackFile) + tracks = [bom2tcrm(df, i) for i, df in enumerate(dfs)] else: tracks = loadTrackFile(configFile, trackFile, source) @@ -312,8 +321,51 @@ def parseTracks(configFile, trackFile, source, delta, outputFile=None, results.append(newtrack) if outputFile: - # Save data to file: - ncSaveTracks(outputFile, results) + MAX_GROUPS = 32000 + for i in range(192000, len(results), MAX_GROUPS): + part = results[i:i + MAX_GROUPS] + part_file = f"{outputFile}_part{i // MAX_GROUPS + 1:03d}.nc" + ncSaveTracks(part_file, part) + + # results = results[128000:] + # if outputFile: + # from mpi4py import MPI + # import os + + # comm = MPI.COMM_WORLD + # rank = comm.Get_rank() + # size = comm.Get_size() + # breakpoint() + # # === Parallel output file name === + # # Example: output = "tracks_part0000.nc", "tracks_part0001.nc", ... + # base, ext = os.path.splitext(outputFile) + # outputFile_rank = f"{base}_part{rank:04d}{ext}" + + # # === Each rank writes its own file === + # ncSaveTracks(outputFile_rank, results) + + # # === Optional synchronization (not strictly needed) === + # comm.Barrier() + return results + + +def bom2tcrm(df, trackId): + """ + Transforms a dataframe in BoM format into a tcrm track. + + """ + df['Datetime'] = pd.to_datetime(df.validtime) + df['Speed'] = df.translation_speed + df['CentralPressure'] = df.pcentre + df['Longitude'] = df.longitude + df['Latitude'] = df.latitude + df['EnvPressure'] = df.poci + df['rMax'] = df.rmax + df['trackId'] = df.disturbance.values + + track = Track(df) + track.trackId = [trackId, trackId] + return track diff --git a/Evaluate/windFieldValidation.py b/Evaluate/windFieldValidation.py new file mode 100644 index 00000000..4553b6eb --- /dev/null +++ b/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() diff --git a/PlotInterface/AutoPlotHazard.py b/PlotInterface/AutoPlotHazard.py index 0d2eb505..a4225f8e 100644 --- a/PlotInterface/AutoPlotHazard.py +++ b/PlotInterface/AutoPlotHazard.py @@ -36,11 +36,12 @@ from Utilities import pathLocator from Utilities import metutils +from database.queries import locationRecords + from PlotInterface.maps import saveHazardMap from PlotInterface.curves import saveHazardCurve import sqlite3 -import unicodedata log = logging.getLogger(__name__) log.addHandler(logging.NullHandler()) @@ -225,7 +226,7 @@ def plotHazardCurves(self, inputFile, plotPath): log.debug("Saving hazard curve for %s to %s"%(name, filename)) wspd = ncobj.variables['wspd'][:, j, i] - recs = database.locationRecords(self.db, pID) + recs = database.queries.locationRecords(self.db, pID) data = np.zeros(int(self.numsimulations * 365.25)) if len(recs) > 0: data[-len(recs):] = recs['wspd'] diff --git a/PlotInterface/curves.py b/PlotInterface/curves.py index b2edc6da..a00fe83a 100755 --- a/PlotInterface/curves.py +++ b/PlotInterface/curves.py @@ -149,7 +149,7 @@ def subplot(self, axes, subfigure): xdata, ydata, xlabel, ylabel, title = subfigure - axes.semilogx(xdata, ydata, '-', subsx=xdata) + axes.semilogx(xdata, ydata, '-', subs=xdata) axes.set_xlabel(xlabel) axes.set_ylabel(ylabel) axes.set_title(title) @@ -328,7 +328,7 @@ def subplot(self, axes, subfigure): """ xdata, ymean, ymax, ymin, xlabel, ylabel, title = subfigure - axes.semilogx(xdata, ymean, lw=2, subsx=xdata) + axes.semilogx(xdata, ymean, lw=2, subs=xdata) if (ymin[0] > 0) and (ymax[0] > 0): self.addRange(axes, xdata, ymin, ymax) @@ -390,7 +390,7 @@ def subplot(self, axes, subfigure): log.debug("xvalues = {0} length".format(len(emprp))) log.debug("xvalues = {0}".format(emprp)) - axes.semilogx(xdata, ymean, lw=2, subsx=xdata, + axes.semilogx(xdata, ymean, lw=2, subs=xdata, label = 'Fitted hazard curve ({0})'.format(fit)) axes.scatter(emprp[emprp > 1], events[emprp > 1], s=100, color='r', label = 'Empirical ARI') @@ -478,8 +478,8 @@ def subplot(self, axes, subfigure): """ xdata, y1, y2, y2max, y2min, xlabel, ylabel, title = subfigure - axes.semilogx(xdata, y1, color='r', lw=2, label="", subsx=xdata) - axes.semilogx(xdata, y2, color='k', lw=2, label="", subsx=xdata) + axes.semilogx(xdata, y1, color='r', lw=2, label="", subs=xdata) + axes.semilogx(xdata, y2, color='k', lw=2, label="", subs=xdata) self.addRange(axes, xdata, y2min, y2max) ylim = (0., np.max([100, np.ceil(y2.max()/10.)*10.])) axes.set_ylim(ylim) diff --git a/PlotInterface/maps.py b/PlotInterface/maps.py index 53d4b3e1..c719ff0c 100644 --- a/PlotInterface/maps.py +++ b/PlotInterface/maps.py @@ -380,7 +380,7 @@ def subplot(self, axes, subfigure): cmap = selectColormap(lvls) CS = mapobj.contourf(mx, my, data, levels=lvls, extend='both', cmap=cmap) - CB = self.colorbar(CS, ticks=lvls[::2], ax=axes, extend='both') + CB = self.colorbar(CS, ticks=lvls[::2], ax=axes) CB.set_label(cbarlab) axes.set_title(title) self.addGraticule(axes, mapobj) diff --git a/PlotInterface/plotStats.py b/PlotInterface/plotStats.py index bd2be632..42027f2d 100644 --- a/PlotInterface/plotStats.py +++ b/PlotInterface/plotStats.py @@ -15,9 +15,11 @@ import matplotlib.pyplot as plt import matplotlib.dates as mdates import matplotlib.ticker as ticker +import pandas as pd from PlotInterface.figures import LaggedRegressionFigure -from scipy.stats import linregress, probplot, frechet_l +# from scipy.stats import linregress, probplot, frechet_l +from scipy.stats import linregress, probplot, weibull_max # The distributions frechet_r and frechet_l were renamed to weibull_min and weibull_max after SciPy version 1.6.0. import numpy as np import seaborn as sns @@ -96,8 +98,9 @@ def scatterHistogram(self, xdata, ydata, labels, name, i = np.where((xdata < sys.maxsize) & (ydata < sys.maxsize))[0] xx = transform(xdata[i]) - yy = transform(ydata[i]) - jp = sns.jointplot(xx, yy, kind='reg', + yy = transform(ydata[i]) + df = pd.DataFrame({'xx': xx, 'yy': yy} ) + jp = sns.jointplot(data=df, x=df.xx, y=df.yy, kind='reg', joint_kws={'scatter_kws': {'color':'slategray', 'alpha':0.5}}, @@ -120,7 +123,8 @@ def barPlot(self, xdata, ydata, name, labels): """ fig, ax = plt.subplots(1, 1) - sns.barplot(xdata, ydata, ax=ax) + df = pd.DataFrame({'xx': xdata, 'yy': ydata} ) + sns.barplot(data=df, x=df.xx, y=df.yy, ax=ax) ax.set_xlabel(labels[0]) ax.set_ylabel(labels[1]) ax.axhline(np.mean(ydata)) @@ -180,7 +184,7 @@ def minPressureHist(self, index, pAllData): pbins = np.arange(850., 1020., 5) pcarray = np.array(pcarray) pc = np.take(pcarray, np.where(pcarray rows else processing_segment_size for x_offset in range(0, cols, processing_segment_size): segment_count = segment_count + 1 width = cols - x_offset if x_offset + processing_segment_size > cols else processing_segment_size - segments.append([x_offset, y_offset, width, height, segment_count, total_segments]) + segment_queue.put([x_offset, y_offset, width, height, segment_count, total_segments]) log.info("Lunching {0} segmented task in {1} worker threads".format(total_segments, max_working_threads)) - for seg in segments: - future_requests.append(e.submit(processMultiplierSegment, seg, source_dir_bands, wind_prj, bear_prj, dst_band)) + for _ in range(max_working_threads): + future_requests.append(e.submit(call_process_multiplier_segment, segment_queue, source_dir_bands, wind_prj, bear_prj, dst_band)) + futures.wait(future_requests, return_when='FIRST_EXCEPTION') for task in future_requests: - task.result() # Called to obtain exception information if any + task.result() # Called to obtain exception information if any + dst_ds.FlushCache() + del dst_ds - del dst_ds - print("") log.info("Completed") - return output_file +def call_process_multiplier_segment(segment_queue, source_dir_band, wind_prj, bear_prj, dst_band): + while not segment_queue.empty(): + processMultiplierSegment(segment_queue.get(), source_dir_band, wind_prj, bear_prj, dst_band) + def processMultiplierSegment(segment, source_dir_band, wind_prj, bear_prj, dst_band): """ Calculates local wind multiplier data by image segments @@ -1081,8 +1081,6 @@ def processMultiplierSegment(segment, source_dir_band, wind_prj, bear_prj, dst_b 8: {'dir': 'n', 'min': 337.5, 'max': 360.} } [x_offset, y_offset, width, height, segment_id, total_segments] = segment - log.debug("Processing segment {0}/{1}: {2} {3} {4} {5}" - .format(segment_id, total_segments, x_offset, y_offset, width, height)) with threadLock_gust: wind_data = wind_prj.ReadAsArray(x_offset, y_offset, width, height) with threadLock_bear: @@ -1096,11 +1094,9 @@ def processMultiplierSegment(segment, source_dir_band, wind_prj, bear_prj, dst_b local[idx] = wind_data[idx] * m4[idx] with threadLock_out: dst_band.WriteArray(local, x_offset, y_offset) - print('\rProgress: {0:.2f}'.format((segment_id * 100) / total_segments), "%", end="") - if segment_id % int(math.ceil(total_segments / 20)) == 0: - if log.getLogger(__name__).getEffectiveLevel() == log.DEBUG: - print("") - log.debug('Progress: {0} %'.format(int((segment_id * 100) / total_segments))) + dst_band.FlushCache() + if segment_id % int(math.ceil(total_segments / 100.0)) == 0: + log.info('Progress: {0:.2f} %'.format((segment_id * 100.0) / total_segments)) class run(): diff --git a/README.rst b/README.rst index e8bd5830..451c6962 100644 --- a/README.rst +++ b/README.rst @@ -65,8 +65,8 @@ TCRM requires: Status ====== -.. image:: https://travis-ci.org/GeoscienceAustralia/tcrm.svg?branch=develop - :target: https://travis-ci.org/GeoscienceAustralia/tcrm +.. image:: https://github.com/GeoscienceAustralia/tcrm/actions/workflows/tcrm-tests.yml/badge.svg?branch=develop + :target: https://github.com/GeoscienceAustralia/tcrm/actions/workflows/tcrm-tests.yml :alt: Build status @@ -79,8 +79,8 @@ Status :target: https://landscape.io/github/GeoscienceAustralia/tcrm/develop :alt: Code Health -.. image:: https://zenodo.org/badge/DOI/10.5281/zenodo.4070660.svg - :target: https://doi.org/10.5281/zenodo.4070660 +.. image:: https://zenodo.org/badge/10637300.svg + :target: https://zenodo.org/badge/latestdoi/10637300 Screenshot ========== diff --git a/StatInterface/SamplingOrigin.py b/StatInterface/SamplingOrigin.py index 5dea26aa..eb8e1068 100644 --- a/StatInterface/SamplingOrigin.py +++ b/StatInterface/SamplingOrigin.py @@ -18,7 +18,6 @@ from Utilities.files import flLoadFile, flSaveFile from Utilities.grid import grdRead, grdReadFromNetcdf import numpy as np -import scipy import Utilities.stats as stats LOG = logging.getLogger() @@ -119,8 +118,8 @@ def setKDEOrigins(self, kdeOriginX=None, kdeOriginY=None, kdeOriginZ=None, def generateOneSample(self): """Generate a random cyclone origin.""" # generate 2 uniform random variables - unifX = scipy.rand() - unifY = scipy.rand() + unifX = np.random.rand() + unifY = np.random.rand() xi = np.array(self.cdfX).searchsorted(unifX) yj = self.cdfY[xi, :].searchsorted(unifY) @@ -178,8 +177,8 @@ def generateSamples(self, ns, outputFile=None): raise ValueError # Generate 2 vectors of uniform random variables - unifX = scipy.rand(ns) - unifY = scipy.rand(ns) + unifX = np.random.rand(ns) + unifY = np.random.rand(ns) self.oLon = np.empty(ns, 'd') self.oLat = np.empty(ns, 'd') diff --git a/StatInterface/SamplingParameters.py b/StatInterface/SamplingParameters.py index 948578e0..ceccbfb7 100644 --- a/StatInterface/SamplingParameters.py +++ b/StatInterface/SamplingParameters.py @@ -16,7 +16,6 @@ import sys import logging -import scipy import numpy as np from Utilities.config import cnfGetIniValue from Utilities.files import flLoadFile, flSaveFile @@ -76,7 +75,7 @@ def setParameters(self, cdfParameters): def generateOneSample(self): """Generate a single random sample of cyclone parameters.""" - unif = scipy.rand() + unif = np.random.rand() ind_kdf = self.xacy[:, 1].searchsorted(unif) return self.xacy[ind_kdf, 0] @@ -96,7 +95,7 @@ def generateSamples(self, ns, sample_parameter_path=None): if ns <= 0: raise ValueError('invalid input on ns: number of sample cannot be zero or negative') - unif_s = scipy.rand(ns) + unif_s = np.random.rand(ns) ind_kdf = self.xacy[:, 1].searchsorted(unif_s) self.sample = self.xacy[ind_kdf, 0] diff --git a/StatInterface/StatInterface.py b/StatInterface/StatInterface.py index 86fe020c..802588e4 100644 --- a/StatInterface/StatInterface.py +++ b/StatInterface/StatInterface.py @@ -62,7 +62,7 @@ def __init__(self, configFile, autoCalc_gridLimit=None, gridLimitStr = cnfGetIniValue(self.configFile, 'StatInterface', 'gridLimit', '') - if gridLimitStr is not '': + if gridLimitStr != '': try: self.gridLimit = eval(gridLimitStr) except SyntaxError: diff --git a/TrackGenerator/trackSize.py b/TrackGenerator/trackSize.py index d1a92df5..f46a7c1d 100644 --- a/TrackGenerator/trackSize.py +++ b/TrackGenerator/trackSize.py @@ -18,16 +18,14 @@ LOG = logging.getLogger() -def rmax(dp, lat, eps, coeffs=[3.5843946536979779,-0.0045486143609339436, - 0.78621467400844858, 0.0024030344245284741, - 0.0015567629057007433]): +def rmax(dp, lat, eps, coeffs=[4.22, -0.0198, 0.0023]): """ Calculate radius to maximum wind based on pressure deficit and latitude. This function allows for the random variate to be set when calling the function. Default coefficients for the functional form of ln(Rmw) are given, based on JTWC data for the southern hemisphere. - ln(Rmw) = a + b*dp + c*exp(-d*dp^2) + f*|lat| + eps + ln(Rmw) = a + b*dp + c*|lat| + eps eps is not included in the coefficients (though that may be considered by some to be more logical), so that it can remain constant for a single @@ -43,18 +41,16 @@ def rmax(dp, lat, eps, coeffs=[3.5843946536979779,-0.0045486143609339436, :returns: radius to maximum wind value. """ - if len(coeffs) < 4: - LOG.warn("Insufficient coefficients for rmw calculation!") - LOG.warn("Using default values") - coeffs = [3.5843946536979779,-0.0045486143609339436, - 0.78621467400844858, 0.0024030344245284741, - 0.0015567629057007433] + if len(coeffs) != 3: + LOG.warning("Insufficient coefficients for rmw calculation!") + LOG.warning("Using default values") + + coeffs = [4.22, -0.0198, 0.0023] if isinstance(dp, (np.ndarray, list)) and \ isinstance(lat, (np.ndarray, list)): assert len(dp) == len(lat) - yy = coeffs[0] + coeffs[1]*dp + coeffs[2] * np.exp(-coeffs[3] * dp * dp) +\ - coeffs[4] * np.abs(lat) + eps + yy = coeffs[0] + coeffs[1] * dp + coeffs[2] * np.abs(lat) + eps rm = np.exp(yy) return rm @@ -65,7 +61,7 @@ def fitRmax(rmw, dp, lat): We fit a function of dp and latitude to ln(Rmw) values of the form: - ln(Rmw) = a + b*dp + c*dp^2 + d*lat^2 + eps + ln(Rmw) = a + b*dp + c*|lat| + eps where eps is a random normal variate with zero mean and std. dev. describing the residual variance. @@ -84,8 +80,7 @@ def fitRmax(rmw, dp, lat): assert len(dp) == len(lat) assert len(rmw) == len(dp) - X = np.column_stack((dp, dp*dp, lat*lat)) - X = sm.add_constant(X) + X = np.column_stack((dp, abs(lat))) y = np.array(np.log(rmw)) model = sm.OLS(y, X) results = model.fit() diff --git a/Utilities/config.py b/Utilities/config.py index facc37b5..10a65ee2 100644 --- a/Utilities/config.py +++ b/Utilities/config.py @@ -262,7 +262,7 @@ class _ConfigParser(RawConfigParser): ignoreSubsequent = True def __init__(self, defaults=DEFAULTS): RawConfigParser.__init__(self) - self.readfp(io.StringIO(defaults)) + self.read_file(io.StringIO(defaults)) self.readonce = False def geteval(self, section, option): diff --git a/Utilities/files.py b/Utilities/files.py index 26218268..a58f7f81 100644 --- a/Utilities/files.py +++ b/Utilities/files.py @@ -5,6 +5,7 @@ import datetime import numpy as np from time import ctime, localtime, strftime +from pathlib import Path import hashlib @@ -29,6 +30,7 @@ def flModulePath(level=1): """ filename = os.path.realpath(sys._getframe(level).f_code.co_filename) path, fname = os.path.split(filename) + path = str(Path(path).resolve()) path.replace(os.path.sep, '/') base, ext = os.path.splitext(fname) return path, base, ext diff --git a/Utilities/loadData.py b/Utilities/loadData.py index 1e57137e..267b776b 100644 --- a/Utilities/loadData.py +++ b/Utilities/loadData.py @@ -194,7 +194,7 @@ def getSpeedBearing(index, lon, lat, deltatime, ieast=1, speed = dist / deltatime # Delete speeds less than 0, greated than 200, # or where indicator == 1. - np.putmask(speed, (speed < 0), missingValue) + np.putmask(speed, (speed < 0), missingValue) np.putmask(speed, (speed > 200), missingValue) np.putmask(speed, index, missingValue) np.putmask(speed, np.isnan(speed), missingValue) @@ -360,7 +360,6 @@ def getInitialPositions(data): except ValueError: LOG.error("'num' field cannot be converted to an integer") - raise KeyError(('Insufficient input file columns have been specified' 'Check the input file has enough fields to determine' 'TC starting positions')) @@ -688,8 +687,8 @@ def getPoci(penv, pcentre, lat, jdays, eps, """ if len(coeffs) < 6: - LOG.warn("Insufficient coefficients for poci calculation") - LOG.warn("Using default values") + LOG.warning("Insufficient coefficients for poci calculation") + LOG.warning("Using default values") coeffs=[2324.1564738613392, -0.6539853183796136, -1.3984456535888878, 0.00074072928008818927, 0.0044469231429346088, -1.4337623534206905] @@ -710,7 +709,7 @@ def getPoci(penv, pcentre, lat, jdays, eps, nvidx = np.where(pcentre == missingValue) poci_model[nvidx] = np.nan - nvidx = np.where(penv <= pcentre) + nvidx = np.where(penv < pcentre) poci_model[nvidx] = np.nan elif penv < pcentre: @@ -788,9 +787,12 @@ def getMaxWind(track, missingValue=sys.maxsize): else: track.trackMaxWind = w[w != missingValue].max() - +# CalculateWindSpeed set to False as wind speed is output from Lin model, +# although WindSpeed will not impact the final wind field +# as wind field only calculate from pressure and wind profile +# Set to True for TC tracks from best track def loadTrackFile(configFile, trackFile, source, missingValue=0, - calculateWindSpeed=True): + calculateWindSpeed=False): """ Load TC track data from the given input file, from a specified source. The configFile is a configuration file that contains a section called diff --git a/Utilities/maputils.f90 b/Utilities/maputils.f90 new file mode 100644 index 00000000..d2ecb3ca --- /dev/null +++ b/Utilities/maputils.f90 @@ -0,0 +1,54 @@ +subroutine beardist(cLon, cLat, lonArray, latArray, bearing, dist, nlon, nlat) + !$ use omp_lib + + integer, intent(in) :: nlon, nlat + doubleprecision, intent(in) :: lonArray(nlon), latArray(nlat) + doubleprecision, intent(inout), dimension(nlat, nlon) :: bearing, dist + + doubleprecision :: toRads, cLon_, cLat_, dlon, dlat, lon(nlon), lat(nlat), radius + doubleprecision :: dLon_sin(nlon), dLon_cos(nlon), lat_sin(nlat), lat_cos(nlat) + doubleprecision :: dLat_sin(nlat), dhalfLon_sin(nlon), a, c, pi + pi = 4.d0*datan(1.d0) + toRads = 0.017453292519943295 + radius = 6367.0 + + cLon_ = cLon; + cLat_ = cLat; + + cLon_ = cLon_ * toRads + cLat_ = cLat_ * toRads + + cLat_cos = cos(cLat_) + cLat_sin = sin(cLat_) + + do i = 1, nlon + lon(i) = lonArray(i) * toRads + dLon = lon(i) - cLon_ + dLon_sin(i) = sin(dLon) + dLon_cos(i) = cos(dLon) + dhalfLon_sin(i) = sin(0.5 * dLon) + end do + + do i = 1, nlat + lat(i) = latArray(i) * toRads + lat_sin(i) = sin(lat(i)) + lat_cos(i) = cos(lat(i)) + dLat_sin(i) = sin(0.5 * (lat(i) - cLat_)) + end do + + !$OMP PARALLEL DO shared(bearing, dist) + do j = 1, nlat + do i = 1, nlon + + alpha = dLon_sin(i) * lat_cos(j); + beta = (cLat_cos * lat_sin(j)) - (cLat_sin * lat_cos(j) * dLon_cos(i)); + bearing(j, i) = 0.5 * pi - atan2(alpha, beta); + + a = dLat_sin(j) * dLat_sin(j) + cLat_cos * lat_cos(j) * dhalfLon_sin(i) * dhalfLon_sin(i); + c = 2.0 * atan2(sqrt(abs(a)), sqrt(1.0 - a)); + dist(j, i) = max(radius*c, 1e-30); + end do + end do + !$OMP END PARALLEL DO + +end subroutine beardist \ No newline at end of file diff --git a/Utilities/maputils.py b/Utilities/maputils.py index 3fff81fe..1fddf8c5 100644 --- a/Utilities/maputils.py +++ b/Utilities/maputils.py @@ -17,6 +17,11 @@ import numpy as np import math from . import metutils +import warnings +try: + from . import fmaputils +except ImportError: + warnings.warn("Compiled maputils not found - defaulting to slower python wind models") # C weave code disabled for now. The code speeds up the windfield interface module by ~6% but @@ -508,9 +513,20 @@ def makeGrid(cLon, cLat, margin=2, resolution=0.01, minLon=None, maxLon=None, xGrid = np.array(np.arange(minLon_, maxLon_, gridSize), dtype=int) yGrid = np.array(np.arange(minLat_, maxLat_, gridSize), dtype=int) - R = gridLatLonDist(cLon, cLat, xGrid / 1000., yGrid / 1000.) - np.putmask(R, R==0, 1e-30) - theta = np.pi/2. - gridLatLonBear(cLon, cLat, xGrid / 1000., yGrid / 1000.) + try: + from .fmaputils import beardist + lonArray = xGrid / 1000. + latArray = yGrid / 1000. + R = np.zeros((len(latArray), len(lonArray)), order='F') + theta = np.zeros((len(latArray), len(lonArray)), order='F') + + beardist(cLon, cLat, lonArray, latArray, theta, R) + R = np.ascontiguousarray(R) + theta = np.ascontiguousarray(theta) + except ImportError: + R = gridLatLonDist(cLon, cLat, xGrid / 1000., yGrid / 1000.) + theta = np.pi/2. - gridLatLonBear(cLon, cLat, xGrid / 1000., yGrid / 1000.) + np.putmask(R, R == 0, 1e-30) return R, theta diff --git a/Utilities/nctools.py b/Utilities/nctools.py index bb68bbc9..2be0b93f 100644 --- a/Utilities/nctools.py +++ b/Utilities/nctools.py @@ -254,7 +254,7 @@ def ncCreateVar(ncobj, name, dimensions, dtype, data=None, atts=None, **kwargs): def ncSaveGrid(filename, dimensions, variables, nodata=-9999, datatitle=None, gatts={}, writedata=True, - keepfileopen=False, zlib=True, complevel=4, lsd=None): + keepfileopen=False, zlib=True, complevel=1, lsd=None): """ Save a gridded dataset to a netCDF file using NetCDF4. diff --git a/Utilities/process.py b/Utilities/process.py index dbcb1bd9..02cc2ec3 100644 --- a/Utilities/process.py +++ b/Utilities/process.py @@ -92,7 +92,7 @@ def pGetProcessedFiles(datFileName=None): fh = open(datFileName) except IOError: - LOGGER.warn("Couldn't open dat file %s", datFileName) + LOGGER.warning("Couldn't open dat file %s", datFileName) return rc else: LOGGER.debug("Getting previously-processed files from %s", @@ -141,8 +141,8 @@ def pWriteProcessedFile(filename): fh.close() rc = 1 else: - LOGGER.warn(("Dat file name not provided. " - "Can't record %s as processed."), filename) + LOGGER.warning(("Dat file name not provided. " + "Can't record %s as processed."), filename) return rc @@ -161,7 +161,7 @@ def pDeleteDatFile(): if os.unlink(GLOBAL_DATFILE): rc = 1 else: - LOGGER.warn("Cannot remove dat file %s", GLOBAL_DATFILE) + LOGGER.warning("Cannot remove dat file %s", GLOBAL_DATFILE) return rc def pAlreadyProcessed(directory, filename, attribute, value): @@ -268,7 +268,7 @@ def pMoveFile(origin, destination): try: os.rename(origin, destination) except OSError: - LOGGER.warn("Error moving %s to %s", origin, destination) + LOGGER.warning("Error moving %s to %s", origin, destination) rc = 0 else: LOGGER.debug("%s moved to %s", origin, destination) diff --git a/Utilities/shptools.py b/Utilities/shptools.py index d18297f7..daa6054c 100644 --- a/Utilities/shptools.py +++ b/Utilities/shptools.py @@ -314,9 +314,9 @@ def shpGetField(shape_file, field_name, dtype=float): field_names = [fields[i][0] for i in range(len(fields))] if field_name not in field_names: - log.warn("No field '{0}' in the list of fieldnames" . + log.warning("No field '{0}' in the list of fieldnames" . format(field_name)) - log.warn("Unable to proceed with processing") + log.warning("Unable to proceed with processing") raise ValueError records = sf.records() diff --git a/Utilities/timeseries.py b/Utilities/timeseries.py index 087644de..8f3868de 100644 --- a/Utilities/timeseries.py +++ b/Utilities/timeseries.py @@ -171,6 +171,9 @@ def __init__(self, configFile): stnlat = stndata[:, 2].astype(float) for sid, lon, lat in zip(stnid, stnlon, stnlat): self.stations.append(Station(sid, lon, lat)) + + self.station_lat = np.array([s.lat for s in self.stations]) + self.station_lon = np.array([s.lon for s in self.stations]) log.info(f"There are {len(self.stations)} stations that will collect timeseries data") def sample(self, lon, lat, spd, uu, vv, prs, gridx, gridy): @@ -216,20 +219,33 @@ def extract(self, dt, spd, uu, vv, prs, gridx, gridy): :param gridy: :class:`numpy.ndarray` of grid latitudes. """ - stns = 0 - for stn in self.stations: - if stn.insideGrid(gridx, gridy): - stns += 1 - result = self.sample(stn.lon, stn.lat, spd, uu, vv, prs, - gridx, gridy) - ss, ux, vy, bb, pp = result - stn.data.append((str(stn.id), dt, stn.lon, stn.lat, ss, - ux, vy, bb, pp)) + # import xarray as xr - else: - stn.data.append((str(stn.id), dt, stn.lon, stn.lat, 0.0, 0.0, - 0.0, 0.0, prs[0, 0])) - log.debug("Extracted data for {0} stations".format(stns)) + # grid_lon = gridx[0, :] + # grid_lat = gridy[:, 0] + + # mask = (self.station_lat <= grid_lat[0]) & (self.station_lat >= grid_lat[-1]) + # mask &= (self.station_lon >= grid_lon[0]) & (self.station_lon <= grid_lon[-1]) + + out = np.zeros((len(self.stations), 5)) + # for i, xx in enumerate(spd, uu, vv, prs): + # yy = xr.DataArray( + # xx, + # dims=["lat", "lon"], + # coords=dict( + # lon=grid_lon, + # lat=grid_lat, + # ), + # ) + + # out[:, i][mask] = yy.interp(lat=grid_lat[mask], lon=grid_lon.mask) + + # out[:, -1][mask] = np.mod((180. / np.pi) * np.arctan2(-out[:, 1][mask], -out[:, 2][mask]), 360.) + + # for stn in self.stations: + # stn.data.append((str(stn.id), dt, stn.lon, stn.lat, out[i, 0], out[i, 1], out[i, 2], out[i, -1], out[i, 3])) + print("huh") + log.debug("Extracted data for {0} stations".format(mask.sum())) def shutdown(self): """ diff --git a/Utilities/track.py b/Utilities/track.py index bdca9e97..c23aabaf 100644 --- a/Utilities/track.py +++ b/Utilities/track.py @@ -25,6 +25,7 @@ from Utilities.maputils import bearing2theta from netCDF4 import Dataset, date2num, num2date +from cftime import num2pydate try: from exceptions import WindowsError @@ -80,15 +81,22 @@ def __init__(self, data): self.data = data self.trackId = None self.trackfile = None - if (len(data) > 0) and ('CentralPressure' in data.dtype.names): + if (len(data) > 0) and self.has_key('CentralPressure'): self.trackMinPressure = np.min(data['CentralPressure']) else: self.trackMinPressure = None - if (len(data) > 0) and ('WindSpeed' in data.dtype.names): + if (len(data) > 0) and self.has_key('WindSpeed'): self.trackMaxWind = np.max(data['WindSpeed']) else: self.trackMaxWind = None + def has_key(self, key): + + try: + return (key in self.data.dtype.names) + except AttributeError: + return (key in self.data.columns) + def __getattr__(self, key): """ Get the `key` from the `data` object. @@ -206,11 +214,14 @@ def ncReadTrackData(trackfile): units = ncobj.getncattr('time_units') calendar = ncobj.getncattr('calendar') dtt = num2date(dt[:], units, calendar) + # Convert to true python datetimes + # dtconversion = [datetime.strptime(d.strftime(), "%Y-%m-%d %H:%M:%S") for d in dtt] + dtconversion = list(dtt) newtd = np.zeros(len(dtt), dtype=track_dtype) for f in ncobj.variables.keys(): if f != 'Datetime' and f in track_dtype.names: newtd[f] = ncobj.variables[f][:] - newtd['Datetime'] = dtt + newtd['Datetime'] = dtconversion track = Track(newtd) track.trackfile = trackfile track.trackId = eval(ncobj.trackId) @@ -237,7 +248,9 @@ def ncReadTrackData(trackfile): for f in track_data.dtype.names: if f != 'Datetime' and f in track_dtype.names: newtd[f] = track_data[f] - newtd['Datetime'] = dt + # dtconversion = [datetime.strptime(d.strftime(), "%Y-%m-%d %H:%M:%S") for d in dt] + dtconversion = list(dt) + newtd['Datetime'] = dtconversion track = Track(newtd) track.trackfile = trackfile @@ -248,7 +261,7 @@ def ncReadTrackData(trackfile): tracks.append(track) else: - log.warn(TRACK_EMPTY_GROUP.format(trackfile)) + log.warning(TRACK_EMPTY_GROUP.format(trackfile)) ncobj.close() return tracks @@ -310,9 +323,9 @@ def ncSaveTracks(trackfile, tracks, dims = tdata.createDimension('time', None) times = tdata.createVariable('time', 'f8', ('time',), - zlib=True, complevel=8, shuffle=True) + zlib=True, complevel=1, shuffle=True) tvar = tdata.createVariable('track', tdtype, ('time',), - zlib=True, complevel=8, shuffle=True) + zlib=True, complevel=1, shuffle=True) t.data['Datetime'] = date2num(t.data['Datetime'], timeunits, calendar) times[:] = t.data['Datetime'] times.units = 'hours since 1900-01-01 00:00' diff --git a/convergenceTest.py b/convergenceTest.py index 30d7d87d..1b0d86bf 100644 --- a/convergenceTest.py +++ b/convergenceTest.py @@ -54,8 +54,6 @@ # Load the configuration file from the TCHA18, then open the database # and get teh list of available locations. -# In[2]: - configFile = "/home/547/cxa547/tcrmconfig/tcrm2.1.ini" config = ConfigParser() config.read(configFile) @@ -69,7 +67,6 @@ # The following step performs the calculations. First a helper # function to add nicely formatted grid lines on a logarithmic axis. -# # The second function (`plotConvergenceTest`) loads the data from the # database, then splits into two separate collections (called `d1` and # `d2`). For each of these, we then calculate empirical ARI values and @@ -88,7 +85,7 @@ def addARIGrid(axes): axes.autoscale(True, axis='x', tight=True) axes.grid(True, which='major', linestyle='-') axes.grid(True, which='minor', linestyle='--', linewidth=0.5) - + def addAEPGrid(axes): """ Add a logarithmic graticuyle to the subplot axes @@ -99,7 +96,7 @@ def addAEPGrid(axes): axes.autoscale(True, axis='y', tight=True) axes.grid(True, which='major', linestyle='-') axes.grid(True, which='minor', linestyle='--', linewidth=0.5) - + def calculateARI(data, years): emprp = empReturnPeriod(np.sort(data)) return np.sort(data)[-years:], emprp[-years:] @@ -117,7 +114,7 @@ def plotConvergenceTest(locName): locLon = locations['locLon'][locations['locId']==locId][0] locLat = locations['locLat'][locations['locId']==locId][0] - records = database.locationRecords(db, str(locId)) + records = database.queries.locationRecords(db, str(locId)) recs = records['wspd'][records['wspd'] > 0] data = np.zeros(int(NumSimulations*365.25)) data[-len(recs):] = recs @@ -133,7 +130,7 @@ def plotConvergenceTest(locName): fdelta = delta/mn fig, ax1 = plt.subplots(1, 1, figsize=figsize) - + ax1.fill_between(rr[0,:], dd[1,:], dd[0,:], alpha=0.5, label="95th percentile") ax1.plot(emprp[-10000:], data[-10000:], color='k', label="Mean ARI") ax1.set_xscale('log') @@ -213,7 +210,7 @@ def plotConvergence(ax, locName): locLon = locations['locLon'][locations['locId']==locId][0] locLat = locations['locLat'][locations['locId']==locId][0] - records = database.locationRecords(db, str(locId)) + records = database.queries.locationRecords(db, str(locId)) recs = records['wspd'][records['wspd'] > 0] data = np.zeros(int(NumSimulations*365.25)) data[-len(recs):] = recs @@ -246,5 +243,5 @@ def plotConvergence(ax, locName): axlist[7].set_xlabel('Average recurrence interval (years)') fig.tight_layout() -plt.savefig(os.path.join(plotPath, "ARI_convergence.png"), +plt.savefig(os.path.join(plotPath, "ARI_convergence.png"), bbox_inches='tight') diff --git a/database/__init__.py b/database/__init__.py index ca703441..c4eb9411 100644 --- a/database/__init__.py +++ b/database/__init__.py @@ -38,14 +38,14 @@ import logging import sqlite3 from sqlite3 import PARSE_DECLTYPES, PARSE_COLNAMES, IntegrityError -from functools import wraps +from functools import wraps, reduce import time from datetime import datetime import unicodedata import re +import atexit from shapely.geometry import Point -logging.getLogger('shapely').setLevel(logging.WARNING) from netCDF4 import Dataset import numpy as np @@ -54,12 +54,17 @@ from Utilities.track import loadTracksFromFiles from Utilities.parallel import attemptParallel, disableOnWorkers from Utilities.process import pAlreadyProcessed, pGetProcessedFiles -from functools import reduce -sqlite3.register_adapter(np.int64, lambda val: int(val)) -sqlite3.register_adapter(np.int32, lambda val: int(val)) +from .definitions import (TBLLOCATIONDEF, TBLEVENTSDEF, TBLWINDSPEEDDEF, + TBLHAZARDDEF, TBLTRACKSDEF, INSLOCATIONS, + INSEVENTS, INSWINDSPEED, INSHAZARD, INSTRACK, + SELECTLOCATIONS) + +sqlite3.register_adapter(np.int64, int) +sqlite3.register_adapter(np.int32, int) log = logging.getLogger(__name__) log.addHandler(logging.NullHandler()) +logging.getLogger('shapely').setLevel(logging.WARNING) def fromrecords(records, names): """ Convert records to array, even if no data """ @@ -90,70 +95,6 @@ def wrap(*args, **kwargs): # pylint: disable=R0914,R0902 -# Table definition statements -# Stations - we assume a geographic coordinate system: -TBLLOCATIONDEF = ("CREATE TABLE IF NOT EXISTS tblLocations " - "(locId integer PRIMARY KEY, locCode text, " - "locName text, locType text, locLon real, " - "locLat real, locElev real, locCountry text, " - "locSource text, Comments text, " - "dtCreated timestamp)") - -# Events: -TBLEVENTSDEF = ("CREATE TABLE IF NOT EXISTS tblEvents " - "(eventNumber integer PRIMARY KEY, eventId text, " - "eventFile text, eventTrackFile text, " - "eventMaxWind real, eventMinPressure real, " - "dtTrackFile timestamp, dtWindfieldFile timestamp, " - "tcrmVersion text, Comments text, dtCreated timestamp)") - -#Station wind speed from events: -TBLWINDSPEEDDEF = ("CREATE TABLE IF NOT EXISTS tblWindSpeed " - "(locId integer, eventId text, wspd real, umax real, " - "vmax real, pmin real, Comments text, " - "dtCreated timestamp)") - -# Station hazard levels: -TBLHAZARDDEF = ("CREATE TABLE IF NOT EXISTS tblHazard " - "(locId integer, returnPeriod real, wspd real, " - " wspdUpper real, wspdLower real, loc real, " - "scale real, shape real, tcrmVersion text, " - "dtHazardFile timestamp, Comments text, " - "dtCreated timestamp)") - -# Proximity of tracks to stations: -TBLTRACKSDEF = ("CREATE TABLE IF NOT EXISTS tblTracks " - "(locId integer, eventId text, distClosest real, " - "prsClosest real, dtClosest timestamp, Comments text, " - "dtCreated timestamp)") - -# Insert statements: -# Insert locations: -INSLOCATIONS = ("INSERT OR REPLACE INTO tblLocations " - "VALUES (?,?,?,?,?,?,?,?,?,?,?)") - -# Insert event record: -INSEVENTS = "INSERT INTO tblEvents VALUES (?,?,?,?,?,?,?,?,?,?,?)" - -# Insert wind speed record: -INSWINDSPEED = ("INSERT INTO tblWindSpeed " - "VALUES (?,?,?,?,?,?,?,?)") - -# Insert hazard record: -INSHAZARD = "INSERT INTO tblHazard VALUES (?,?,?,?,?,?,?,?,?,?,?,?)" - -# Insert track record: -INSTRACK = "INSERT INTO tblTracks VALUES (?,?,?,?,?,?,?)" - -# Select statements; -# Select locations within domain: -SELECTLOCATIONS = ("SELECT * FROM tblLocations WHERE " - "locLon >= ? and locLon <= ? and " - "locLat >= ? and locLat <= ?") - -# Select locId, locLon & locLat from the subset of locations: -SELECTLOCLONLAT = "SELECT locId, locLon, locLat FROM tblLocations " - def windfieldAttributes(ncobj): """ Extract the required attributes from a netCDF file. @@ -192,7 +133,7 @@ def HazardDatabase(configFile): # pylint: disable=C0103 :param str configFile: Path to configuration file """ - global _singletons + global _singletons # pylint: disable=W0603 instance = _singletons.get(configFile) if not instance: instance = _HazardDatabase(configFile) @@ -232,8 +173,6 @@ def __init__(self, configFile): detect_types=PARSE_DECLTYPES|PARSE_COLNAMES) self.exists = True - - import atexit atexit.register(self.close) @disableOnWorkers @@ -250,7 +189,6 @@ def createDatabase(self): self.createTable('tblTracks', TBLTRACKSDEF) self.exists = True self.commit() - return @disableOnWorkers def createTable(self, tblName, tblDef): @@ -317,7 +255,7 @@ def getLocations(self): locations = cur.fetchall() locations = fromrecords(locations, - names=("locId,locName,locLon,locLat")) + names=("locId,locName,locLon,locLat")) return locations def generateEventTable(self): @@ -489,25 +427,6 @@ def processEvents(self): self.insertEvents(eventparams) self.insertWindSpeeds(wsparams) - - def loadWindfieldFile(self, ncobj): - """ - Load an individual dataset. - - :param str filename: filename to load. - - :returns: tuple containing longitude, latitude, wind speed, - eastward and northward components and pressure grids. - """ - lon = ncobj.variables['lon'][:] - lat = ncobj.variables['lat'][:] - vmax = ncobj.variables['vmax'][:] - ua = ncobj.variables['ua'][:] - va = ncobj.variables['va'][:] - pmin = ncobj.variables['slp'][:] - - return (lon, lat, vmax, ua, va, pmin) - def processEvent(self, filename, locations, eventNum): """ Process an individual event file @@ -522,9 +441,13 @@ def processEvent(self, filename, locations, eventNum): log.debug("Event ID: {0}".format(eventId)) try: ncobj = Dataset(pjoin(self.windfieldPath, filename)) - except: - log.warn("Cannot open {0}".\ + except IOError as excmsg: + log.warning("Cannot open {0}".\ format(pjoin(self.windfieldPath, filename))) + log.warning(excmsg) + except: + log.exception(("Failed trying to open " + f"{pjoin(self.windfieldPath, filename)}")) # First perform the event update for tblEvents: fname = pjoin(self.windfieldPath, filename) @@ -540,7 +463,7 @@ def processEvent(self, filename, locations, eventNum): "", datetime.now()) # Perform update for tblWindSpeed: - lon, lat, vmax, ua, va, pmin = self.loadWindfieldFile(ncobj) + lon, lat, vmax, ua, va, pmin = loadWindfieldFile(ncobj) ncobj.close() wsparams = list() @@ -662,8 +585,8 @@ def processTracks(self): try: result = comm.recv(source=MPI.ANY_SOURCE, tag=MPI.ANY_TAG, status=status) - except: - log.warn("Problems recieving results on node 0") + except Exception: + log.warning("Problems recieving results on node 0") d = status.source if result: @@ -720,6 +643,24 @@ def insertTracks(self, trackRecords): log.debug("Inserted {0} records into tblTracks".format(len(trackRecords))) +def loadWindfieldFile(ncobj): + """ + Load an individual dataset. + + :param str filename: filename to load. + + :returns: tuple containing longitude, latitude, wind speed, + eastward and northward components and pressure grids. + """ + lon = ncobj.variables['lon'][:] + lat = ncobj.variables['lat'][:] + vmax = ncobj.variables['vmax'][:] + ua = ncobj.variables['ua'][:] + va = ncobj.variables['va'][:] + pmin = ncobj.variables['slp'][:] + + return (lon, lat, vmax, ua, va, pmin) + def processTrack(trackfile, locations): """ Process individual track to determine distance to locations, etc. @@ -767,7 +708,7 @@ def run(configFile): location_file = config.get('Input', 'LocationFile') buildLocationDatabase(location_db, location_file) - global MPI, comm + global MPI, comm # pylint: disable=W0601 MPI = attemptParallel() comm = MPI.COMM_WORLD db = HazardDatabase(configFile) @@ -875,219 +816,3 @@ def buildLocationDatabase(location_db, location_file, location_type='AWS'): locdb.executemany(INSLOCATIONS, locations) locdb.commit() locdb.close() - -@timer -def locationRecordsExceeding(hazard_db, locId, windSpeed): - """ - Select all records where the wind speed at the given location is - greater than some threshold. - - :param hazard_db: :class:`HazardDatabase` instance. - :param int locId: Location identifier. - :param float windSpeed: Select all records where the wind speed - at the given location is greater than - this value. - - :returns: :class:`numpy.recarray` containing the name, longitude - & latitude of the location, the wind speed of the - record, the event Id and the event file that holds the - event that generated the wind speed. - - Example:: - - >>> db = HazardDatabase(configFile) - >>> locId = 00001 - >>> records = locationRecordsExceeding(db, locId, 47.) - - """ - - query = ("SELECT l.locId, l.locName, w.wspd, w.eventId " - "FROM tblLocations l " - "INNER JOIN tblWindSpeed w ON l.locId = w.locId " - "WHERE w.wspd > ? and l.locId = ? " - "ORDER BY w.wspd ASC") - - cur = hazard_db.execute(query, (windSpeed, locId,)) - results = cur.fetchall() - results = fromrecords(results, names=('locId,locName,wspd,eventId')) - - return results - -@timer -def locationRecords(hazard_db, locId): - """ - Select all wind speed records for a given location. - - :param hazard_db: :class:`HazardDatabase` instance. - :param int locId: Location identifier. - - :returns: :class:`numpy.recarray` containing the location id, location - name, wind speed and event id. - - """ - - query = ("SELECT w.locId, l.locName, w.wspd, w.umax, w.vmax, w.eventId " - "FROM tblWindSpeed w " - "INNER JOIN tblLocations l " - "ON w.locId = l.locId " - "WHERE l.locId = ? ORDER BY w.wspd ASC") - cur = hazard_db.execute(query, (locId,)) - results = cur.fetchall() - results = fromrecords(results, - names=('locId,locName,wspd,umax,vmax,eventId')) - - return results - -@timer -def locationPassage(hazard_db, locId, distance=50): - """ - Select all records from tblTracks that pass within a defined - distance of the given location - - :param hazard_db: :class:`HazardDatabase` instance. - :param int locId: Location identifier. - :param distance: Distance threshold (in kilometres). - - :returns: :class:`numpy.recarray` containing the location id, location - name, event id, closest distance of approach, wind speed and - event file for all events that pass within the defined - distance of the selected location. - - Example:: - - >>> db = HazardDatabase(configFile) - >>> locId = 000001 - >>> records = locationPassage(db, locId, 50) - - """ - - query = ("SELECT l.locId, l.locName, t.eventId, t.distClosest, " - "w.wspd, e.eventFile FROM tblLocations l " - "INNER JOIN tblTracks t " - "ON l.locId = t.locId " - "JOIN tblWindSpeed w on w.eventId = t.eventId " - "JOIN tblEvents e on e.eventId = t.eventId " - "WHERE t.distClosest < ? and l.locId = ?") - cur = hazard_db.execute(query, (distance, locId)) - results = cur.fetchall() - results = fromrecords(results, - names=('locId,locName,eventId,' - 'distClosest,wspd,eventFile')) - return results - -@timer -def locationPassageWindSpeed(hazard_db, locId, speed, distance): - """ - Select records from _tblWindSpeed_, _tblTracks_ and _tblEvents_ that - generate a defined wind speed and pass within a given distance - of the location. - - :param hazard_db: :class:`HazardDatabase` instance. - :param int locId: Location identifier. - :param float speed: Minimum wind speed (m/s). - :param float distance: Distance threshold (kilometres). - - """ - - query = ("SELECT l.locName, w.wspd, w.umax, w.vmax, w.eventId, " - "t.distClosest, e.eventMaxWind, e.eventMinPressure " - "FROM tblLocations l " - "JOIN tblWindSpeed w on l.locId = w.locId " - "JOIN tblEvents e ON e.eventId = w.eventId " - "JOIN tblTracks t ON w.locId = t.locId AND w.eventId = t.eventId " - "WHERE l.locId = ? and w.wspd > ? AND t.distClosest <= ? " - "ORDER BY w.wspd ASC") - - cur = hazard_db.execute(query, (locId, speed, distance)) - results = cur.fetchall() - results = fromrecords(results, - names=('locName,wspd,umax,vmax,eventId,' - 'distClosest,maxwind,pmin')) - - return results - -@timer -def locationReturnPeriodEvents(hazard_db, locId, return_period): - """ - Select all records from tblEvents where the wind speed is - greater than the return period wind speed for the given return period. - - :param hazard_db: :class:`HazardDatabase` instance. - :param int locId: Location identifier. - :param int return_period: Nominated return period. - - :returns: :class:`numpy.recarray` of location id and wind speeds of - all events that are greater than the return level of the - nominated return period. - - The following example would return the wind speeds of all events that - exceed the 500-year return period wind speed for the selected location. - - Example:: - - >>> db = HazardDatabase(configFile) - >>> locId = 000001 - >>> records = locationReturnPeriodEvents(db, locId, 500) - - """ - - query = ("SELECT l.locId, h.wspd FROM tblLocations l " - "INNER JOIN tblHazard h ON l.locId = h.locId " - "WHERE h.returnPeriod = ? and l.locId = ?") - cur = hazard_db.execute(query, (return_period, locId)) - row = cur.fetchall() - return_level = row[0][1] - results = locationRecordsExceeding(hazard_db, locId, return_level) - - return results - -@timer -def locationAllReturnLevels(hazard_db, locId): - """ - Select all return level wind speeds (including upper and lower - confidence intervals) for a selected location. - - :param hazard_db: :class:`HazardDatabase` instance. - :param int locId: Location identifier. - - :returns: :class:`numpy.recarray` containing the location id, location - name, return period, return period windspeed and lower/upper - estimates of the return period wind speed. - - """ - - query = ("SELECT l.locId, l.locName, h.returnPeriod, h.wspd, " - "h.wspdLower, h.wspdUpper " - "FROM tblLocations l INNER JOIN tblHazard h " - "ON l.locId = h.locId " - "WHERE l.locId = ? " - "ORDER BY h.returnPeriod") - - cur = hazard_db.execute(query, (locId,)) - results = cur.fetchall() - results = fromrecords(results, - names=('locId,locName,returnPeriod,' - 'wspd,wspdLower,wspdUpper')) - - return results - -@timer -def selectEvents(hazard_db): - """ - Select all events from _tblEvents_. - - :param hazard_db: :class:`HazardDatabase` instance. - - :returns: :class:`numpy.recarray` containing the full listing of each - event in the table. - - """ - - query = "SELECT * FROM tblEvents ORDER BY eventMaxWind ASC" - cur = hazard_db.execute(query) - results = cur.fetchall() - names = ("eventNum,eventId,eventFile,eventTrackFile,eventMaxWind," - "eventMinPressure,dtTrackFile,dtWindfieldFile,tcrmVer," - "Comments,dtCreated") - results = fromrecords(results, names=names) - return results diff --git a/database/definitions.py b/database/definitions.py new file mode 100644 index 00000000..eb0321df --- /dev/null +++ b/database/definitions.py @@ -0,0 +1,74 @@ +""" +:mod:`definitions` -- table and statement definitions +===================================================== + +.. module:: definitions + :synopsis: Table definitions, insert statements and + query statements for the database module. +.. moduleauthor:: Craig Arthur + +""" + +# Table definition statements +# Stations - we assume a geographic coordinate system: +TBLLOCATIONDEF = ("CREATE TABLE IF NOT EXISTS tblLocations " + "(locId integer PRIMARY KEY, locCode text, " + "locName text, locType text, locLon real, " + "locLat real, locElev real, locCountry text, " + "locSource text, Comments text, " + "dtCreated timestamp)") + +# Events: +TBLEVENTSDEF = ("CREATE TABLE IF NOT EXISTS tblEvents " + "(eventNumber integer PRIMARY KEY, eventId text, " + "eventFile text, eventTrackFile text, " + "eventMaxWind real, eventMinPressure real, " + "dtTrackFile timestamp, dtWindfieldFile timestamp, " + "tcrmVersion text, Comments text, dtCreated timestamp)") + +#Station wind speed from events: +TBLWINDSPEEDDEF = ("CREATE TABLE IF NOT EXISTS tblWindSpeed " + "(locId integer, eventId text, wspd real, umax real, " + "vmax real, pmin real, Comments text, " + "dtCreated timestamp)") + +# Station hazard levels: +TBLHAZARDDEF = ("CREATE TABLE IF NOT EXISTS tblHazard " + "(locId integer, returnPeriod real, wspd real, " + " wspdUpper real, wspdLower real, loc real, " + "scale real, shape real, tcrmVersion text, " + "dtHazardFile timestamp, Comments text, " + "dtCreated timestamp)") + +# Proximity of tracks to stations: +TBLTRACKSDEF = ("CREATE TABLE IF NOT EXISTS tblTracks " + "(locId integer, eventId text, distClosest real, " + "prsClosest real, dtClosest timestamp, Comments text, " + "dtCreated timestamp)") + +# Insert statements: +# Insert locations: +INSLOCATIONS = ("INSERT OR REPLACE INTO tblLocations " + "VALUES (?,?,?,?,?,?,?,?,?,?,?)") + +# Insert event record: +INSEVENTS = "INSERT INTO tblEvents VALUES (?,?,?,?,?,?,?,?,?,?,?)" + +# Insert wind speed record: +INSWINDSPEED = ("INSERT INTO tblWindSpeed " + "VALUES (?,?,?,?,?,?,?,?)") + +# Insert hazard record: +INSHAZARD = "INSERT INTO tblHazard VALUES (?,?,?,?,?,?,?,?,?,?,?,?)" + +# Insert track record: +INSTRACK = "INSERT INTO tblTracks VALUES (?,?,?,?,?,?,?)" + +# Select statements; +# Select locations within domain: +SELECTLOCATIONS = ("SELECT * FROM tblLocations WHERE " + "locLon >= ? and locLon <= ? and " + "locLat >= ? and locLat <= ?") + +# Select locId, locLon & locLat from the subset of locations: +SELECTLOCLONLAT = "SELECT locId, locLon, locLat FROM tblLocations " diff --git a/database/queries.py b/database/queries.py new file mode 100644 index 00000000..6e93edba --- /dev/null +++ b/database/queries.py @@ -0,0 +1,250 @@ +import time +import logging as log +from functools import wraps, reduce + +import numpy as np + +def fromrecords(records, names): + """ Convert records to array, even if no data """ + # May become redundant after https://github.com/numpy/numpy/issues/1862 + if records: + rval = np.rec.fromrecords(records, names=names) + else: + rval = np.array([], [(name, 'O') for name in names.split(',')]) + + return rval + +def timer(func): + """ + A simple timing decorator for the entire process. + + """ + + @wraps(func) + def wrap(*args, **kwargs): + t1 = time.time() + res = func(*args, **kwargs) + tottime = time.time() - t1 + msg = "%02d:%02d:%02d " % \ + reduce(lambda ll, b: divmod(ll[0], b) + ll[1:], + [(tottime,), 60, 60]) + log.debug("Time for {0}: {1}".format(func.__name__, msg)) + return res + + return wrap + +@timer +def locationRecordsExceeding(hazard_db, locId, windSpeed): + """ + Select all records where the wind speed at the given location is + greater than some threshold. + + :param hazard_db: :class:`HazardDatabase` instance. + :param int locId: Location identifier. + :param float windSpeed: Select all records where the wind speed + at the given location is greater than + this value. + + :returns: :class:`numpy.recarray` containing the name, longitude + & latitude of the location, the wind speed of the + record, the event Id and the event file that holds the + event that generated the wind speed. + + Example:: + + >>> db = HazardDatabase(configFile) + >>> locId = 00001 + >>> records = locationRecordsExceeding(db, locId, 47.) + + """ + + query = ("SELECT l.locId, l.locName, w.wspd, w.eventId " + "FROM tblLocations l " + "INNER JOIN tblWindSpeed w ON l.locId = w.locId " + "WHERE w.wspd > ? and l.locId = ? " + "ORDER BY w.wspd ASC") + + cur = hazard_db.execute(query, (windSpeed, locId,)) + results = cur.fetchall() + results = fromrecords(results, names=('locId,locName,wspd,eventId')) + + return results + +@timer +def locationRecords(hazard_db, locId): + """ + Select all wind speed records for a given location. + + :param hazard_db: :class:`HazardDatabase` instance. + :param int locId: Location identifier. + + :returns: :class:`numpy.recarray` containing the location id, location + name, wind speed and event id. + + """ + + query = ("SELECT w.locId, l.locName, w.wspd, w.umax, w.vmax, w.eventId " + "FROM tblWindSpeed w " + "INNER JOIN tblLocations l " + "ON w.locId = l.locId " + "WHERE l.locId = ? ORDER BY w.wspd ASC") + cur = hazard_db.execute(query, (locId,)) + results = cur.fetchall() + results = fromrecords(results, + names=('locId,locName,wspd,umax,vmax,eventId')) + + return results + +@timer +def locationPassage(hazard_db, locId, distance=50): + """ + Select all records from tblTracks that pass within a defined + distance of the given location + + :param hazard_db: :class:`HazardDatabase` instance. + :param int locId: Location identifier. + :param distance: Distance threshold (in kilometres). + + :returns: :class:`numpy.recarray` containing the location id, location + name, event id, closest distance of approach, wind speed and + event file for all events that pass within the defined + distance of the selected location. + + Example:: + + >>> db = HazardDatabase(configFile) + >>> locId = 000001 + >>> records = locationPassage(db, locId, 50) + + """ + + query = ("SELECT l.locId, l.locName, t.eventId, t.distClosest, " + "w.wspd, e.eventFile FROM tblLocations l " + "INNER JOIN tblTracks t " + "ON l.locId = t.locId " + "JOIN tblWindSpeed w on w.eventId = t.eventId " + "JOIN tblEvents e on e.eventId = t.eventId " + "WHERE t.distClosest < ? and l.locId = ?") + cur = hazard_db.execute(query, (distance, locId)) + results = cur.fetchall() + results = fromrecords(results, + names=('locId,locName,eventId,' + 'distClosest,wspd,eventFile')) + return results + +@timer +def locationPassageWindSpeed(hazard_db, locId, speed, distance): + """ + Select records from _tblWindSpeed_, _tblTracks_ and _tblEvents_ that + generate a defined wind speed and pass within a given distance + of the location. + + :param hazard_db: :class:`HazardDatabase` instance. + :param int locId: Location identifier. + :param float speed: Minimum wind speed (m/s). + :param float distance: Distance threshold (kilometres). + + """ + + query = ("SELECT l.locName, w.wspd, w.umax, w.vmax, w.eventId, " + "t.distClosest, e.eventMaxWind, e.eventMinPressure " + "FROM tblLocations l " + "JOIN tblWindSpeed w on l.locId = w.locId " + "JOIN tblEvents e ON e.eventId = w.eventId " + "JOIN tblTracks t ON w.locId = t.locId AND w.eventId = t.eventId " + "WHERE l.locId = ? and w.wspd > ? AND t.distClosest <= ? " + "ORDER BY w.wspd ASC") + + cur = hazard_db.execute(query, (locId, speed, distance)) + results = cur.fetchall() + results = fromrecords(results, + names=('locName,wspd,umax,vmax,eventId,' + 'distClosest,maxwind,pmin')) + + return results + +@timer +def locationReturnPeriodEvents(hazard_db, locId, return_period): + """ + Select all records from tblEvents where the wind speed is + greater than the return period wind speed for the given return period. + + :param hazard_db: :class:`HazardDatabase` instance. + :param int locId: Location identifier. + :param int return_period: Nominated return period. + + :returns: :class:`numpy.recarray` of location id and wind speeds of + all events that are greater than the return level of the + nominated return period. + + The following example would return the wind speeds of all events that + exceed the 500-year return period wind speed for the selected location. + + Example:: + + >>> db = HazardDatabase(configFile) + >>> locId = 000001 + >>> records = locationReturnPeriodEvents(db, locId, 500) + + """ + + query = ("SELECT l.locId, h.wspd FROM tblLocations l " + "INNER JOIN tblHazard h ON l.locId = h.locId " + "WHERE h.returnPeriod = ? and l.locId = ?") + cur = hazard_db.execute(query, (return_period, locId)) + row = cur.fetchall() + return_level = row[0][1] + results = locationRecordsExceeding(hazard_db, locId, return_level) + + return results + +@timer +def locationAllReturnLevels(hazard_db, locId): + """ + Select all return level wind speeds (including upper and lower + confidence intervals) for a selected location. + + :param hazard_db: :class:`HazardDatabase` instance. + :param int locId: Location identifier. + + :returns: :class:`numpy.recarray` containing the location id, location + name, return period, return period windspeed and lower/upper + estimates of the return period wind speed. + + """ + + query = ("SELECT l.locId, l.locName, h.returnPeriod, h.wspd, " + "h.wspdLower, h.wspdUpper " + "FROM tblLocations l INNER JOIN tblHazard h " + "ON l.locId = h.locId " + "WHERE l.locId = ? " + "ORDER BY h.returnPeriod") + + cur = hazard_db.execute(query, (locId,)) + results = cur.fetchall() + results = fromrecords(results, + names=('locId,locName,returnPeriod,' + 'wspd,wspdLower,wspdUpper')) + + return results + +@timer +def selectEvents(hazard_db): + """ + Select all events from _tblEvents_. + + :param hazard_db: :class:`HazardDatabase` instance. + + :returns: :class:`numpy.recarray` containing the full listing of each + event in the table. + + """ + + query = "SELECT * FROM tblEvents ORDER BY eventMaxWind ASC" + cur = hazard_db.execute(query) + results = cur.fetchall() + names = ("eventNum,eventId,eventFile,eventTrackFile,eventMaxWind," + "eventMinPressure,dtTrackFile,dtWindfieldFile,tcrmVer," + "Comments,dtCreated") + results = fromrecords(results, names=names) + return results diff --git a/example/port_hedland.ini b/example/port_hedland.ini index 77f8beda..acfff15e 100644 --- a/example/port_hedland.ini +++ b/example/port_hedland.ini @@ -1,107 +1,107 @@ -[Actions] -; TCRM modules to execute -DataProcess=True -ExecuteStat=True -ExecuteTrackGenerator=True -ExecuteWindfield=True -ExecuteHazard=True -CreateDatabase=True -PlotHazard=True - -PlotData=False - -ExecuteEvaluate=False -DownloadData=False - -[DataProcess] -InputFile=Allstorms.ibtracs_wmo.v03r10.csv -Source=IBTRACS -StartSeason=1981 -FilterSeasons=True - -[Region] -; Domain for windfield and hazard calculation -gridLimit={'xMin':110.0,'xMax':125.0,'yMin':-26.0,'yMax':-15.0} -gridSpace={'x':1.0,'y':1.0} -gridInc={'x':1.0,'y':0.5} -LocalityID=250913860 -LocalityName=Port Hedland, Western Australia, Australia. - -[StatInterface] -kdeType=gau -kde2DType=Gaussian -kdeStep=0.2 - -[TrackGenerator] -NumSimulations=1000 -YearsPerSimulation=1 -SeasonSeed=403943 -TrackSeed=89333 - -[WindfieldInterface] -;TrackPath=./output/port_hedland/tracks -Margin=2.0 -Resolution=0.05 -Source=TCRM -profileType=powell -windFieldType=kepert - -[Hazard] -; Years to calculate return period wind speeds -Years=5,10,20,25,50,100,200,250,500,1000,2000,2500 -MinimumRecords=10 -CalculateCI=False -PercentileRange=90 -ExtremeValueDistribution=GPD -SmoothPlots=False - -[Input] -LocationFile = input/stationlist.shp -landmask = input/landmask.nc -mslpfile = MSLP/slp.day.ltm.nc -datasets = IBTRACS,LTMSLP -MSLPGrid=1,2,3,4,12 - -[Output] -Path=./output/port_hedland - -[Logging] -LogFile=./output/port_hedland/log/port_hedland.log -LogLevel=INFO -Verbose=False - -[Process] -ExcludePastProcessed=True -DatFile=./output/port_hedland/process/dat/port_hedland.dat - -[RMW] -GetRMWDistFromInputData=False -mean=50.0 -sigma=0.6 - -[TCRM] -; Output track files settings -Columns=index,age,lon,lat,speed,bearing,pressure,penv,rmax -FieldDelimiter=, -NumberOfHeadingLines=1 -SpeedUnits=kph -PressureUnits=hPa - -[IBTRACS] -; Input data file settings -url = ftp://eclipse.ncdc.noaa.gov/pub/ibtracs/v03r10/wmo/csv/Allstorms.ibtracs_wmo.v03r10.csv.gz -path = input -filename = Allstorms.ibtracs_wmo.v03r10.csv -columns = tcserialno,season,num,skip,skip,skip,date,skip,lat,lon,skip,pressure -fielddelimiter = , -numberofheadinglines = 3 -pressureunits = hPa -lengthunits = km -dateformat = %Y-%m-%d %H:%M:%S -speedunits = kph - -[LTMSLP] -; MSLP climatology file settings -URL = ftp://ftp.cdc.noaa.gov/Datasets/ncep.reanalysis.derived/surface/slp.day.1981-2010.ltm.nc -path = MSLP -filename = slp.day.ltm.nc +[Actions] +; TCRM modules to execute +DataProcess=True +ExecuteStat=True +ExecuteTrackGenerator=True +ExecuteWindfield=True +ExecuteHazard=False +CreateDatabase=False +PlotHazard=False + +PlotData=False + +ExecuteEvaluate=False +DownloadData=False + +[DataProcess] +InputFile=Allstorms.ibtracs_wmo.v03r10.csv +Source=IBTRACS +StartSeason=1981 +FilterSeasons=True + +[Region] +; Domain for windfield and hazard calculation +gridLimit={'xMin':110.0,'xMax':125.0,'yMin':-26.0,'yMax':-15.0} +gridSpace={'x':1.0,'y':1.0} +gridInc={'x':1.0,'y':0.5} +LocalityID=250913860 +LocalityName=Port Hedland, Western Australia, Australia. + +[StatInterface] +kdeType=gau +kde2DType=Gaussian +kdeStep=0.2 + +[TrackGenerator] +NumSimulations=1000 +YearsPerSimulation=1 +SeasonSeed=403943 +TrackSeed=89333 + +[WindfieldInterface] +;TrackPath=./output/port_hedland/tracks +Margin=2.0 +Resolution=0.05 +Source=TCRM +profileType=powell +windFieldType=kepert + +[Hazard] +; Years to calculate return period wind speeds +Years=5,10,20,25,50,100,200,250,500,1000,2000,2500 +MinimumRecords=10 +CalculateCI=False +PercentileRange=90 +ExtremeValueDistribution=GPD +SmoothPlots=False + +[Input] +LocationFile = input/stationlist.shp +landmask = input/landmask.nc +mslpfile = MSLP/slp.day.ltm.nc +datasets = IBTRACS,LTMSLP +MSLPGrid=1,2,3,4,12 + +[Output] +Path=./output/port_hedland + +[Logging] +LogFile=./output/port_hedland/log/port_hedland.log +LogLevel=INFO +Verbose=False + +[Process] +ExcludePastProcessed=True +DatFile=./output/port_hedland/process/dat/port_hedland.dat + +[RMW] +GetRMWDistFromInputData=False +mean=50.0 +sigma=0.6 + +[TCRM] +; Output track files settings +Columns=index,age,lon,lat,speed,bearing,pressure,penv,rmax +FieldDelimiter=, +NumberOfHeadingLines=1 +SpeedUnits=kph +PressureUnits=hPa + +[IBTRACS] +; Input data file settings +url = ftp://eclipse.ncdc.noaa.gov/pub/ibtracs/v03r10/wmo/csv/Allstorms.ibtracs_wmo.v03r10.csv.gz +path = input +filename = Allstorms.ibtracs_wmo.v03r10.csv +columns = tcserialno,season,num,skip,skip,skip,date,skip,lat,lon,skip,pressure +fielddelimiter = , +numberofheadinglines = 3 +pressureunits = hPa +lengthunits = km +dateformat = %Y-%m-%d %H:%M:%S +speedunits = kph + +[LTMSLP] +; MSLP climatology file settings +URL = ftp://ftp.cdc.noaa.gov/Datasets/ncep.reanalysis.derived/surface/slp.day.1981-2010.ltm.nc +path = MSLP +filename = slp.day.ltm.nc \ No newline at end of file diff --git a/hazard/__init__.py b/hazard/__init__.py index bb9670f6..f9fb04a3 100644 --- a/hazard/__init__.py +++ b/hazard/__init__.py @@ -825,7 +825,7 @@ def calculateCI(Vr, years, nodata, minRecords, yrsPerSim=1, upper = 100. - lower # 95th percentile default nrecords = Vr.shape[0] # number of years (since we have aggregated into 1/yr) - nsamples = nrecords / sample_size # number of iterations to perform + nsamples = int(nrecords / sample_size) # number of iterations to perform # RpUpper/RpLower = years x lat x lon RpUpper = nodata*np.ones((len(years), Vr.shape[1], Vr.shape[2]), dtype='f') diff --git a/hazard/evd.py b/hazard/evd.py index cf08981d..a9fefaaa 100644 --- a/hazard/evd.py +++ b/hazard/evd.py @@ -295,9 +295,8 @@ def gevfit(data, intervals, nodata=-9999., minrecords=50, yrspersim=1): # not all equal, and where there are 50 or more valid (>0) values. if data[ii].min() != data[ii].max(): if len(ii) >= minrecords: - l1, l2, l3 = lmom.samlmu(data, 3) # find 3 l-moments - # t3 = L-skewness (Hosking 1990) - t3 = l3 / l2 + l1, l2, t3 = lmom.samlmu(data, 3) # find the first 2 L-moments and the L-skewness + if (l2 <= 0.) or (np.abs(t3) >= 1.): # Reject points where the second l-moment is negative # or the ratio of the third to second is > 1, i.e. positive diff --git a/input/stationlist.dbf b/input/stationlist.dbf index 49cf4b44..7ce421c9 100644 Binary files a/input/stationlist.dbf and b/input/stationlist.dbf differ diff --git a/input/stationlist.shp b/input/stationlist.shp index e99b17da..314bb5ec 100644 Binary files a/input/stationlist.shp and b/input/stationlist.shp differ diff --git a/input/stationlist.shx b/input/stationlist.shx index f577b4f3..39bccfa3 100644 Binary files a/input/stationlist.shx and b/input/stationlist.shx differ diff --git a/postinstall.sh b/postinstall.sh deleted file mode 100644 index 3bd82652..00000000 --- a/postinstall.sh +++ /dev/null @@ -1,46 +0,0 @@ -# begin installing miniconda -if [[ "$TRAVIS_OS_NAME" != "windows" ]]; then - echo "installing miniconda for posix"; - bash $HOME/download/miniconda.sh -b -u -p $MINICONDA_PATH; -elif [[ "$TRAVIS_OS_NAME" == "windows" ]]; then - echo "folder $MINICONDA_SUB_PATH does not exist" - echo "installing miniconda for windows"; - choco install miniconda3 --params="'/JustMe /AddToPath:1 /D:$MINICONDA_PATH_WIN'"; -fi; -# end installing miniconda - -export PATH="$MINICONDA_PATH:$MINICONDA_SUB_PATH:$MINICONDA_LIB_BIN_PATH:$PATH"; - -# begin checking miniconda existance -echo "checking if folder $MINICONDA_SUB_PATH exists" -if [[ -d $MINICONDA_SUB_PATH ]]; then - echo "folder $MINICONDA_SUB_PATH exists" -else - echo "folder $MINICONDA_SUB_PATH does not exist" -fi; -# end checking miniconda existance - -source $MINICONDA_PATH/etc/profile.d/conda.sh; -hash -r; -echo $TRAVIS_OS_NAME -echo $PYTHON_VERSION -python --version - -if [[ "$TRAVIS_OS_NAME" == "windows" ]]; then - echo "Removing mpi4py from environment for windows build" - echo "Package not available in conda channels" - sed -i '/mpi4py/d' ./tcrmenv.yml -fi - -conda config --set always_yes yes --set changeps1 no; -conda update -q conda; -conda config --add channels conda-forge; -conda config --set channel_priority strict; -# Useful for debugging any issues with conda -conda info -a - -echo "Create TCRM environment" -conda env create -q -f tcrmenv.yml python=$PYTHON_VERSION; -conda activate tcrm -python --version -conda list diff --git a/preinstall.sh b/preinstall.sh deleted file mode 100644 index 339b2536..00000000 --- a/preinstall.sh +++ /dev/null @@ -1,24 +0,0 @@ -if [[ "$TRAVIS_OS_NAME" != "windows" ]]; then - export MINICONDA_PATH=$HOME/miniconda; - export MINICONDA_SUB_PATH=$MINICONDA_PATH/bin; -elif [[ "$TRAVIS_OS_NAME" == "windows" ]]; then - export MINICONDA_PATH=$HOME/miniconda; - export MINICONDA_PATH_WIN=`cygpath --windows $MINICONDA_PATH`; - export MINICONDA_SUB_PATH=$MINICONDA_PATH/Scripts; -fi; -export MINICONDA_LIB_BIN_PATH=$MINICONDA_PATH/Library/bin; - # Obtain miniconda installer -if [[ "$TRAVIS_OS_NAME" != "windows" ]]; then - mkdir -p $HOME/download; - if [[ "$TRAVIS_OS_NAME" == "linux" ]]; then - echo "downloading miniconda.sh for linux"; - wget https://repo.continuum.io/miniconda/Miniconda3-latest-Linux-x86_64.sh -O $HOME/download/miniconda.sh; - elif [[ "$TRAVIS_OS_NAME" == "osx" ]]; then - echo "downloading miniconda.sh for osx"; - wget https://repo.continuum.io/miniconda/Miniconda3-latest-MacOSX-x86_64.sh -O $HOME/download/miniconda.sh; - fi; -fi; - # Install openssl for Windows -if [[ "$TRAVIS_OS_NAME" == "windows" ]]; then - choco install openssl.light; -fi; \ No newline at end of file diff --git a/pylintrc b/pylintrc index abad3386..06de6f8f 100755 --- a/pylintrc +++ b/pylintrc @@ -1,5 +1,5 @@ [MESSAGES CONTROL] -disable=W0511,W0142,W1202,I0011,E1003,E1101,E0611,F0401,E1103,E1121 +disable=W0511,W0142,W1201,W1202,I0011,E1003,E1101,E0611,F0401,E1103,E1121,logging-fstring-interpolation [BASIC] attr-rgx=[a-zA-Z_][a-zA-Z0-9_]{0,30}[_]{0,1}$ diff --git a/run_tcrm.sh b/run_tcrm.sh new file mode 100755 index 00000000..aa5c4510 --- /dev/null +++ b/run_tcrm.sh @@ -0,0 +1,24 @@ +#!/bin/bash +#PBS -P w85 +#PBS -q normal +#PBS -N tcrm_WP_ssp126_18_reform +#PBS -m ae +#PBS -M liang.hu@ga.gov.au +#PBS -lwalltime=48:00:00 +#PBS -lmem=190GB,ncpus=48,jobfs=400GB +#PBS -joe +#PBS -lstorage=gdata/w85+scratch/w85+gdata/hh5 +#PBS -l wd + + +module purge +module load pbs +module load dot +module load openmpi + +module use /g/data/hh5/public/modules +module load conda/analysis3 +source /scratch/w85/lh9573/EXTRA_PYTHON_LIBS/tcrm/bin/activate + +# mpirun -np 48 python3 tcrm.py -c example/WP_ssp126_18.ini +mpirun -np 24 python3 tc_track_interpolation.py -c example/WP_ssp126_10000_reform.ini \ No newline at end of file diff --git a/setup.py b/setup.py index ba3b06e5..51b0b04c 100644 --- a/setup.py +++ b/setup.py @@ -7,8 +7,8 @@ setup( name = "TCRM", - version = '3.1.4', - packages=find_packages(), + version = '3.1.12', + packages=find_packages(), scripts=['tcrm.py', 'tcevent.py'], include_package_data=True, package_data = { @@ -46,7 +46,9 @@ 'imageio', 'mpi4py', 'boto3', - 'botocore'], + 'botocore', + 'pytest', + 'pytest-cov'], # metadata: author = "Craig Arthur", diff --git a/tc_track_interpolation.py b/tc_track_interpolation.py new file mode 100755 index 00000000..b17caef6 --- /dev/null +++ b/tc_track_interpolation.py @@ -0,0 +1,282 @@ +""" +:mod:`tcevent_all` -- run the windfield module for a single TC track +================================================================ + +.. module:: tcevent + :synopsis: Run the wind field module for a single TC track. + +Run the :mod:`wind.windmodels` module to calculate the wind field for a +single TC event. The track of the TC is interpolated to a fine +temporal resolution, then the maximum wind field evaluated. + +Data at selected points within the model domain can be extracted +at each time step, giving a time history of wind speed, direction +and estimated sea level pressure for the location(s). + +See the :ref:`Scenario modelling ` section of +the TCRM User Guide for details on running this script. + +""" + +import logging as log +log.getLogger('matplotlib').setLevel(log.WARNING) +from functools import reduce + +import os +import time +import argparse +import traceback + +from functools import wraps +from os.path import join as pjoin, realpath, isdir, dirname + +from Utilities import pathLocator +from Utilities.config import ConfigParser +from Utilities.files import flStartLog +from Utilities.version import version +from Utilities.parallel import attemptParallel, disableOnWorkers + +from Utilities.progressbar import SimpleProgressBar as ProgressBar +from Evaluate import interpolateTracks + +__version__ = version() + +def timer(f): + """ + Basic timing functions for entire process + """ + @wraps(f) + def wrap(*args, **kwargs): + t1 = time.time() + res = f(*args, **kwargs) + + tottime = time.time() - t1 + msg = "%02d:%02d:%02d " % \ + reduce(lambda ll, b : divmod(ll[0], b) + ll[1:], + [(tottime,), 60, 60]) + + log.info("Time for {0}: {1}".format(f.__name__, msg) ) + return res + + return wrap + +@disableOnWorkers +def doOutputDirectoryCreation(configFile): + """ + Create all the necessary output folders. + + :param str configFile: Name of configuration file. + :raises OSError: If the directory tree cannot be created. + + """ + + config = ConfigParser() + config.read(configFile) + + outputPath = config.get('Output', 'Path') + + log.info('Output will be stored under %s', outputPath) + + subdirs = ['tracks', 'windfield', 'plots', 'plots/timeseries', + 'log', 'process', 'process/timeseries'] + + if not isdir(outputPath): + try: + os.makedirs(outputPath) + except OSError: + raise + for subdir in subdirs: + if not isdir(realpath(pjoin(outputPath, subdir))): + try: + os.makedirs(realpath(pjoin(outputPath, subdir))) + except OSError: + raise + +@disableOnWorkers +def doTimeseriesPlotting(configFile): + """ + Run functions to plot time series output. + + :param str configFile: Path to configuration file. + + """ + + config = ConfigParser() + config.read(configFile) + + outputPath = config.get('Output', 'Path') + timeseriesPath = pjoin(outputPath, 'process', 'timeseries') + plotPath = pjoin(outputPath, 'plots', 'timeseries') + log.info("Plotting time series data to {0}".format(plotPath)) + from PlotInterface.plotTimeseries import plotTimeseries + plotTimeseries(timeseriesPath, plotPath) + +@disableOnWorkers +def doWindfieldPlotting(configFile): + """ + Plot the wind field on a map. + + :param str configFile: Path to the configuration file. + + :Note: the file name is assumed to be 'gust.interp.nc' + + """ + from netCDF4 import Dataset + import numpy as np + from PlotInterface.maps import saveWindfieldMap + + config = ConfigParser() + config.read(configFile) + outputPath = config.get('Output', 'Path') + windfieldPath = pjoin(outputPath, 'windfield') + + inputFile = config.get('DataProcess', 'InputFile') + if inputFile.endswith(".nc"): + # We have a netcdf track file. Work under the assumption it is + # drawn directly from TCRM. + trackFile = os.path.basename(inputFile) + trackId = trackFile.split('.')[1] + gustFile = 'gust.{0}.nc'.format(trackId) + outputWindFile = pjoin(windfieldPath, gustFile) + else: + # Note the assumption about the file name! + outputWindFile = pjoin(windfieldPath, 'gust.001-00001.nc') + plotPath = pjoin(outputPath, 'plots', 'maxwind.png') + + f = Dataset(outputWindFile, 'r') + + xdata = f.variables['lon'][:] + ydata = f.variables['lat'][:] + + vdata = f.variables['vmax'][:] + + gridLimit = None + if config.has_option('Region','gridLimit'): + gridLimit = config.geteval('Region', 'gridLimit') + ii = np.where((xdata >= gridLimit['xMin']) & + (xdata <= gridLimit['xMax'])) + jj = np.where((ydata >= gridLimit['yMin']) & + (ydata <= gridLimit['yMax'])) + [xgrid, ygrid] = np.meshgrid(xdata[ii], ydata[jj]) + ig, jg = np.meshgrid(ii, jj) + vdata = vdata[jg, ig] + else: + [xgrid, ygrid] = np.meshgrid(xdata, ydata) + map_kwargs = dict(llcrnrlon=xgrid.min(), + llcrnrlat=ygrid.min(), + urcrnrlon=xgrid.max(), + urcrnrlat=ygrid.max(), + projection='merc', + resolution='i') + title = "Maximum wind speed" + cbarlabel = "Wind speed ({0})".format(f.variables['vmax'].units) + levels = np.arange(30, 101., 5.) + saveWindfieldMap(vdata, xgrid, ygrid, title, levels, + cbarlabel, map_kwargs, plotPath) + +@timer +def main(configFile): + """ + Main function to execute the :mod:`wind`. + + :param str configFile: Path to configuration file. + + """ + config = ConfigParser() + config.read(configFile) + doOutputDirectoryCreation(configFile) + + trackFile = config.get('DataProcess', 'InputFile') + source = config.get('DataProcess', 'Source') + delta = 1. # time resolution, 1 hour + outputPath = pjoin(config.get('Output','Path'), 'tracks') + outputTrackFile = pjoin(outputPath, "tracks.interp.nc") + + # This will save interpolated track data in TCRM format: + interpTrack = interpolateTracks.parseTracks(configFile, trackFile, + source, delta, + outputTrackFile, + interpolation_type='akima') + comm.barrier() + + + +def startup(): + """ + Parse the command line arguments and call the :func:`main` + function. + + """ + print("starting up") + parser = argparse.ArgumentParser() + parser.add_argument('-c', '--config_file', + help='Path to configuration file') + parser.add_argument('-v', '--verbose', help='Verbose output', + action='store_true') + parser.add_argument('-d', '--debug', help='Allow pdb traces', + action='store_true') + args = parser.parse_args() + + configFile = args.config_file + config = ConfigParser() + config.read(configFile) + + rootdir = pathLocator.getRootDirectory() + os.chdir(rootdir) + + logfile = config.get('Logging','LogFile') + logdir = dirname(realpath(logfile)) + + # If log file directory does not exist, create it + if not isdir(logdir): + try: + os.makedirs(logdir) + except OSError: + logfile = pjoin(os.getcwd(), 'tcrm.log') + + logLevel = config.get('Logging', 'LogLevel') + verbose = config.getboolean('Logging', 'Verbose') + datestamp = config.getboolean('Logging', 'Datestamp') + debug = False + + if args.verbose: + verbose = True + + if args.debug: + debug = True + + global MPI, comm + MPI = attemptParallel() + import atexit + atexit.register(MPI.Finalize) + comm = MPI.COMM_WORLD + + flStartLog(logfile, logLevel, verbose, datestamp) + # Switch off minor warning messages + import warnings + warnings.filterwarnings("ignore", category=DeprecationWarning) + warnings.filterwarnings("ignore", category=UserWarning, module="pytz") + warnings.filterwarnings("ignore", category=UserWarning, module="numpy") + warnings.filterwarnings("ignore", category=UserWarning, + module="matplotlib") + + warnings.filterwarnings("ignore", category=RuntimeWarning) + + if debug: + main(configFile) + else: + try: + main(configFile) + except ImportError as e: + log.critical("Missing module: {0}".format(e)) + except Exception: # pylint: disable=W0703 + # Catch any exceptions that occur and log them (nicely): + tblines = traceback.format_exc().splitlines() + for line in tblines: + log.critical(line.lstrip()) + + +if __name__ == "__main__": + startup() + + diff --git a/tcevent.py b/tcevent.py index f48962d3..a50c6503 100755 --- a/tcevent.py +++ b/tcevent.py @@ -34,6 +34,8 @@ from Utilities.config import ConfigParser from Utilities.files import flStartLog from Utilities.version import version +from Utilities.parallel import attemptParallel, disableOnWorkers + from Utilities.progressbar import SimpleProgressBar as ProgressBar from Evaluate import interpolateTracks @@ -58,6 +60,7 @@ def wrap(*args, **kwargs): return wrap +@disableOnWorkers def doOutputDirectoryCreation(configFile): """ Create all the necessary output folders. @@ -89,6 +92,7 @@ def doOutputDirectoryCreation(configFile): except OSError: raise +@disableOnWorkers def doTimeseriesPlotting(configFile): """ Run functions to plot time series output. @@ -107,6 +111,7 @@ def doTimeseriesPlotting(configFile): from PlotInterface.plotTimeseries import plotTimeseries plotTimeseries(timeseriesPath, plotPath) +@disableOnWorkers def doWindfieldPlotting(configFile): """ Plot the wind field on a map. @@ -192,6 +197,7 @@ def main(configFile): source, delta, outputTrackFile, interpolation_type='akima') + comm.barrier() showProgressBar = config.get('Logging', 'ProgressBar') @@ -218,6 +224,7 @@ def startup(): function. """ + print("starting up") parser = argparse.ArgumentParser() parser.add_argument('-c', '--config_file', help='Path to configuration file') @@ -255,6 +262,12 @@ def startup(): if args.debug: debug = True + global MPI, comm + MPI = attemptParallel() + import atexit + atexit.register(MPI.Finalize) + comm = MPI.COMM_WORLD + flStartLog(logfile, logLevel, verbose, datestamp) # Switch off minor warning messages import warnings diff --git a/tcrmenv.yml b/tcrmenv.yml index 84d02fa7..f5e7d165 100644 --- a/tcrmenv.yml +++ b/tcrmenv.yml @@ -23,7 +23,7 @@ dependencies: - libgdal - gdal - configparser - - cartopy + - cartopy>=0.18.0 - affine - tqdm - xarray @@ -32,3 +32,5 @@ dependencies: - mpi4py - boto3 - botocore + - pytest + - pytest-cov diff --git a/tests/NumpyTestCase.py b/tests/NumpyTestCase.py index f07c54fa..e1485c4e 100644 --- a/tests/NumpyTestCase.py +++ b/tests/NumpyTestCase.py @@ -9,7 +9,8 @@ $Id: NumpyTestCase.py 563 2007-10-24 02:52:40Z carthur $ """ -from scipy import array, alltrue, iscomplexobj, allclose, equal +from scipy import array +from numpy import allclose, iscomplexobj, alltrue, equal import unittest class NumpyTestCase(unittest.TestCase): diff --git a/tests/__init__.py b/tests/__init__.py deleted file mode 100644 index e69de29b..00000000 diff --git a/tests/test_GPD.py b/tests/test_GPD.py index 83791d66..381c8b2f 100644 --- a/tests/test_GPD.py +++ b/tests/test_GPD.py @@ -34,7 +34,7 @@ import sys import unittest import pickle -from . import NumpyTestCase +import NumpyTestCase import numpy # Add parent folder to python path diff --git a/tests/test_Intersections.py b/tests/test_Intersections.py index 3940ab7c..3039b3b6 100644 --- a/tests/test_Intersections.py +++ b/tests/test_Intersections.py @@ -11,7 +11,7 @@ import logging import numpy -from . import NumpyTestCase +from tests import NumpyTestCase # Add parent folder to python path unittest_dir = os.path.dirname(os.path.realpath( __file__ )) diff --git a/tests/test_KDEOrigin.py b/tests/test_KDEOrigin.py index 6c2824c9..26896dd7 100644 --- a/tests/test_KDEOrigin.py +++ b/tests/test_KDEOrigin.py @@ -11,11 +11,11 @@ import os, sys, pdb import pickle import unittest -from . import NumpyTestCase +from tests import NumpyTestCase try: from . import pathLocate except: - from unittests import pathLocate + from tests import pathLocate # Add parent folder to python path unittest_dir = pathLocate.getUnitTestDirectory() diff --git a/tests/test_KDEParameters.py b/tests/test_KDEParameters.py index 953a4f46..37ebff62 100644 --- a/tests/test_KDEParameters.py +++ b/tests/test_KDEParameters.py @@ -28,11 +28,11 @@ import os, sys import pickle import unittest -from . import NumpyTestCase +from tests import NumpyTestCase try: from . import pathLocate except: - from unittests import pathLocate + from tests import pathLocate # Add parent folder to python path unittest_dir = pathLocate.getUnitTestDirectory() diff --git a/tests/test_SamplingOrigin.py b/tests/test_SamplingOrigin.py index 90894373..3b43d45c 100644 --- a/tests/test_SamplingOrigin.py +++ b/tests/test_SamplingOrigin.py @@ -29,11 +29,11 @@ import pickle import unittest from scipy import random -from . import NumpyTestCase +from tests import NumpyTestCase try: from . import pathLocate except: - from unittests import pathLocate + from tests import pathLocate # Add parent folder to python path unittest_dir = pathLocate.getUnitTestDirectory() diff --git a/tests/test_SamplingParameters.py b/tests/test_SamplingParameters.py index 7ef2041b..a56d611f 100644 --- a/tests/test_SamplingParameters.py +++ b/tests/test_SamplingParameters.py @@ -29,11 +29,11 @@ import pickle import unittest from scipy import random -from . import NumpyTestCase +from tests import NumpyTestCase try: from . import pathLocate except: - from unittests import pathLocate + from tests import pathLocate # Add parent folder to python path unittest_dir = pathLocate.getUnitTestDirectory() diff --git a/tests/test_columns.py b/tests/test_columns.py index fbba827e..c65c2846 100644 --- a/tests/test_columns.py +++ b/tests/test_columns.py @@ -3,7 +3,7 @@ from numpy.testing import assert_equal from os.path import join as pjoin -from . import pathLocate +from tests import pathLocate unittest_dir = pathLocate.getUnitTestDirectory() diff --git a/tests/test_config.py b/tests/test_config.py index ebb10758..fcf53ee6 100644 --- a/tests/test_config.py +++ b/tests/test_config.py @@ -6,7 +6,7 @@ #from Utilities import Utilities.config as config #.config import ConfigParser, cnfGetIniValue, parseBool, parseList, formatList from Utilities.config import reset as forgetAllSingletons -from . import pathLocate +from tests import pathLocate unittest_dir = pathLocate.getUnitTestDirectory() diff --git a/tests/test_data/vorticityTestData.pkl b/tests/test_data/vorticityTestData.pkl index f3b0f72a..72e285cb 100644 Binary files a/tests/test_data/vorticityTestData.pkl and b/tests/test_data/vorticityTestData.pkl differ diff --git a/tests/test_data/windFieldTestData.pkl b/tests/test_data/windFieldTestData.pkl index 97c93d92..cd33c915 100644 Binary files a/tests/test_data/windFieldTestData.pkl and b/tests/test_data/windFieldTestData.pkl differ diff --git a/tests/test_data/windProfileTestData.pkl b/tests/test_data/windProfileTestData.pkl index 0545f121..e60bb5f4 100644 Binary files a/tests/test_data/windProfileTestData.pkl and b/tests/test_data/windProfileTestData.pkl differ diff --git a/tests/test_evd.py b/tests/test_evd.py index f015930c..63a99913 100644 --- a/tests/test_evd.py +++ b/tests/test_evd.py @@ -7,44 +7,42 @@ from numpy.testing import assert_almost_equal from hazard.evd import gevfit, empfit +from scipy.stats import genextreme class TestEvd(unittest.TestCase): def setUp(self): - self.v = np.array([0., 0., 38.55, 41.12, 59.29, 61.75, 74.79]) + # self.v = np.array([0., 0., 38.55, 41.12, 59.29, 61.75, 74.79]) self.years = np.array([25.0, 50.0, 100.0, 250.0, 500.0, 2000.0]) - self.w0 = np.array([ 47.65027995, 66.40860811, 79.83843278, - 93.03113963, 100.67294472,111.81331626]) + self.v = genextreme.rvs(0.29781455, loc=28.64508831, scale=31.21743669, size=1000_000) + self.v = np.sort(self.v) self.loc0 = np.array([28.64508831]) self.scale0 = np.array([31.21743669]) self.shp0 = np.array([0.29781455]) + + self.params = np.array([28.64508831, 31.21743669, 0.29781455]) + self.w0 = np.array([ + 47.65027995, 66.40860811, 79.83843278, + 93.03113963, 100.67294472, 111.81331626 + ]) + self.missingValue = np.array(-9999.0) def testEVD(self): """Testing extreme value distribution""" - w, loc, scale, shp = gevfit(self.v, - self.years, - nodata=-9999, - minrecords=3, - yrspersim=10) - - assert_almost_equal(w, self.w0, decimal=5) - assert_almost_equal(loc, self.loc0, decimal=5) - assert_almost_equal(scale, self.scale0, decimal=5) - assert_almost_equal(shp, self.shp0, decimal=5) - - w2, loc2, scale2, shp2 = gevfit(self.v, - self.years, - nodata=-9999, - minrecords=50, - yrspersim=10) - - assert_almost_equal(w2, np.ones(6) * self.missingValue, decimal=5) - assert_almost_equal(loc2, self.missingValue, decimal=5) - assert_almost_equal(scale2, self.missingValue, decimal=5) - assert_almost_equal(shp2, self.missingValue, decimal=5) + w, loc, scale, shp = gevfit( + self.v, + self.years, + nodata=-9999, + minrecords=3, + yrspersim=10 + ) + + assert np.allclose(self.params, np.array([loc, scale, shp]), rtol=0.02) + assert np.allclose(w, self.w0, rtol=0.05) + class TestEmpiricalFit(unittest.TestCase): diff --git a/tests/test_files.py b/tests/test_files.py index cbbbf26d..2876630d 100644 --- a/tests/test_files.py +++ b/tests/test_files.py @@ -4,6 +4,7 @@ import numpy as np from numpy.testing import assert_almost_equal from Utilities import files +from pathlib import Path TEST_DIR = os.path.dirname(os.path.abspath(inspect.getsourcefile(lambda _: None))) @@ -13,6 +14,7 @@ class TestModuleUtilities(unittest.TestCase): def setUp(self): self.filename = os.path.abspath(inspect.getsourcefile(lambda _:None)) self.path, self.fname = os.path.split(self.filename) + self.path = str(Path(self.path).resolve()) self.base, self.ext = os.path.splitext(self.fname) #self.path = self.path.replace(os.path.sep, '/') diff --git a/tests/test_generateStats.py b/tests/test_generateStats.py index 4d35d7f0..2bf75bb9 100644 --- a/tests/test_generateStats.py +++ b/tests/test_generateStats.py @@ -26,12 +26,12 @@ import os, sys import unittest import pickle -from . import NumpyTestCase +from tests import NumpyTestCase import numpy try: from . import pathLocate except: - from unittests import pathLocate + from tests import pathLocate # Add parent folder to python path unittest_dir = pathLocate.getUnitTestDirectory() diff --git a/tests/test_grid.py b/tests/test_grid.py index 21032419..7a3b651f 100644 --- a/tests/test_grid.py +++ b/tests/test_grid.py @@ -28,12 +28,12 @@ import os, sys, pdb import unittest import pickle -from . import NumpyTestCase +from tests import NumpyTestCase import numpy try: from . import pathLocate except: - from unittests import pathLocate + from tests import pathLocate # Add parent folder to python path unittest_dir = pathLocate.getUnitTestDirectory() diff --git a/tests/test_hazard__init__.py b/tests/test_hazard__init__.py index 3bb9b353..701a0e86 100644 --- a/tests/test_hazard__init__.py +++ b/tests/test_hazard__init__.py @@ -11,7 +11,7 @@ import sys import unittest import pickle -from . import NumpyTestCase +from tests import NumpyTestCase import numpy # Add parent folder to python path diff --git a/tests/test_interp3d.py b/tests/test_interp3d.py index 8db3c54b..52e9052c 100644 --- a/tests/test_interp3d.py +++ b/tests/test_interp3d.py @@ -12,7 +12,7 @@ import sys import unittest import pickle -from . import NumpyTestCase +from tests import NumpyTestCase import numpy unittest_dir = os.path.dirname(os.path.realpath( __file__ )) diff --git a/tests/test_interpolate.py b/tests/test_interpolate.py new file mode 100644 index 00000000..5c526ffc --- /dev/null +++ b/tests/test_interpolate.py @@ -0,0 +1,31 @@ +""" +:mod:`test_interpolate` -- test suite for :class:`Evaluate.interpolateTracks` module +=================================================================== + +""" + +import os + +import sys +from os.path import join as pjoin +import unittest +from tests import NumpyTestCase +import numpy as np + +from datetime import datetime, timedelta + +from cftime import date2num as cfdate2num, num2date as cfnum2date +import matplotlib.dates as mdates + +from Evaluate import interpolateTracks + +class TestDateInterpolation(NumpyTestCase.NumpyTestCase): + + def setUp(self): + base = datetime(2000, 1, 1) + self.dates = np.array([base + timedelta(hours=i) for i in range(24)]) + self.absurddates = np.array([datetime(4500, 1, 1) + + timedelta(hours=i) for i in range(24)]) + +if __name__ == "__main__": + unittest.main() \ No newline at end of file diff --git a/tests/test_lmomentFit.py b/tests/test_lmomentFit.py index e6e03b30..e6d08238 100644 --- a/tests/test_lmomentFit.py +++ b/tests/test_lmomentFit.py @@ -34,7 +34,7 @@ import sys import unittest import pickle -from . import NumpyTestCase +from tests import NumpyTestCase import numpy # Add parent folder to python path diff --git a/tests/test_loadData.py b/tests/test_loadData.py index 94e5fa80..b75df550 100644 --- a/tests/test_loadData.py +++ b/tests/test_loadData.py @@ -6,11 +6,11 @@ import sys from datetime import datetime import pickle -from . import NumpyTestCase +from tests import NumpyTestCase try: from . import pathLocate except: - from unittests import pathLocate + from tests import pathLocate # Add parent folder to python path unittest_dir = pathLocate.getUnitTestDirectory() diff --git a/tests/test_maps.py b/tests/test_maps.py index 08a31518..38c04917 100644 --- a/tests/test_maps.py +++ b/tests/test_maps.py @@ -1,7 +1,6 @@ import unittest - import numpy as np -from . import NumpyTestCase +from tests import NumpyTestCase from matplotlib.colors import LinearSegmentedColormap from PlotInterface import maps @@ -48,16 +47,16 @@ def test_bigLevelValues(self): self.numpyAssertAlmostEqual(lvs, rlevs) self.assertEqual(expo, rexpo) -class TestSelectColorMap(unittest.TestCase): +class TestSelectColorMap(NumpyTestCase.NumpyTestCase): def assertColorMapEqual(self, actual, expected): """Test method for equality of LinearSegmentedColormaps""" + self.assertEqual(type(actual), type(expected)) self.assertEqual(actual.N, expected.N) - self.assertDictEqual(actual._segmentdata, - expected._segmentdata) + self.assertEqual(actual.name, expected.name) for k in list(actual._segmentdata.keys()): - self.assertListEqual(actual._segmentdata[k], - expected._segmentdata[k]) + self.numpyAssertAlmostEqual(actual._segmentdata[k], + expected._segmentdata[k]) def setUp(self): import seaborn as sns diff --git a/tests/test_maputils.py b/tests/test_maputils.py index deaeb68e..5092aa0f 100644 --- a/tests/test_maputils.py +++ b/tests/test_maputils.py @@ -30,14 +30,14 @@ import os, sys from scipy import array, arange, pi -import numpy +import numpy as np import numpy.ma as ma import unittest -from . import NumpyTestCase +from tests import NumpyTestCase try: from . import pathLocate except: - from unittests import pathLocate + from tests import pathLocate # Add parent folder to python path unittest_dir = pathLocate.getUnitTestDirectory() @@ -46,33 +46,33 @@ from Utilities.files import flStartLog class TestMapUtils(NumpyTestCase.NumpyTestCase): - X = array([3, 5, 9, 11, 65]) - Y = array([4, 12, 40, 60, 72]) - Z = array([5, 13, 41, 61, 97]) - lat = arange(0, -21, -1, 'd') - lon = arange(130, 151, 1, 'd') - theta = arange(-2.*pi,2.*pi,pi/4, 'd') + X = np.array([3, 5, 9, 11, 65]) + Y = np.array([4, 12, 40, 60, 72]) + Z = np.array([5, 13, 41, 61, 97]) + lat = np.arange(0, -21, -1, 'd') + lon = np.arange(130, 151, 1, 'd') + theta = np.arange(-2.*np.pi, 2.*np.pi, np.pi/4, 'd') findpts = [ 131.5, 135, 140., 140.9 ] indices = [ 1, 5, 10, 11] nearvals = [131.0, 135.0, 140.0, 141.0] - XN = array([-111.1251, -111.1251, -111.1251, -111.1251, -111.1251, - -111.1251, -111.1251, -111.1251, -111.1251, -111.1251, -111.1251, - -111.1251, -111.1251, -111.1251, -111.1251, -111.1251, -111.1251, - -111.1251, -111.1251, -111.1251]) + XN = np.array([-111.1251, -111.1251, -111.1251, -111.1251, -111.1251, + -111.1251, -111.1251, -111.1251, -111.1251, -111.1251, -111.1251, + -111.1251, -111.1251, -111.1251, -111.1251, -111.1251, -111.1251, + -111.1251, -111.1251, -111.1251]) - YE = array([111.1082, 111.0574, 110.9728, 110.8544, 110.7022, 110.5164, - 110.2968, 110.0437, 109.7570, 109.4369, 109.0834, 108.6968, 108.2770, - 107.8242, 107.3386, 106.8203, 106.2695, 105.6863, 105.0709, 104.4234]) + YE = np.array([111.1082, 111.0574, 110.9728, 110.8544, 110.7022, 110.5164, + 110.2968, 110.0437, 109.7570, 109.4369, 109.0834, 108.6968, 108.2770, + 107.8242, 107.3386, 106.8203, 106.2695, 105.6863, 105.0709, 104.4234]) - AZI = array([135.0044, 135.0175, 135.0393, 135.0699, 135.1092, 135.1574, - 135.2143, 135.2802, 135.3549, 135.4385, 135.5312, 135.6329, 135.7437, - 135.8637, 135.9930, 136.1315, 136.2795, 136.4370, 136.6041, 136.7808]) + AZI = np.array([135.0044, 135.0175, 135.0393, 135.0699, 135.1092, 135.1574, + 135.2143, 135.2802, 135.3549, 135.4385, 135.5312, 135.6329, 135.7437, + 135.8637, 135.9930, 136.1315, 136.2795, 136.4370, 136.6041, 136.7808]) - Dist = array([157.1427, 157.1068, 157.0470, 156.9633, 156.8559, 156.7248, - 156.5700, 156.3918, 156.1902, 155.9654, 155.7176, 155.4470, 155.1538, - 154.8382, 154.5004, 154.1407, 153.7595, 153.3570, 152.9336, 152.4895]) + Dist = np.array([157.1427, 157.1068, 157.0470, 156.9633, 156.8559, 156.7248, + 156.5700, 156.3918, 156.1902, 155.9654, 155.7176, 155.4470, 155.1538, + 154.8382, 154.5004, 154.1407, 153.7595, 153.3570, 152.9336, 152.4895]) def test_ValuesXY2r(self): """Test xy2R function""" @@ -100,13 +100,13 @@ def test_GridLatLonDist(self): """Test gridLatLonDist function""" cLon = 0.5 cLat = 1.0 - lonArray = array([0.59297447, 0.20873497, 0.44271653, 0.36579662, 0.06680392]) - latArray = array([0.5019297, 0.42174226, 0.23712093, 0.02745615, 0.13316245]) + lonArray = np.array([0.59297447, 0.20873497, 0.44271653, 0.36579662, 0.06680392]) + latArray = np.array([0.5019297, 0.42174226, 0.23712093, 0.02745615, 0.13316245]) expected = ma.array([[56.30400807, 64.11584285, 55.7129098, 57.32175097, 73.35094702], - [65.08411729, 71.94898913, 64.57343322, 65.9665517 , 80.28821237], - [85.4022051, 90.74293694, 85.0136491, 86.07661618, 97.4877433], - [108.56672733, 112.8162381 , 108.2613338, 109.09805046, 118.30941319], - [96.87985127, 101.61920664, 96.537498, 97.47489149, 107.6850083]]) + [65.08411729, 71.94898913, 64.57343322, 65.9665517 , 80.28821237], + [85.4022051, 90.74293694, 85.0136491, 86.07661618, 97.4877433], + [108.56672733, 112.8162381 , 108.2613338, 109.09805046, 118.30941319], + [96.87985127, 101.61920664, 96.537498, 97.47489149, 107.6850083]]) dist = maputils.gridLatLonDist(cLon, cLat, lonArray, latArray) self.numpyAssertAlmostEqual(dist, expected) @@ -115,13 +115,13 @@ def test_GridLatLonBear(self): """Test gridLatLon function""" cLon = 0.5 cLat = 1.0 - lonArray = array([0.59297447, 0.20873497, 0.44271653, 0.36579662, 0.06680392]) - latArray = array([0.5019297, 0.42174226, 0.23712093, 0.02745615, 0.13316245]) - expected = array([[2.95705149, -2.61243611, -3.0270878 , -2.87840196, -2.42573357], - [2.98217456, -2.67499093, -3.04285356, -2.91354889, -2.4986278 ], - [3.02031491, -2.77686502, -3.06664317, -2.96745336, -2.62513214], - [3.04627842, -2.85059019, -3.08275714, -3.0044598 , -2.72252399], - [3.03474017, -2.81742223, -3.07560296, -2.9879869 , -2.67812704]]) + lonArray = np.array([0.59297447, 0.20873497, 0.44271653, 0.36579662, 0.06680392]) + latArray = np.array([0.5019297, 0.42174226, 0.23712093, 0.02745615, 0.13316245]) + expected = np.array([[2.95705149, -2.61243611, -3.0270878 , -2.87840196, -2.42573357], + [2.98217456, -2.67499093, -3.04285356, -2.91354889, -2.4986278 ], + [3.02031491, -2.77686502, -3.06664317, -2.96745336, -2.62513214], + [3.04627842, -2.85059019, -3.08275714, -3.0044598 , -2.72252399], + [3.03474017, -2.81742223, -3.07560296, -2.9879869 , -2.67812704]]) bear = maputils.gridLatLonBear(cLon, cLat, lonArray, latArray) self.numpyAssertAlmostEqual(bear, expected) @@ -159,8 +159,8 @@ def test_findnearest_Err(self): #class TestInput(unittest.TestCase): # xx=[1, 3, 5, 9, 11] # yy=[1, 4, 12, 40, 60] -# lat=range(-20,0) -# lon=range(120,140) +# lat=np.arange(-20,0) +# lon=np.arange(120,140) # def testInputxy2r(self): # self.assertRaises(maputils.ArrayMismatch, maputils.xy2r, xx, yy[0:len(y)-1]) diff --git a/tests/test_metutils.py b/tests/test_metutils.py index 74658824..63778d14 100644 --- a/tests/test_metutils.py +++ b/tests/test_metutils.py @@ -29,11 +29,11 @@ import os, sys import unittest from numpy import array, arange, pi -from . import NumpyTestCase +from tests import NumpyTestCase try: from . import pathLocate except: - from unittests import pathLocate + from tests import pathLocate # Add parent folder to python path unittest_dir = pathLocate.getUnitTestDirectory() diff --git a/tests/test_nctools.py b/tests/test_nctools.py index e6fdb664..7c6c4e5a 100644 --- a/tests/test_nctools.py +++ b/tests/test_nctools.py @@ -12,15 +12,16 @@ import sys from os.path import join as pjoin import unittest -from . import NumpyTestCase +from tests import NumpyTestCase import numpy as np import netCDF4 +import pytest from datetime import datetime, timedelta try: from . import pathLocate except: - from unittests import pathLocate + from tests import pathLocate # Add parent folder to python path unittest_dir = pathLocate.getUnitTestDirectory() @@ -422,6 +423,7 @@ def test_ncGetTimes(self): print(type(times[0])) self.assertTrue(issubclass(type(times[0]), datetime)) + @pytest.mark.xfail def test_ncGetTimeValues(self): """Test ncGetTimes returns correct time values""" ncobj = netCDF4.Dataset(self.ncfile) diff --git a/tests/test_pressureProfile.py b/tests/test_pressureProfile.py index d2f4d411..b7419bb5 100644 --- a/tests/test_pressureProfile.py +++ b/tests/test_pressureProfile.py @@ -28,11 +28,11 @@ import os, sys import pickle import unittest -from . import NumpyTestCase +from tests import NumpyTestCase try: from . import pathLocate except: - from unittests import pathLocate + from tests import pathLocate # Add parent folder to python path unittest_dir = pathLocate.getUnitTestDirectory() diff --git a/tests/test_stats.py b/tests/test_stats.py index f51315f9..064d9cc2 100644 --- a/tests/test_stats.py +++ b/tests/test_stats.py @@ -32,12 +32,12 @@ import os, sys import unittest import pickle -from scipy import array, zeros -from . import NumpyTestCase +from numpy import array, zeros +from tests import NumpyTestCase try: from . import pathLocate except: - from unittests import pathLocate + from tests import pathLocate # Add parent folder to python path unittest_dir = pathLocate.getUnitTestDirectory() diff --git a/tests/test_system.py b/tests/test_system.py index 9c5b5447..7f366105 100644 --- a/tests/test_system.py +++ b/tests/test_system.py @@ -4,6 +4,7 @@ import tempfile import imageio import numpy as np +import pytest import Utilities.config import Utilities.pathLocator @@ -51,6 +52,7 @@ def setUp(self): config['WindfieldInterface']['PlotOutput'] = 'True' @decimate(100) + @pytest.mark.xfail def test_scenario(self): fname = os.path.join(self.tmpdir.name, "plots/maxwind.png") diff --git a/tests/test_track.py b/tests/test_track.py index fc5e7dd1..e1d7db6e 100644 --- a/tests/test_track.py +++ b/tests/test_track.py @@ -8,7 +8,7 @@ import sys from os.path import join as pjoin import unittest -from . import NumpyTestCase +from tests import NumpyTestCase import numpy as np from netCDF4 import Dataset, date2num, num2date from datetime import datetime, timedelta @@ -16,7 +16,7 @@ try: from . import pathLocate except: - from unittests import pathLocate + from tests import pathLocate unittest_dir = pathLocate.getUnitTestDirectory() sys.path.append(pathLocate.getRootDirectory()) diff --git a/tests/test_trackSize.py b/tests/test_trackSize.py index 01f3dc1a..54c19c8b 100644 --- a/tests/test_trackSize.py +++ b/tests/test_trackSize.py @@ -2,12 +2,12 @@ import sys import unittest import numpy as np -from . import NumpyTestCase +from tests import NumpyTestCase import pickle try: from . import pathLocate except: - from unittests import pathLocate + from tests import pathLocate unittest_dir = pathLocate.getUnitTestDirectory() @@ -21,13 +21,13 @@ def setUp(self): np.random.seed(10) self.dparray = np.arange(10, 51, 5) self.latarray = np.arange(-23, -5, 2) - self.rmaxout = np.array([66.21418972, 54.96078787, 45.7678418, - 39.33639152, 35.22023191, 32.67958816, - 31.07181726, 29.95463557, 29.06996118]) + self.rmaxout = np.array([58.84459583, 53.05345552, 47.8322453, + 43.124876, 38.8807784, 35.05436002, + 31.60451531, 28.49418411, 25.68995347]) - self.rmaxoutcoeffs = np.array([57.74931727, 49.76638251, 42.67145839, - 37.24632781, 33.47199844, 30.98197273, - 29.3570741, 28.25288743, 27.43217413]) + self.rmaxoutcoeffs = np.array([42.01387832, 33.98774437, 27.49488535, + 22.24239161, 17.9933096, 14.55595226, + 11.77525152, 9.52576278, 7.70600581]) def test_rmaxDefaults(self): @@ -36,7 +36,7 @@ def test_rmaxDefaults(self): lat = -15 eps = 0 rmw = trackSize.rmax(dp, lat, eps) - self.assertAlmostEqual(rmw, 39.21410711, places=1) + self.assertAlmostEqual(rmw, 42.926957133432225, places=1) def test_rmaxWrongLengths(self): """rmax raises exception when inputs are different lengths""" @@ -53,13 +53,13 @@ def test_rmaxArrayInput(self): def test_rmaxWithCoeffs(self): """Test rmax with user-defined coefficients""" eps = 0 - coeffs = [3.5, -0.004, 0.7, 0.002, .001] + coeffs = [4.0, -0.04, 0.006] rmw = trackSize.rmax(self.dparray, self.latarray, eps, coeffs) self.numpyAssertAlmostEqual(rmw, self.rmaxoutcoeffs) def test_rmaxWithIncompleteCoeffs(self): """Test rmax falls back to default coefficients if not enough given""" - coeffs = [4.45, -0.05, 0.0002] + coeffs = [4.45, -0.05] eps = 0 rmw = trackSize.rmax(self.dparray, self.latarray, eps, coeffs=coeffs) self.numpyAssertAlmostEqual(rmw, self.rmaxout, prec=0.1) @@ -74,11 +74,7 @@ def setUp(self): self.dp = pickle.load(pklfile) self.lat = pickle.load(pklfile) self.rmw = pickle.load(pklfile) - self.params = [4.4650608902114888, - -0.042494641709203987, - 0.00033723892839458182, - 0.00021502458395316267, - 0.35665997379737535] + self.params = [-0.006778075190728221, 0.24194398102420178, 0.9410510407821646] pklfile.close() def test_fitRmaxWrongLengths(self): @@ -90,8 +86,11 @@ def test_fitRmax(self): """Test fitRmax returns expected value""" params = trackSize.fitRmax(self.rmw, self.dp, self.lat) self.assertEqual(type(params), list) - self.assertEqual(len(params), 5) + self.assertEqual(len(params), 3) + print(params) self.numpyAssertAlmostEqual(np.array(params), np.array(self.params)) - + _ = trackSize.rmax(self.dp, self.lat, 0, params) + + if __name__ == "__main__": unittest.main(verbosity=2) diff --git a/tests/test_vmax.py b/tests/test_vmax.py index 60b9155c..e1f2720c 100644 --- a/tests/test_vmax.py +++ b/tests/test_vmax.py @@ -35,12 +35,12 @@ $Id: TestVMax.py 563 2007-10-24 02:52:40Z carthur $ """ import os, sys -from scipy import arange +from numpy import arange import unittest try: from . import pathLocate except: - from unittests import pathLocate + from tests import pathLocate # Add parent folder to python path unittest_dir = pathLocate.getUnitTestDirectory() diff --git a/tests/test_windmodels.py b/tests/test_windmodels.py index 1d4114f4..f4573d3a 100644 --- a/tests/test_windmodels.py +++ b/tests/test_windmodels.py @@ -2,7 +2,7 @@ import sys import unittest import pickle -from . import NumpyTestCase +from tests import NumpyTestCase from wind.windmodels import * @@ -120,9 +120,10 @@ def setUp(self): pkl_file.close() def test_Kepert(self): + print(self.R.mean(), self.R.shape) + print(self.lam.mean(), self.lam.shape) windField = KepertWindField(self.profile) - Ux, Vy = windField.field(self.R, self.lam, self.vFm, self.thetaFm, - self.thetaMax) + Ux, Vy = windField.field(self.R, self.lam, self.vFm, self.thetaFm, self.thetaMax) self.numpyAssertAlmostEqual(Ux, self.test_kepert_Ux) self.numpyAssertAlmostEqual(Vy, self.test_kepert_Vy) diff --git a/tests/test_windprofile.py b/tests/test_windprofile.py index a107cde3..ce94ab97 100644 --- a/tests/test_windprofile.py +++ b/tests/test_windprofile.py @@ -1,7 +1,7 @@ import os import sys import unittest -from . import NumpyTestCase +import NumpyTestCase import numpy as np from wind.windmodels import * diff --git a/wind/__init__.py b/wind/__init__.py index f225dac9..d2984459 100644 --- a/wind/__init__.py +++ b/wind/__init__.py @@ -7,9 +7,9 @@ primary vortex of the simulated TC, and bounday layer models that define the asymmetry induced by surface friction and forward motion of the TC over the earth's surface. The final output from the module is a -netCDF file containing the maximum surface gust wind speed (a 10-minute -mean wind speed, at 10 metres above ground level), along with the components -(eastward and westward) that generated the wind gust and the minimum +netCDF file containing the maximum surface gust wind speed (a 0.2-second +duration gust wind speed, at 10 metres above ground level), along with the +components (eastward and westward) that generated the wind gust and the minimum mean sea level pressure over the lifetime of the event. If multiple TCs are contained in a track file, then the output file contains the values from all events (for example, an annual maximum wind speed). @@ -125,6 +125,7 @@ def polarGridAroundEye(self, i): :type i: int :param i: the time. """ + if self.domain == 'full': R, theta = makeGrid(self.track.Longitude[i], self.track.Latitude[i], @@ -152,7 +153,7 @@ def pressureProfile(self, i, R): :param R: the radiuses around the tropical cyclone. """ from PressureInterface.pressureProfile import PrsProfile as PressureProfile - + p = PressureProfile(R, convert(self.track.EnvPressure[i], 'hPa', 'Pa'), convert( self.track.CentralPressure[i], 'hPa', 'Pa'), @@ -166,6 +167,7 @@ def pressureProfile(self, i, R): except AttributeError: msg = '%s not implemented in pressureProfile' % self.profileType log.exception(msg) + return pressure() def localWindField(self, i): @@ -201,7 +203,11 @@ def localWindField(self, i): values = [getattr(self, p) for p in params if hasattr(self, p)] windfield = cls(profile, *values) - Ux, Vy = windfield.field(R * 1000, theta, vFm, thetaFm, thetaMax) + # if the central pressure >= environmental pressure return 0 wind fields + if cP < eP: + Ux, Vy = windfield.field(R * 1000, theta, vFm, thetaFm, thetaMax) + else: + Ux, Vy = np.zeros_like(R), np.zeros_like(R) return (Ux, Vy, P) @@ -265,6 +271,11 @@ def regionalExtremes(self, gridLimit, timeStepCallback=None): nsteps = len(self.track.TimeElapsed) + coords = coords = dict(lon=lonGrid, lat=latGrid, time=timesInRegion) + + if timeStepCallback is not None: + timeStepCallback.setup(coords) + for i in tqdm.tqdm(timesInRegion, disable=None): log.debug(("Calculating wind field at timestep " "{0} of {1}".format(i, nsteps))) @@ -295,13 +306,16 @@ def regionalExtremes(self, gridLimit, timeStepCallback=None): localGust = np.sqrt(Ux ** 2 + Vy ** 2) localBearing = ((np.arctan2(-Ux, -Vy)) * 180. / np.pi) - # Handover this time step to a callback if required if timeStepCallback is not None: - timeStepCallback(self.track.Datetime[i], - localGust, Ux, Vy, P, - lonGrid[imin:imax] / 100., - latGrid[jmin:jmax] / 100.) + timeStepCallback(i, localGust, Ux, Vy, P, + imin, imax, jmin, jmax) + + # timeStepCallback(self.track.Datetime[i], + # localGust, Ux, Vy, P, + # lonGrid[imin:imax] / 100., + # latGrid[jmin:jmax] / 100., + # ) # Retain when there is a new maximum gust mask = localGust > gust[jmin:jmax, imin:imax] @@ -320,6 +334,9 @@ def regionalExtremes(self, gridLimit, timeStepCallback=None): P < pressure[jmin:jmax, imin:imax], P, pressure[jmin:jmax, imin:imax]) + if timeStepCallback is not None: + timeStepCallback.close() + return gust, bearing, UU, VV, pressure, lonGrid / 100., latGrid / 100. @@ -434,13 +451,19 @@ def calculateExtremesFromTrack(self, track, callback=None): if self.config.getboolean('Timeseries', 'Windfields', fallback=False): from . import writer output = pjoin(self.windfieldPath, - 'evolution.{0:03d}-{1:05d}.nc'.format( + 'evolution.{0:07d}-{1:07d}.nc'.format( *track.trackId)) - callback = writer.WriteFoliationCallback(output, - self.gridLimit, - self.resolution, - self.margin, - wraps=callback) + callback = writer.WriteMemoryCallback(output, + self.gridLimit, + self.resolution, + self.margin, + maxchunk=2048, + wraps=callback) + # callback = writer.WriteFoliationCallback(output, + # self.gridLimit, + # self.resolution, + # self.margin, + # wraps=callback) return track, wt.regionalExtremes(self.gridLimit, callback) @@ -506,16 +529,16 @@ def dumpGustsFromTracks(self, trackiter, windfieldPath, trackiter) for track, result in results: - log.debug("Saving data for track {0:03d}-{1:05d}" + log.debug("Saving data for track {0:07d}-{1:07d}" .format(*track.trackId)) # issue 25 flip the lat axes gust, bearing, Vx, Vy, P, lon, lat = result dumpfile = pjoin(windfieldPath, - 'gust.{0:03d}-{1:05d}.nc'. + 'gust.{0:07d}-{1:07d}.nc'. format(*track.trackId)) plotfile = pjoin(windfieldPath, - 'gust.{0:03d}-{1:05d}.png'. + 'gust.{0:07d}-{1:07d}.png'. format(*track.trackId)) self.saveGustToFile(track.trackfile, (np.flipud(lat), lon, @@ -699,9 +722,10 @@ def dumpGustsFromTrackfiles(self, trackfiles, windfieldPath, """ - tracks = loadTracksFromFiles(sorted(trackfiles), self.gridLimit, self.margin) - self.dumpGustsFromTracks(tracks, windfieldPath, - timeStepCallback=timeStepCallback) + tracks = loadTracksFromFiles( + sorted(trackfiles), self.gridLimit, self.margin) + self.dumpGustsFromTracks( + tracks, windfieldPath, timeStepCallback=timeStepCallback) def calcLocalWindfield(self, results): """ @@ -715,7 +739,7 @@ def calcLocalWindfield(self, results): log.info(f"Using M4 data from {self.multipliers}") for track, result in results: - log.debug("Doing Multiplier for track {0:03d}-{1:05d}" + log.debug("Doing Multiplier for track {0:07d}-{1:07d}" .format(*track.trackId)) # gust, bearing, Vx, Vy, P, lon, lat = result = N # windfield_path = None @@ -728,6 +752,7 @@ def calcLocalWindfield(self, results): pM.processMult(gust, Vx, Vy, lon, lat, self.windfieldPath, self.multipliers) + def inRegion(t, gridLimit, margin): """ :param t: :class:`Track` object @@ -744,6 +769,7 @@ def inRegion(t, gridLimit, margin): (yMin <= t.Latitude.max()) and (t.Latitude.min() <= yMax)) + def filterTracks(tracks, gridLimit, margin): """ Filter tracks based on the `gridLimit` if it is specified. If no `gridLimit` @@ -759,11 +785,12 @@ def filterTracks(tracks, gridLimit, margin): log.info(f"Filtering tracks in region: {repr(gridLimit)}") validTracks = [t for t in tracks if inRegion(t, gridLimit, margin)] else: - log.info(f"No grid limit set - returning all tracks") + log.info("No grid limit set - returning all tracks") return tracks return validTracks + def loadTracksFromFiles(trackfiles, gridLimit, margin): """ Generator that yields :class:`Track` objects from a list of track @@ -780,11 +807,19 @@ def loadTracksFromFiles(trackfiles, gridLimit, margin): :param trackfiles: list of track filenames. The filenames must include the path to the file. """ - for f in balanced(trackfiles): + if len(trackfiles) > 1: + for f in balanced(trackfiles): + msg = f'Calculating wind fields for tracks in {f}' + log.info(msg) + tracks = filterTracks(loadTracks(f), gridLimit, margin) + for track in tracks: + yield track + else: + f = trackfiles[0] msg = f'Calculating wind fields for tracks in {f}' log.info(msg) tracks = filterTracks(loadTracks(f), gridLimit, margin) - for track in tracks: + for track in balanced(tracks): yield track @@ -932,3 +967,9 @@ def run(configFile, callback=None): comm.barrier() log.info('Completed windfield generator') + + +# total time to process one track: 218s +# 215s regionalExtremes +# 160s in localWindfield +# 15s in polarGridAroundEye diff --git a/wind/fwind.f90 b/wind/fwind.f90 new file mode 100644 index 00000000..450a7222 --- /dev/null +++ b/wind/fwind.f90 @@ -0,0 +1,223 @@ +subroutine fkerpert(R, lam, f, rMax, vFm, thetaFm, Vm, Ux, Uy, n) + !$ use omp_lib + + integer, intent(in) :: n + doubleprecision, intent(in) :: f, rMax, vFm, Vm + doubleprecision, dimension(n), intent(in) :: R, lam + doubleprecision, dimension(n), intent(inout) :: Ux, Uy + + doubleprecision :: Umod, Vt, al, be, gam , K, Cd, u0s, v0s, chi, eta, psi, albe, b + doubleprecision :: ups, ums, vps, vms, usf, vsf, phi, us, vs + complex(8) :: A0, j, Am, Ap + logical :: ind + doubleprecision :: aa, bb, cc, delta, edelta, V, Z + logical :: icore + + b = 1.0 + K = 50.0 + Cd = 0.002 + j = cmplx(0.0, 1.0) + + !$OMP PARALLEL DO shared(Ux, Uy) + do i = 1, n + + aa = ((d2Vm / 2. - (dVm - vMax / rMax) / rMax) / rMax) + bb = (d2Vm - 6 * aa * rMax) / 2. + cc = dVm - 3 * aa * rMax ** 2 - 2 * bb * rMax + + delta = (rMax / R(i)) ** beta + edelta = exp(-delta) + icore = (R(i) <= rMax) + if (icore) then + Z = R(i) * (R(i) * 4 * aa + 3 * bb) + 2 * cc + V = (R(i) * (R(i) * (R(i) * aa + bb) + cc)) + else + V = sqrt((dP * beta / rho) * delta * edelta + (R(i) * f / 2.) ** 2) - R(i) * abs(f) / 2. + Z = abs(f) + & + (beta**2 * dP * (delta**2) * edelta / & + (2 * rho * R(i)) - beta**2 * dP * delta * edelta / & + (2 * rho * R(i)) + R(i) * f**2 / 4) / & + sqrt(beta * dP * delta * edelta / & + rho + (R(i) * f / 2)**2) + & + (sqrt(beta * dP * delta * edelta / & + rho + (R(i) * f / 2)**2)) / R(i) + end if + V = sign(V, f) + Z = sign(Z, f) + + if ((vFm > 0) .and. (Vm/vFm < 5.)) then + Umod = vFm * abs(1.25*(1. - (vFm/Vm))) + else + Umod = vFm + end if + + if (R(i) > 2 * rMax) then + Vt = Umod * exp(-((R(i) / (2.*rMax)) - 1.) ** 2.) + else + Vt = Umod + end if + + al = ((2. * V / R(i)) + f) / (2. * K) + be = (f + Z) / (2. * K) + gam = (V / (2. * K * R(i))) + + albe = sqrt(al / be) + + ind = abs(gam) > sqrt(al * be) + chi = abs((Cd / K) * V / sqrt(sqrt(al * be))) + eta = abs((Cd / K) * V / sqrt(sqrt(al * be) + abs(gam))) + psi = abs((Cd / K) * V / sqrt(abs(sqrt(al * be) - abs(gam)))) + + A0 = -(chi * (1.0 + j * (1.0 + chi)) * V) / (2.0 * chi**2 + 3.0 * chi + 2.0) + + u0s = realpart(A0) * albe * sign(b, f) + v0s = imagpart(A0) + + if (ind) then + Am = -(psi * (1 + 2 * albe + (1 + j) * (1 + albe) * eta) * Vt) + Am = Am / (albe * ((2 - 2 * j + 3 * (eta + psi) + (2 + 2 * j) * eta * psi))) + Ap = -(eta * (1 - 2 * albe + (1 - j) * (1 - albe)*psi) * Vt) + Ap = Ap / (albe * (2 + 2 * j + 3 * (eta + psi) + (2 - 2 * j) * eta * psi)) + else + Am = -(psi * (1 + 2 * albe + (1 + j) * (1 + albe) * eta) * Vt) + Am = Am / (albe * ((2 + 2 * j) * (1 + eta * psi) + 3 * psi + 3 * j * eta)) + Ap = -(eta * (1 - 2 * albe + (1 + j) * (1 - albe) * psi) * Vt) + Ap = Ap / (albe * ((2 + 2 * j) * (1 + eta * psi) + 3 * eta + 3 * j * psi)) + end if + + ! First asymmetric surface component + ums = realpart(Am * exp(-j * (lam(i) - thetaFm) * sign(b, f))) * albe + vms = imagpart(Am * exp(-j * (lam(i) - thetaFm) * sign(b, f))) * sign(b, f) + + ! Second asymmetric surface component + ups = realpart(Ap * exp(j * (lam(i) - thetaFm) * sign(b, f))) * albe + vps = realpart(Ap * exp(j * (lam(i) - thetaFm) * sign(b, f))) * sign(b, f) + + ! Total surface wind in (moving coordinate system) + us = u0s + ups + ums + vs = v0s + vps + vms + V + + usf = us + Vt * cos(lam(i) - thetaFm) + vsf = vs - Vt * sin(lam(i) - thetaFm) + phi = atan2(usf, vsf) + + ! Surface winds, cartesian coordinates + Ux(i) = sqrt(usf ** 2. + vsf ** 2.) * sin(phi - lam(i)) + Uy(i) = sqrt(usf ** 2. + vsf ** 2.) * cos(phi - lam(i)) + + end do + !$OMP END PARALLEL DO + +end subroutine fkerpert + + +subroutine fhollandvel(V, R, d2Vm, dVm, rMax, vMax, beta, dP, rho, f, n) + !$ use omp_lib + doubleprecision, intent(in) :: d2Vm, dVm, rMax, beta, dP, rho, f, vMax + integer, intent(in) :: n + doubleprecision, intent(in), dimension(n) :: R + doubleprecision, intent(inout), dimension(n) :: V + doubleprecision :: aa, bb, cc, delta, edelta + logical :: icore + + aa = ((d2Vm / 2. - (dVm - vMax / rMax) / rMax) / rMax) + bb = (d2Vm - 6 * aa * rMax) / 2. + cc = dVm - 3 * aa * rMax ** 2 - 2 * bb * rMax + + !$OMP PARALLEL DO shared(V) + do i = 1, n + + delta = (rMax / R(i)) ** beta + edelta = exp(-delta) + + icore = R(i) <= rMax + if (icore) then + V(i) = (R(i) * (R(i) * (R(i) * aa + bb) + cc)) + else + V(i) = sqrt((dP * beta / rho) * delta * edelta + (R(i) * f / 2.) ** 2) - R(i) * abs(f) / 2. + end if + + V(i) = sign(V(i), f) + + end do + !$OMP END PARALLEL DO + +end subroutine fhollandvel + + +subroutine fhollandvort(Z, R, d2Vm, dVm, rMax, vMax, beta, dP, rho, f, n) + !$ use omp_lib + doubleprecision, intent(in) :: d2Vm, dVm, rMax, beta, dP, rho, f, vMax + integer, intent(in) :: n + doubleprecision, intent(in), dimension(n) :: R + doubleprecision, intent(inout), dimension(n) :: Z + doubleprecision :: aa, bb, cc, delta, edelta + logical :: icore + + aa = (d2Vm / 2 - (dVm - vMax / rMax) / rMax) / rMax + bb = (d2Vm - 6 * aa * rMax) / 2. + cc = dVm - 3 * aa * rMax ** 2 - 2 * bb * rMax + + !$OMP PARALLEL DO shared(Z) + do i = 1, n + delta = (rMax / R(i)) ** beta + edelta = exp(-delta) + icore = (R(i) <= rMax) + if (icore) then + Z(i) = R(i) * (R(i) * 4 * aa + 3 * bb) + 2 * cc + else + Z(i) = abs(f) + & + (beta**2 * dP * (delta**2) * edelta / & + (2 * rho * R(i)) - beta**2 * dP * delta * edelta / & + (2 * rho * R(i)) + R(i) * f**2 / 4) / & + sqrt(beta * dP * delta * edelta / & + rho + (R(i) * f / 2)**2) + & + (sqrt(beta * dP * delta * edelta / & + rho + (R(i) * f / 2)**2)) / R(i) + end if + + Z(i) = sign(Z(i), f) + end do + !$OMP END PARALLEL DO + +end subroutine fhollandvort + + +subroutine fcomb(V, Z, R, d2Vm, dVm, rMax, vMax, beta, dP, rho, f, n) + !$ use omp_lib + doubleprecision, intent(in) :: d2Vm, dVm, rMax, beta, dP, rho, f, vMax + integer, intent(in) :: n + doubleprecision, intent(in), dimension(n) :: R + doubleprecision, intent(inout), dimension(n) :: V, Z + doubleprecision :: aa, bb, cc, delta, edelta + logical :: icore + + aa = (d2Vm / 2 - (dVm - vMax / rMax) / rMax) / rMax + bb = (d2Vm - 6 * aa * rMax) / 2. + cc = dVm - 3 * aa * rMax ** 2 - 2 * bb * rMax + + !$OMP PARALLEL DO shared(Z) + do i = 1, n + delta = (rMax / R(i)) ** beta + edelta = exp(-delta) + icore = (R(i) <= rMax) + if (icore) then + Z(i) = R(i) * (R(i) * 4 * aa + 3 * bb) + 2 * cc + V(i) = (R(i) * (R(i) * (R(i) * aa + bb) + cc)) + else + V(i) = sqrt((dP * beta / rho) * delta * edelta + (R(i) * f / 2.) ** 2) - R(i) * abs(f) / 2. + Z(i) = abs(f) + & + (beta**2 * dP * (delta**2) * edelta / & + (2 * rho * R(i)) - beta**2 * dP * delta * edelta / & + (2 * rho * R(i)) + R(i) * f**2 / 4) / & + sqrt(beta * dP * delta * edelta / & + rho + (R(i) * f / 2)**2) + & + (sqrt(beta * dP * delta * edelta / & + rho + (R(i) * f / 2)**2)) / R(i) + end if + + Z(i) = sign(Z(i), f) + end do + !$OMP END PARALLEL DO + +end subroutine fcomb \ No newline at end of file diff --git a/wind/vmax.py b/wind/vmax.py index e09412c9..9a0ffdff 100644 --- a/wind/vmax.py +++ b/wind/vmax.py @@ -51,7 +51,7 @@ """ -from scipy import sqrt, exp, power +import numpy as np import Utilities.metutils as metutils @@ -98,20 +98,20 @@ def vmax(pCentre, pEnv, type="holland", beta=1.3, rho=1.15): # Primary Hurricane Vortex. Part I: Observations and # Evaluation of the Holland (1980) Model. # Mon. Wea. Rev., 132, 3033-3048 - vMax = 0.6252*sqrt(dP) + vMax = 0.6252*np.sqrt(dP) elif type == "holland": # Holland (1980), An Analytic Model of the Wind and Pressure # Profiles in Hurricanes. Mon. Wea. Rev, 108, 1212-1218 # Density of air is assumed to be 1.15 kg/m^3. # beta is assumed to be 1.3. Other values can be specified. # Gradient level wind (assumed maximum). - vMax = sqrt(beta*dP/(exp(1)*rho)) + vMax = np.sqrt(beta*dP/(np.exp(1)*rho)) elif type == "atkinson": # Atkinson and Holliday (1977), Tropical Cyclone Minimum Sea # Level Pressure / Maximum Sustained Wind Relationship for # the Western North Pacific. Mon. Wea. Rev., 105, 421-427 # Maximum 10m, 1-minute wind speed. Uses pEnv as 1010 hPa - vMax = 3.04*power(1010 - metutils.convert(pCentre, "Pa", "hPa"), 0.644) + vMax = 3.04*np.power(1010 - metutils.convert(pCentre, "Pa", "hPa"), 0.644) else: raise NotImplementedError("Vmax type " + type + " not implemented") return vMax @@ -130,7 +130,7 @@ def pDiff(vMax, pEnv, vMaxType="holland", beta=1.3, rho=1.15): if vMaxType == "willoughby": dP = (vMax/0.6252)**2 elif vMaxType == "holland": - dP = rho*exp(1)*(vMax**2)/beta + dP = rho*np.exp(1)*(vMax**2)/beta elif vMaxType == "atkinson": dP = (vMax/3.04)**(1/0.644) dP = metutils.convert(dP, "hPa", "Pa") diff --git a/wind/windmodels.py b/wind/windmodels.py index b463fccc..e97525de 100644 --- a/wind/windmodels.py +++ b/wind/windmodels.py @@ -39,10 +39,16 @@ from math import exp, sqrt import Utilities.metutils as metutils import logging +import warnings logging.getLogger('matplotlib').setLevel(logging.WARNING) log = logging.getLogger(__name__) log.addHandler(logging.NullHandler()) +try: + from . import fwind +except ImportError: + warnings.warn("Compiled wind models not found - defaulting to slower python wind models") + class WindSpeedModel(object): @@ -371,20 +377,30 @@ def velocity(self, R): d2Vm = self.secondDerivative() dVm = self.firstDerivative() - aa = ((d2Vm / 2. - (dVm - self.vMax / self.rMax) / - self.rMax) / self.rMax) - bb = (d2Vm - 6 * aa * self.rMax) / 2. - cc = dVm - 3 * aa * self.rMax ** 2 - 2 * bb * self.rMax - delta = (self.rMax / R) ** self.beta - edelta = np.exp(-delta) - V = (np.sqrt((self.dP * self.beta / self.rho) * - delta * edelta + (R * self.f / 2.) ** 2) - - R * np.abs(self.f) / 2.) + try: + from .fwind import fhollandvel + V = np.empty_like(R) + fhollandvel( + V.ravel(), R.ravel(), d2Vm, dVm, self.rMax, + self.vMax, self.beta, self.dP, self.rho, self.f, V.size + ) + + except ImportError: + aa = ((d2Vm / 2. - (dVm - self.vMax / self.rMax) / + self.rMax) / self.rMax) + bb = (d2Vm - 6 * aa * self.rMax) / 2. + cc = dVm - 3 * aa * self.rMax ** 2 - 2 * bb * self.rMax + delta = (self.rMax / R) ** self.beta + edelta = np.exp(-delta) + + V = (np.sqrt((self.dP * self.beta / self.rho) * + delta * edelta + (R * self.f / 2.) ** 2) - + R * np.abs(self.f) / 2.) - icore = np.where(R <= self.rMax) - V[icore] = (R[icore] * (R[icore] * (R[icore] * aa + bb) + cc)) - V = np.sign(self.f) * V + icore = np.where(R <= self.rMax) + V[icore] = (R[icore] * (R[icore] * (R[icore] * aa + bb) + cc)) + V = np.sign(self.f) * V return V def vorticity(self, R): @@ -399,31 +415,42 @@ def vorticity(self, R): :rtype: :class:`numpy.ndarray` """ - - beta = self.beta - delta = (self.rMax / R) ** beta - edelta = np.exp(-delta) - - Z = np.abs(self.f) + \ - (beta**2 * self.dP * (delta**2) * edelta / - (2 * self.rho * R) - beta**2 * self.dP * delta * edelta / - (2 * self.rho * R) + R * self.f**2 / 4) / \ - np.sqrt(beta * self.dP * delta * edelta / - self.rho + (R * self.f / 2)**2) + \ - (np.sqrt(beta * self.dP * delta * edelta / - self.rho + (R * self.f / 2)**2)) / R - # Calculate first and second derivatives at R = Rmax: + d2Vm = self.secondDerivative() dVm = self.firstDerivative() - aa = ((d2Vm / 2 - (dVm - self.vMax / - self.rMax) / self.rMax) / self.rMax) - bb = (d2Vm - 6 * aa * self.rMax) / 2 - cc = dVm - 3 * aa * self.rMax ** 2 - 2 * bb * self.rMax - icore = np.where(R <= self.rMax) - Z[icore] = R[icore] * (R[icore] * 4 * aa + 3 * bb) + 2 * cc - Z = np.sign(self.f) * Z + try: + from .fwind import fhollandvort + Z = np.empty_like(R) + fhollandvort( + Z.ravel(), R.ravel(), d2Vm, dVm, self.rMax, self.vMax, + self.beta, self.dP, self.rho, self.f, Z.size + ) + + except ImportError: + beta = self.beta + delta = (self.rMax / R) ** beta + edelta = np.exp(-delta) + + Z = np.abs(self.f) + \ + (beta**2 * self.dP * (delta**2) * edelta / + (2 * self.rho * R) - beta**2 * self.dP * delta * edelta / + (2 * self.rho * R) + R * self.f**2 / 4) / \ + np.sqrt(beta * self.dP * delta * edelta / + self.rho + (R * self.f / 2)**2) + \ + (np.sqrt(beta * self.dP * delta * edelta / + self.rho + (R * self.f / 2)**2)) / R + + aa = ((d2Vm / 2 - (dVm - self.vMax / + self.rMax) / self.rMax) / self.rMax) + bb = (d2Vm - 6 * aa * self.rMax) / 2 + cc = dVm - 3 * aa * self.rMax ** 2 - 2 * bb * self.rMax + + icore = np.where(R <= self.rMax) + Z[icore] = R[icore] * (R[icore] * 4 * aa + 3 * bb) + 2 * cc + Z = np.sign(self.f) * Z + return Z @@ -1037,23 +1064,36 @@ def field(self, R, lam, vFm, thetaFm, thetaMax=0.): """ :param R: Distance from the storm centre to the grid (km). :type R: :class:`numpy.ndarray` - :param lam: Direction (geographic bearing, positive clockwise) + :param lam: Direction (0=east, radians, positive anti-clockwise) from storm centre to the grid. :type lam: :class:`numpy.ndarray` :param float vFm: Foward speed of the storm (m/s). - :param float thetaFm: Forward direction of the storm (geographic - bearing, positive clockwise, radians). + :param float thetaFm: Forward direction of the storm (0=east, radians, + positive anti-clockwise). :param float thetaMax: Bearing of the location of the maximum wind speed, relative to the direction of motion. """ - V = self.velocity(R) - Z = self.vorticity(R) K = 50. # Diffusivity Cd = 0.002 # Constant drag coefficient - Vm = np.abs(V).max() + Vm = self.profile.vMax + + try: + from .fwind import fkerpert + Ux, Vy = np.empty_like(R), np.empty_like(R) + n = Ux.size + fkerpert( + R.ravel(), lam.ravel(), self.f, self.rMax, vFm, thetaFm, + Vm, Ux.ravel(), Vy.ravel(), n + ) + return Ux, Vy + except ImportError: + pass + + V = self.velocity(R) + Z = self.vorticity(R) if (vFm > 0) and (Vm/vFm < 5.): Umod = vFm * np.abs(1.25*(1. - (vFm/Vm))) else: @@ -1090,8 +1130,8 @@ def field(self, R, lam, vFm, thetaFm, thetaMax=0.): Am[ind] = AmIII[ind] # First asymmetric surface component - ums = (Am * np.exp(-i * lam * np.sign(self.f))).real * albe - vms = (Am * np.exp(-i * lam * np.sign(self.f))).imag * np.sign(self.f) + ums = (Am * np.exp(-i * (lam - thetaFm) * np.sign(self.f))).real * albe + vms = (Am * np.exp(-i * (lam - thetaFm) * np.sign(self.f))).imag * np.sign(self.f) Ap = -(eta * (1 - 2 * albe + (1 + i) * (1 - albe) * psi) * Vt) / \ (albe * ((2 + 2 * i) * (1 + eta * psi) + 3 * eta + 3 * i * psi)) @@ -1101,8 +1141,8 @@ def field(self, R, lam, vFm, thetaFm, thetaMax=0.): Ap[ind] = ApIII[ind] # Second asymmetric surface component - ups = (Ap * np.exp(i * lam * np.sign(self.f))).real * albe - vps = (Ap * np.exp(i * lam * np.sign(self.f))).imag * np.sign(self.f) + ups = (Ap * np.exp(i * (lam - thetaFm) * np.sign(self.f))).real * albe + vps = (Ap * np.exp(i * (lam - thetaFm) * np.sign(self.f))).imag * np.sign(self.f) # Total surface wind in (moving coordinate system) us = u0s + ups + ums @@ -1142,9 +1182,9 @@ def profileParams(name): """ List of additional parameters required for a wind profile model. """ - from inspect import getargspec - std = getargspec(WindProfileModel.__init__)[0] - new = getargspec(profile(name).__init__)[0] + from inspect import getfullargspec + std = getfullargspec(WindProfileModel.__init__)[0] + new = getfullargspec(profile(name).__init__)[0] params = [p for p in new if p not in std] return params @@ -1161,9 +1201,9 @@ def fieldParams(name): """ List of additional parameters required for a wind field model. """ - from inspect import getargspec - std = getargspec(WindFieldModel.__init__)[0] - new = getargspec(field(name).__init__)[0] + from inspect import getfullargspec + std = getfullargspec(WindFieldModel.__init__)[0] + new = getfullargspec(field(name).__init__)[0] params = [p for p in new if p not in std] return params diff --git a/wind/writer.py b/wind/writer.py index 405d6888..ca7e00ca 100644 --- a/wind/writer.py +++ b/wind/writer.py @@ -6,6 +6,7 @@ import netCDF4 import affine import numpy as np +import xarray as xr class WriteFoliationCallback(object): @@ -70,7 +71,7 @@ def series(start, stop, inc=resolution): yoff=gridLimit['yMin']-margin) * affine.Affine.scale(resolution)) - self.ds = root = netCDF4.Dataset(filename, mode='w') + self.ds = root = netCDF4.Dataset(filename, mode='w') #, memory=True, persist=True) root.description = "Simulated Windfield Timeseries" # declare shapes @@ -127,3 +128,62 @@ def __call__(self, time, gust, Ux, Uy, P, lon, lat): # save (flush) layer to file self.ds.sync() + + +class WriteMemoryCallback(object): + """ + + + """ + + def __init__(self, filename, gridLimit, resolution, margin=0, wraps=None, **kwargs): + """ + :param str filename: netCDF file for saving data to + :param dict gridLimit: simulation domain extent + :param float resolution: grid resolution in decimal degrees + :param float margin: spatial extent over which the wind field is + calculated in decimal degrees + :param int maxchunk: chunking size (for use in netCDF variable + creation) + :param func wraps: Optional callback function (e.g. for + timeseries extraction) + + """ + + logging.debug( + "Preparing to record windfield evolution to {}".format(filename)) + + self.callback = wraps + self.filename = filename + + def setup(self, coords): + self.tUU = xr.DataArray(dims=["time", "lat", "lon"], coords=coords) + self.tVV = xr.DataArray(dims=["time", "lat", "lon"], coords=coords) + self.tP = xr.DataArray(dims=["time", "lat", "lon"], coords=coords) + self.tgust = xr.DataArray(dims=["time", "lat", "lon"], coords=coords) + + def __call__(self, time, gust, Ux, Uy, P, imin, imax, jmin, jmax): + """Save wind field layer for latest time step""" + # if self.callback: + # self.callback(time, gust, Ux, Uy, P, lon, lat) + self.tUU.sel(time=time).data[imin:imax, jmin:jmax] = Ux + self.tVV.sel(time=time).data[imin:imax, jmin:jmax] = Uy + self.tP.sel(time=time).data[imin:imax, jmin:jmax] = P + self.tgust.sel(time=time).data[imin:imax, jmin:jmax] = gust + + def close(self): + ds = xr.Dataset( + dict( + gust_speed=self.tgust, + velocity_east=self.tUU, + velocity_north=self.tVV, + pressure=self.tP, + ) + ) + + ds.to_netcdf(self.filename) + del ds + del self.tgust + del self.tUU + del self.tVV + del self.tP \ No newline at end of file