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 924fd2b1..77ca7dea 100644 --- a/Evaluate/interpolateTracks.py +++ b/Evaluate/interpolateTracks.py @@ -8,7 +8,10 @@ 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 +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 @@ -101,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] @@ -140,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), @@ -262,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): """ @@ -299,6 +303,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) @@ -317,3 +324,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/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 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/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 + + 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/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/wind/__init__.py b/wind/__init__.py index b6d24b33..437d7951 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], @@ -201,7 +202,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 +270,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 +305,17 @@ 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. @@ -436,11 +453,17 @@ def calculateExtremesFromTrack(self, track, callback=None): output = pjoin(self.windfieldPath, 'evolution.{0:03d}-{1:05d}.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) @@ -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): """ @@ -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` @@ -764,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 @@ -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/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 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