From b97be82b9fb03db173c50cb390b3bf9b6f6f1a80 Mon Sep 17 00:00:00 2001 From: Kieran Ricardo Date: Wed, 16 Mar 2022 14:50:07 +1100 Subject: [PATCH 01/16] run tcrm from a pycxml file --- Evaluate/evaluate.py | 24 ++++++++++++++++++++++-- Evaluate/interpolateTracks.py | 5 +++++ 2 files changed, 27 insertions(+), 2 deletions(-) diff --git a/Evaluate/evaluate.py b/Evaluate/evaluate.py index c306c8c1..9f3cf1bc 100644 --- a/Evaluate/evaluate.py +++ b/Evaluate/evaluate.py @@ -32,7 +32,7 @@ from matplotlib import pyplot, cm from matplotlib.dates import date2num -from mpl_toolkits.basemap import Basemap +#from mpl_toolkits.basemap import Basemap from scipy.stats import scoreatpercentile as percentile from datetime import datetime @@ -45,7 +45,7 @@ from Utilities.maputils import bearing2theta import Utilities.Intersections as Int -import Utilities.colours as colours +#ßimport Utilities.colours as colours import warnings warnings.filterwarnings("ignore", category=RuntimeWarning) @@ -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 924fd2b1..280d78d8 100644 --- a/Evaluate/interpolateTracks.py +++ b/Evaluate/interpolateTracks.py @@ -9,6 +9,8 @@ from Utilities.maputils import latLon2Azi from Utilities.loadData import loadTrackFile, maxWindSpeed from Utilities.track import Track, ncSaveTracks +from .evaluate import bom2tcrm +from pycxml.pycxml import loadfile LOG = logging.getLogger(__name__) LOG.addHandler(logging.NullHandler()) @@ -299,6 +301,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(fp) + tracks = [bom2tcrm(df, i) for i, df in enumerate(dfs)] else: tracks = loadTrackFile(configFile, trackFile, source) From cf3d91af56f71b5a4a996283cbe74f894a4b3837 Mon Sep 17 00:00:00 2001 From: Kieran Ricardo Date: Wed, 16 Mar 2022 14:50:24 +1100 Subject: [PATCH 02/16] run tcrm from a pycxml file --- Evaluate/evaluate.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Evaluate/evaluate.py b/Evaluate/evaluate.py index 9f3cf1bc..9b8b92c6 100644 --- a/Evaluate/evaluate.py +++ b/Evaluate/evaluate.py @@ -32,7 +32,7 @@ from matplotlib import pyplot, cm from matplotlib.dates import date2num -#from mpl_toolkits.basemap import Basemap +from mpl_toolkits.basemap import Basemap from scipy.stats import scoreatpercentile as percentile from datetime import datetime @@ -45,7 +45,7 @@ from Utilities.maputils import bearing2theta import Utilities.Intersections as Int -#ßimport Utilities.colours as colours +import Utilities.colours as colours import warnings warnings.filterwarnings("ignore", category=RuntimeWarning) From b46f9805854081d9ad9eeefba944f6dc69e1fce5 Mon Sep 17 00:00:00 2001 From: Kieran Ricardo Date: Wed, 16 Mar 2022 17:18:22 +1100 Subject: [PATCH 03/16] fixes --- Evaluate/interpolateTracks.py | 23 +++++++++++++++++++++-- Utilities/track.py | 11 +++++++++-- example/yasi.ini | 3 ++- 3 files changed, 32 insertions(+), 5 deletions(-) diff --git a/Evaluate/interpolateTracks.py b/Evaluate/interpolateTracks.py index 280d78d8..cec2878d 100644 --- a/Evaluate/interpolateTracks.py +++ b/Evaluate/interpolateTracks.py @@ -9,8 +9,8 @@ from Utilities.maputils import latLon2Azi from Utilities.loadData import loadTrackFile, maxWindSpeed from Utilities.track import Track, ncSaveTracks -from .evaluate import bom2tcrm from pycxml.pycxml import loadfile +import pandas as pd LOG = logging.getLogger(__name__) LOG.addHandler(logging.NullHandler()) @@ -302,7 +302,7 @@ def parseTracks(configFile, trackFile, source, delta, outputFile=None, from Utilities.track import ncReadTrackData tracks = ncReadTrackData(trackFile) elif trackFile.endswith("xml"): - dfs = loadfile(fp) + dfs = loadfile(trackFile) tracks = [bom2tcrm(df, i) for i, df in enumerate(dfs)] else: tracks = loadTrackFile(configFile, trackFile, source) @@ -322,3 +322,22 @@ def parseTracks(configFile, trackFile, source, delta, outputFile=None, 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/Utilities/track.py b/Utilities/track.py index 7f9bfa07..84b6b98f 100644 --- a/Utilities/track.py +++ b/Utilities/track.py @@ -81,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. diff --git a/example/yasi.ini b/example/yasi.ini index f8fcb732..49677c3f 100644 --- a/example/yasi.ini +++ b/example/yasi.ini @@ -1,5 +1,6 @@ +#./example/yasi.csv [DataProcess] -InputFile = ./example/yasi.csv +InputFile = /scratch/w85/cxml/2011013100.xml Source = BDECK FilterSeasons = False From 0479d314fa0467d35bbceae84387c217e11be4e8 Mon Sep 17 00:00:00 2001 From: Kieran Ricardo Date: Thu, 17 Mar 2022 10:54:11 +1100 Subject: [PATCH 04/16] handle negative central pressure deficits --- wind/__init__.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/wind/__init__.py b/wind/__init__.py index b6d24b33..43240abc 100644 --- a/wind/__init__.py +++ b/wind/__init__.py @@ -201,7 +201,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) From 384c07247e71cc8307045331d0dfde810bbcc16f Mon Sep 17 00:00:00 2001 From: Kieran Ricardo Date: Thu, 17 Mar 2022 12:52:00 +1100 Subject: [PATCH 05/16] revert example back --- example/port_hedland.ini | 107 --------------------------------------- example/yasi.ini | 3 +- 2 files changed, 1 insertion(+), 109 deletions(-) diff --git a/example/port_hedland.ini b/example/port_hedland.ini index 77f8beda..e69de29b 100644 --- a/example/port_hedland.ini +++ b/example/port_hedland.ini @@ -1,107 +0,0 @@ -[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 diff --git a/example/yasi.ini b/example/yasi.ini index 49677c3f..f8fcb732 100644 --- a/example/yasi.ini +++ b/example/yasi.ini @@ -1,6 +1,5 @@ -#./example/yasi.csv [DataProcess] -InputFile = /scratch/w85/cxml/2011013100.xml +InputFile = ./example/yasi.csv Source = BDECK FilterSeasons = False From 432f4c2a451378be936b82bdf7d14b44e0bd0c6e Mon Sep 17 00:00:00 2001 From: Kieran Ricardo Date: Fri, 18 Mar 2022 15:22:34 +1100 Subject: [PATCH 06/16] paralell tracks --- wind/__init__.py | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/wind/__init__.py b/wind/__init__.py index 43240abc..eff2d9ab 100644 --- a/wind/__init__.py +++ b/wind/__init__.py @@ -784,14 +784,23 @@ 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 + def loadTracks(trackfile): """ Read tracks from a track .nc file and return a list of :class:`Track` From ffdd8fb5c7f522bc21d846b0833fd62e9be033fc Mon Sep 17 00:00:00 2001 From: Kieran Ricardo Date: Mon, 21 Mar 2022 17:56:51 +1100 Subject: [PATCH 07/16] xarray performance enhancements - 25% faster overall and 20x faster data writing --- Utilities/timeseries.py | 42 ++++++++++++++++++++--------- tcevent.py | 9 ++++--- wind/__init__.py | 41 +++++++++++++++++++++------- wind/writer.py | 60 ++++++++++++++++++++++++++++++++++++++++- 4 files changed, 125 insertions(+), 27 deletions(-) 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/tcevent.py b/tcevent.py index f48962d3..2306c746 100755 --- a/tcevent.py +++ b/tcevent.py @@ -188,10 +188,10 @@ def main(configFile): 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') + # interpTrack = interpolateTracks.parseTracks(configFile, trackFile, + # source, delta, + # outputTrackFile, + # interpolation_type='akima') showProgressBar = config.get('Logging', 'ProgressBar') @@ -218,6 +218,7 @@ def startup(): function. """ + print("starting up") parser = argparse.ArgumentParser() parser.add_argument('-c', '--config_file', help='Path to configuration file') diff --git a/wind/__init__.py b/wind/__init__.py index eff2d9ab..b4370844 100644 --- a/wind/__init__.py +++ b/wind/__init__.py @@ -268,6 +268,13 @@ def regionalExtremes(self, gridLimit, timeStepCallback=None): (self.track.Latitude <= yMax))[0] nsteps = len(self.track.TimeElapsed) + + timesInRegion = timesInRegion[:5] + + 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 " @@ -291,21 +298,25 @@ def regionalExtremes(self, gridLimit, timeStepCallback=None): gridMargin) / gridStep) + 1 # Calculate the local wind speeds and pressure at time i - Ux, Vy, P = self.localWindField(i) + Ux, Vy, P = self.localWindField(i) # Calculate the local wind gust and bearing Ux *= self.gustFactor Vy *= self.gustFactor localGust = np.sqrt(Ux ** 2 + Vy ** 2) localBearing = ((np.arctan2(-Ux, -Vy)) * 180. / np.pi) - + t2 = time.time() # 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] @@ -324,6 +335,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. @@ -440,10 +454,11 @@ def calculateExtremesFromTrack(self, track, callback=None): output = pjoin(self.windfieldPath, 'evolution.{0:03d}-{1:05d}.nc'.format( *track.trackId)) - callback = writer.WriteFoliationCallback(output, + callback = writer.WriteMemoryCallback(output, self.gridLimit, self.resolution, self.margin, + maxchunk=2048, wraps=callback) return track, wt.regionalExtremes(self.gridLimit, callback) @@ -531,6 +546,7 @@ def dumpGustsFromTracks(self, trackiter, windfieldPath, if self.multipliers is not None: self.calcLocalWindfield(results) + def plotGustToFile(self, result, filename): """ @@ -704,8 +720,7 @@ def dumpGustsFromTrackfiles(self, trackfiles, windfieldPath, """ tracks = loadTracksFromFiles(sorted(trackfiles), self.gridLimit, self.margin) - self.dumpGustsFromTracks(tracks, windfieldPath, - timeStepCallback=timeStepCallback) + self.dumpGustsFromTracks(tracks, windfieldPath, timeStepCallback=timeStepCallback) def calcLocalWindfield(self, results): """ @@ -847,6 +862,7 @@ def balanced(iterable): scattering. """ P, p = MPI.COMM_WORLD.size, MPI.COMM_WORLD.rank + iterable = itertools.islice(iterable, 0, 2) return itertools.islice(iterable, p, None, P) @@ -894,6 +910,7 @@ def run(configFile, callback=None): else: timestepCallback = None + print("Timestep callback:", timestepCallback) multipliers = None if config.has_option('Input', 'Multipliers'): multipliers = config.get('Input', 'Multipliers') @@ -945,3 +962,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 \ No newline at end of file diff --git a/wind/writer.py b/wind/writer.py index 405d6888..750241ce 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,60 @@ 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) + + + From c157dcbee8072976823b54916e5239fd328a30ff Mon Sep 17 00:00:00 2001 From: Kieran Ricardo Date: Wed, 23 Mar 2022 10:21:22 +1100 Subject: [PATCH 08/16] wip --- wind/__init__.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/wind/__init__.py b/wind/__init__.py index b4370844..2dac7893 100644 --- a/wind/__init__.py +++ b/wind/__init__.py @@ -268,9 +268,7 @@ def regionalExtremes(self, gridLimit, timeStepCallback=None): (self.track.Latitude <= yMax))[0] nsteps = len(self.track.TimeElapsed) - - timesInRegion = timesInRegion[:5] - + coords = coords=dict(lon=lonGrid, lat=latGrid, time=timesInRegion) if timeStepCallback is not None: @@ -306,7 +304,6 @@ def regionalExtremes(self, gridLimit, timeStepCallback=None): localGust = np.sqrt(Ux ** 2 + Vy ** 2) localBearing = ((np.arctan2(-Ux, -Vy)) * 180. / np.pi) - t2 = time.time() # Handover this time step to a callback if required if timeStepCallback is not None: @@ -460,6 +457,11 @@ def calculateExtremesFromTrack(self, track, callback=None): 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) @@ -862,7 +864,7 @@ def balanced(iterable): scattering. """ P, p = MPI.COMM_WORLD.size, MPI.COMM_WORLD.rank - iterable = itertools.islice(iterable, 0, 2) + iterable = itertools.islice(iterable, 0, 100) return itertools.islice(iterable, p, None, P) From 5e6c1bc5566cbd2146794795f33087de3e958ded Mon Sep 17 00:00:00 2001 From: Kieran Ricardo Date: Wed, 23 Mar 2022 14:21:50 +1100 Subject: [PATCH 09/16] wip --- PressureInterface/pressureProfile.f90 | 2 +- wind/__init__.py | 8 ++++++++ 2 files changed, 9 insertions(+), 1 deletion(-) diff --git a/PressureInterface/pressureProfile.f90 b/PressureInterface/pressureProfile.f90 index 429b53ad..bd50ca79 100644 --- a/PressureInterface/pressureProfile.f90 +++ b/PressureInterface/pressureProfile.f90 @@ -7,7 +7,7 @@ subroutine fhollandpressure(P, R, rMax, pc, dP, beta, n) !$OMP PARALLEL DO shared(P) do i = 1, n - P(i) = pCentre + dP * exp(-(rMax / R) ** beta) + P(i) = pCentre + dP * exp(-(rMax / R(i)) ** beta) end do !$OMP END PARALLEL DO diff --git a/wind/__init__.py b/wind/__init__.py index 2dac7893..72e94017 100644 --- a/wind/__init__.py +++ b/wind/__init__.py @@ -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], @@ -298,10 +299,17 @@ def regionalExtremes(self, gridLimit, timeStepCallback=None): # Calculate the local wind speeds and pressure at time i Ux, Vy, P = self.localWindField(i) + from cProfile import runctx + runctx("_ = self.localWindField(i)", locals(), globals(), sort='cumtime') + print("####\n" * 10) + runctx("_ = self.localWindField(i)", locals(), globals(), sort='tottime') + print("####\n" * 10) + raise ZeroDivisionError # Calculate the local wind gust and bearing Ux *= self.gustFactor Vy *= self.gustFactor + localGust = np.sqrt(Ux ** 2 + Vy ** 2) localBearing = ((np.arctan2(-Ux, -Vy)) * 180. / np.pi) # Handover this time step to a callback if required From 0adba4977960330d86f1cde0fc36ae18b4cf6533 Mon Sep 17 00:00:00 2001 From: Kieran Ricardo Date: Wed, 23 Mar 2022 16:19:37 +1100 Subject: [PATCH 10/16] remove debug --- wind/__init__.py | 10 ++-------- 1 file changed, 2 insertions(+), 8 deletions(-) diff --git a/wind/__init__.py b/wind/__init__.py index 72e94017..1cb271f4 100644 --- a/wind/__init__.py +++ b/wind/__init__.py @@ -297,14 +297,9 @@ def regionalExtremes(self, gridLimit, timeStepCallback=None): gridMargin) / gridStep) + 1 # Calculate the local wind speeds and pressure at time i - Ux, Vy, P = self.localWindField(i) - from cProfile import runctx - runctx("_ = self.localWindField(i)", locals(), globals(), sort='cumtime') - print("####\n" * 10) - runctx("_ = self.localWindField(i)", locals(), globals(), sort='tottime') - print("####\n" * 10) - raise ZeroDivisionError + + # Calculate the local wind gust and bearing Ux *= self.gustFactor Vy *= self.gustFactor @@ -920,7 +915,6 @@ def run(configFile, callback=None): else: timestepCallback = None - print("Timestep callback:", timestepCallback) multipliers = None if config.has_option('Input', 'Multipliers'): multipliers = config.get('Input', 'Multipliers') From 6c0eaac34fe1bb906eeb80fde01f69e35eba860e Mon Sep 17 00:00:00 2001 From: Kieran Ricardo Date: Thu, 24 Mar 2022 09:55:31 +1100 Subject: [PATCH 11/16] cleanup memory --- wind/writer.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/wind/writer.py b/wind/writer.py index 750241ce..ca7e00ca 100644 --- a/wind/writer.py +++ b/wind/writer.py @@ -182,6 +182,8 @@ def close(self): ) 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 From c6185abf0c5c75f68f9cc0acf4f86b4f9b0bb48a Mon Sep 17 00:00:00 2001 From: wcarthur Date: Fri, 20 Sep 2024 16:03:49 +1000 Subject: [PATCH 12/16] Remove double definition of iterable --- wind/__init__.py | 1 - 1 file changed, 1 deletion(-) diff --git a/wind/__init__.py b/wind/__init__.py index 1cb271f4..495f027b 100644 --- a/wind/__init__.py +++ b/wind/__init__.py @@ -867,7 +867,6 @@ def balanced(iterable): scattering. """ P, p = MPI.COMM_WORLD.size, MPI.COMM_WORLD.rank - iterable = itertools.islice(iterable, 0, 100) return itertools.islice(iterable, p, None, P) From dcbe470dd6eb1c4c2ae66641737cf5e4f3114aa5 Mon Sep 17 00:00:00 2001 From: wcarthur Date: Fri, 20 Sep 2024 16:04:39 +1000 Subject: [PATCH 13/16] Remove old scipy imports --- wind/vmax.py | 1 - 1 file changed, 1 deletion(-) diff --git a/wind/vmax.py b/wind/vmax.py index 7e7eb51d..9a0ffdff 100644 --- a/wind/vmax.py +++ b/wind/vmax.py @@ -51,7 +51,6 @@ """ -from scipy import sqrt, exp, power import numpy as np import Utilities.metutils as metutils From c1b01ee914db0e570fcb20dbca4d8e4526e5b6aa Mon Sep 17 00:00:00 2001 From: wcarthur Date: Fri, 20 Sep 2024 16:06:58 +1000 Subject: [PATCH 14/16] Enable MPI execution --- Evaluate/interpolateTracks.py | 8 +++++--- tcevent.py | 20 ++++++++++++++++---- 2 files changed, 21 insertions(+), 7 deletions(-) diff --git a/Evaluate/interpolateTracks.py b/Evaluate/interpolateTracks.py index cec2878d..77ca7dea 100644 --- a/Evaluate/interpolateTracks.py +++ b/Evaluate/interpolateTracks.py @@ -8,6 +8,7 @@ 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 @@ -72,7 +73,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 @@ -103,7 +104,7 @@ def interpolate(track, delta, interpolation_type=None): 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] @@ -142,7 +143,7 @@ def interpolate(track, delta, interpolation_type=None): # Use the Akima interpolation method: try: import akima - except ImportError: + except ModuleNotFoundError: LOG.exception(("Akima interpolation module unavailable " " - default to scipy.interpolate")) nLon = splev(newtime, splrep(timestep, track.Longitude, s=0), @@ -264,6 +265,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): """ diff --git a/tcevent.py b/tcevent.py index 2306c746..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. @@ -188,10 +193,11 @@ def main(configFile): 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') + interpTrack = interpolateTracks.parseTracks(configFile, trackFile, + source, delta, + outputTrackFile, + interpolation_type='akima') + comm.barrier() showProgressBar = config.get('Logging', 'ProgressBar') @@ -256,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 From 24e954c7453a999d22d5a7af9aaca7c53ca576c8 Mon Sep 17 00:00:00 2001 From: Craig Arthur Date: Tue, 24 Sep 2024 10:38:39 +1000 Subject: [PATCH 15/16] Keep Port Hedland example, formatting changes --- example/port_hedland.ini | 107 +++++++++++++++++++++++++++++++++++++++ wind/__init__.py | 42 +++++++-------- 2 files changed, 129 insertions(+), 20 deletions(-) diff --git a/example/port_hedland.ini b/example/port_hedland.ini index e69de29b..acfff15e 100644 --- a/example/port_hedland.ini +++ b/example/port_hedland.ini @@ -0,0 +1,107 @@ +[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/wind/__init__.py b/wind/__init__.py index 495f027b..437d7951 100644 --- a/wind/__init__.py +++ b/wind/__init__.py @@ -269,9 +269,9 @@ def regionalExtremes(self, gridLimit, timeStepCallback=None): (self.track.Latitude <= yMax))[0] nsteps = len(self.track.TimeElapsed) - - coords = coords=dict(lon=lonGrid, lat=latGrid, time=timesInRegion) - + + coords = coords = dict(lon=lonGrid, lat=latGrid, time=timesInRegion) + if timeStepCallback is not None: timeStepCallback.setup(coords) @@ -299,19 +299,18 @@ def regionalExtremes(self, gridLimit, timeStepCallback=None): # Calculate the local wind speeds and pressure at time i Ux, Vy, P = self.localWindField(i) - # Calculate the local wind gust and bearing Ux *= self.gustFactor Vy *= self.gustFactor - 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(i, localGust, Ux, Vy, P, imin, imax, jmin, jmax) - + + timeStepCallback(i, localGust, Ux, Vy, P, + imin, imax, jmin, jmax) + # timeStepCallback(self.track.Datetime[i], # localGust, Ux, Vy, P, # lonGrid[imin:imax] / 100., @@ -455,11 +454,11 @@ def calculateExtremesFromTrack(self, track, callback=None): 'evolution.{0:03d}-{1:05d}.nc'.format( *track.trackId)) callback = writer.WriteMemoryCallback(output, - self.gridLimit, - self.resolution, - self.margin, - maxchunk=2048, - wraps=callback) + self.gridLimit, + self.resolution, + self.margin, + maxchunk=2048, + wraps=callback) # callback = writer.WriteFoliationCallback(output, # self.gridLimit, # self.resolution, @@ -551,7 +550,6 @@ def dumpGustsFromTracks(self, trackiter, windfieldPath, if self.multipliers is not None: self.calcLocalWindfield(results) - def plotGustToFile(self, result, filename): """ @@ -724,8 +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): """ @@ -752,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 @@ -768,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` @@ -788,6 +790,7 @@ def filterTracks(tracks, gridLimit, margin): return validTracks + def loadTracksFromFiles(trackfiles, gridLimit, margin): """ Generator that yields :class:`Track` objects from a list of track @@ -820,7 +823,6 @@ def loadTracksFromFiles(trackfiles, gridLimit, margin): yield track - def loadTracks(trackfile): """ Read tracks from a track .nc file and return a list of :class:`Track` @@ -968,6 +970,6 @@ def run(configFile, callback=None): # total time to process one track: 218s -## 215s regionalExtremes -### 160s in localWindfield -### 15s in polarGridAroundEye \ No newline at end of file +# 215s regionalExtremes +# 160s in localWindfield +# 15s in polarGridAroundEye From d15b6483f95a964c30ff092bccb7e9ad7e067804 Mon Sep 17 00:00:00 2001 From: wcarthur Date: Thu, 17 Oct 2024 13:28:58 +1100 Subject: [PATCH 16/16] Add ensemble example to README --- docs/examples.rst | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/docs/examples.rst b/docs/examples.rst index 900fbc5a..bb032b80 100644 --- a/docs/examples.rst +++ b/docs/examples.rst @@ -69,3 +69,22 @@ the Bureau of Meteorology. See the :ref:`scenariomodelling` section for more details. +Ensemble simulation +------------------- + +Some forecasting agencies provide ensemble track forecasts. We can run +`tcevent.py` in parallel to simulate the wind fields from such a forecast. + +An example is provided in `examples/ensemble.ini` + +Requirements:: + +* `pycxml` is a prototype library that reads CycloneXML files, used for +distributing ensemble track forecasts (e.g. by the Australian Bureau of +Meteorology). + +To run:: + + $ mpirun -np 16 python tcevent.py -c example/ensemble.ini + +