diff --git a/.travis.yml b/.travis.yml index 3356e207..e58c08b2 100644 --- a/.travis.yml +++ b/.travis.yml @@ -2,6 +2,7 @@ language: python python: - '3.6' - '3.7' + - '3.8' install: - sudo apt-get update - wget https://repo.continuum.io/miniconda/Miniconda3-latest-Linux-x86_64.sh @@ -15,7 +16,7 @@ install: - conda update -q conda - conda config --add channels conda-forge - conda config --set channel_priority strict - - conda create -q -n test-env python=$TRAVIS_PYTHON_VERSION pip numpy scipy matplotlib basemap shapely nose netcdf4 cftime coverage coveralls pycurl pyproj seaborn simplejson sqlite statsmodels libgdal gdal configparser cartopy affine tqdm xarray gxx_linux-64 pthread-stubs imageio + - conda create -q -n test-env python=$TRAVIS_PYTHON_VERSION pip numpy scipy matplotlib basemap shapely nose netcdf4 cftime coverage coveralls pycurl pyproj seaborn simplejson sqlite statsmodels libgdal gdal configparser cartopy affine tqdm xarray gxx_linux-64 pthread-stubs imageio boto3 botocore - source activate test-env branches: except: diff --git a/Evaluate/landfallIntensity.py b/Evaluate/landfallIntensity.py new file mode 100644 index 00000000..eeab38f4 --- /dev/null +++ b/Evaluate/landfallIntensity.py @@ -0,0 +1,335 @@ +#!/usr/bin/env python +# coding: utf-8 + +# # Landfall rates for simulated tropical cyclones +# +# This notebook looks at the rate of landfalling tropical cyclones simulated in +# the TCHA. "Landfall" is a loose definition, where we have set out a series of +# gates around the coastline, which are set 50 km off the coast, and each gate +# is 200 km wide. +# +# We look at landfall as this is where the impacts of TCs are felt. We will +# compare to historical rates at a later point. +# +# The gates are defined in a vector shapefile as line segments. We'll look at the +# counts, and the mean intensity at landfall to see how well it corresponds to +# historical landfall rates and intensity. +# +# Note however, we'd expect the historical record to be somewhat different to the +# mean values determined here. Arguably, the historical record should be (statistically) +# indistiguishable from the range of scenarios. It should appear to be like any other +# member of the distribution. +# +# The simulations used here are 35-year simulations of TC activity in the Australian region. +# This corresponds to the length of historical record used as the input record +# to the TCHA (1981-2016). The TCs are simulated with a 3-hour timestep. + +import os +import sys +from os import walk +from os.path import join as pjoin +import matplotlib.pyplot as plt +from Utilities import track + +from shapely.geometry import LineString, Point + +import numpy as np +import pandas as pd +import geopandas as gpd + +import logging +import seaborn as sns +sns.set_style('whitegrid') + +LOGGER = logging.getLogger(__name__) + +# Start with reading in the gates into a `GeoDataFrame`, and adding some additional +# attributes. This `GeoDataFrame` will be duplicated for each simulation, then +# aggregated for the summary statistics. + +def loadGateFile(gateFile): + try: + gates = gpd.read_file(gateFile) + except: + LOGGER.error(f"Cannot open {gateFile}") + LOGGER.error("Check the file exists and is a valid shapefile") + sys.exit() + + # Add some extra attributes + gates['sim'] = 0 + gates['count'] = 0 + gates['meanlfintensity'] = np.nan + gates['minlfintensity'] = np.nan + gates['cat1'] = 0 + gates['cat2'] = 0 + gates['cat3'] = 0 + gates['cat4'] = 0 + gates['cat5'] = 0 + + return gates + +# Define a function to read in the track files. The track files are the +# TCRM format files, which are netCDF files, with a heirarchical structure. +# The `tracks.ncReadTrackData` function does this, but we still need to convert +# each track to a `GeoDataFrame` and add a `geometry` attribute that +# represents the line of each time step in the track history. We also assign +# a category attribute to each segment, which at this point in time is simply +# based on categorising the central pressure according to the Bureau of +# Meteorology's TC intensity scale. + +def readTracks(trackFile): + """ + Read a track file and create a `GeoPandas.GeoDataFrame` from the data, with + separate polyline features for each unique track in the file. Also adds a + nominal intensity category based on the central pressure value. + + :param str trackFile: Path to a TCRM-format track file + + :returns: `GeoPandas.GeoDataFrame` of tracks + """ + + LOGGER.debug(f"Loading track data from {trackFile}") + tracks = track.ncReadTrackData(trackFile) + trackgdf = [] + for t in tracks: + segments = [] + for n in range(len(t.data) - 1): + segment = LineString([[t.Longitude[n], t.Latitude[n]],[t.Longitude[n+1], t.Latitude[n+1]]]) + segments.append(segment) + gdf = gpd.GeoDataFrame.from_records(t.data[:-1]) + gdf['geometry'] = segments + gdf['category'] = pd.cut(gdf['CentralPressure'], + bins=[0, 930, 955, 970, 985, 990, 1020], + labels=[5,4,3,2,1,0]) + trackgdf.append(gdf) + trackgdf = pd.concat(trackgdf) + return trackgdf + +def isLeft(line, point): + """ + Test whether a point is to the left of a (directed) line segment. + + :param line: :class:`Shapely.geometry.LineString` of the line feature being tested + :param point: :class:`Shapely.geometry.Point` being tested + + :returns: `True` if the point is to the left of the line segment, `False` otherwise + """ + start = Point(line.coords[0]) + end = Point(line.coords[1]) + + det = (end.x - start.x) * (point.y - start.y) - (end.y - start.y) * (point.x - start.x) + if det > 0: return True + if det <= 0: return False + +def isLandfall(gate, tracks): + """ + Determine the count of tracks crossing a gate segment. + + :param gate: `shapely.Geometry.LineSegment` + :param tracks: `GeoPandas.GeoDataFrame` of tracks, where the `geometry` field is a + `LineSegment` + + """ + LOGGER.debug(f"Determining landfalls for track collection") + + crossings = tracks.crosses(gate.geometry) + landfall = [] + for t in tracks[crossings].itertuples(): + if isLeft(gate.geometry, Point(t.geometry.coords[0])): + landfall.append(True) + else: + landfall.append(False) + + return tracks[crossings][landfall] + + +# This function counts the number of track segments that cross the coastal gates. + +def countCrossings(gates, tracks, sim): + """ + Count the crossing rate of all gates for all tracks in a given simulation. + + :param gates: `GeoDataFrame` containing the landfall gates + :param tracks: `GeoDataFrame` containing `Track` objects + :param int sim: Ordinal simulation number + + """ + gates['sim'] = sim + for i, gate in enumerate(gates.itertuples(index=False)): + ncrossings = 0 + l = isLandfall(gate, tracks) + ncrossings = len(l) + if ncrossings > 0: + gates['count'].iloc[i] = ncrossings + gates['meanlfintensity'].iloc[i] = l['CentralPressure'].mean() + gates['minlfintensity'].iloc[i] = l['CentralPressure'].min() + cathist, bins = np.histogram(l['category'].values, bins=[0,1,2,3,4,5, 6]) + gates['cat1'].iloc[i] = cathist[0] + gates['cat2'].iloc[i] = cathist[1] + gates['cat3'].iloc[i] = cathist[2] + gates['cat4'].iloc[i] = cathist[3] + gates['cat5'].iloc[i] = cathist[4] + else: + gates['count'].iloc[i] = 0 + gates['meanlfintensity'].iloc[i] = np.nan + gates['minlfintensity'].iloc[i] = np.nan + gates['cat1'].iloc[i] = 0 + gates['cat2'].iloc[i] = 0 + gates['cat3'].iloc[i] = 0 + gates['cat4'].iloc[i] = 0 + gates['cat5'].iloc[i] = 0 + + return gates + +# A bunch of helper functions to return quantiles +def q10(x): return x.quantile(0.1) +def q90(x): return x.quantile(0.9) +def q25(x): return x.quantile(0.25) +def q75(x): return x.quantile(0.75) + + + +# Get the list of track files from the input data path (this'll become a +# config option in a scripted version of theis notebook). Just find all +# files in the directory that end with "nc" - assume that you're pointing +# to a directory that only has track files. + +def loadLandfallRates(datapath, gates): + + if not os.path.isdir(datapath): + LOGGER.error(f"{datapath} is not a valid directory") + raise IOError(f"{datapath} is not a valid directory") + filelist = [] + for (dirpath, dirnames, filenames) in walk(datapath): + filelist.extend([fn for fn in filenames if fn.endswith('nc')]) + break + nfiles = len(filelist) + if nfiles==0: + LOGGER.error(f"There are no valid trackfiles in the input path {datapath}") + raise IOError(f"There are no valid trackfiles in the input path {datapath}") + LOGGER.info(f"There are {nfiles} track files in {datapath}") + + + # Now we loop through all the gates and determine the landfall rates for each simulation. + + gatedflist = [] + for sim, f in enumerate(filelist): + q, r = np.divmod(sim*10, len(filelist)) + if r==0: + LOGGER.info(f"Loaded {sim} events ({q*10}%)") + gatedf = gates.copy() + tracks = readTracks(pjoin(datapath,f)) + gatedf = countCrossings(gatedf, tracks, sim) + gatedflist.append(gatedf) + + gatesummary = pd.concat(gatedflist) + LOGGER.info("Grouping data") + + + gs = gatesummary.groupby('gate').agg({'count':['sum',np.nanmean, np.nanstd, 'min', 'max', q10, q90], + 'cat1': ['sum',np.nanmean, np.nanstd], + 'cat2': ['sum',np.nanmean, np.nanstd], + 'cat3': ['sum',np.nanmean, np.nanstd], + 'cat4': ['sum',np.nanmean, np.nanstd], + 'cat5': ['sum',np.nanmean, np.nanstd], + 'meanlfintensity':[np.nanmean, q10, q25, q75, q90], + 'minlfintensity':[np.nanmean,'min', np.nanstd]}, as_index=False) + + gs.columns = ['_'.join(col).strip() for col in gs.columns.values] + gs.reset_index(col_level=1) + gs.columns = gs.columns.get_level_values(0) + + + gatedata = gates[['gate', 'longitude', 'latitude', 'label', 'geometry']].join(gs) + + return gatedata + +def plotLandfallIntensity(df, datapath): + + LOGGER.info("Plotting landfall intensity") + width=0.4 + fig, ax = plt.subplots(3,1,figsize=(12,16),sharex=True) + cat12 = np.add(df['cat1_nanmean'], df['cat2_nanmean']).tolist() + cat123 = np.add(cat12, df['cat3_nanmean']).tolist() + cat1234 = np.add(cat123, df['cat4_nanmean']).tolist() + ax[0].bar(df['gate'], df['cat1_nanmean'], color='b', label="Cat 1") + ax[0].bar(df['gate'], df['cat2_nanmean'], bottom=df['cat1_nanmean'], color='g', label='Cat 2') + ax[0].bar(df['gate'], df['cat3_nanmean'], bottom=cat12, color='y', label='Cat 3') + ax[0].bar(df['gate'], df['cat4_nanmean'], bottom=cat123, color='orange', label='Cat 4') + ax[0].bar(df['gate'], df['cat5_nanmean'], bottom=cat1234, color='r', label='Cat 5') + + ax[0].legend() + ax[0].set_ylabel("Mean number of TCs") + ax[1].plot(df['gate'], df['minlfintensity_nanmean'], + label='Minimum landfall intensity') + ax[1].plot(df['gate'], df['meanlfintensity_nanmean'], color='r', + label='Mean landfall intensity') + ax[1].fill_between(df['gate'], df['meanlfintensity_q10'], + df['meanlfintensity_q90'], color='r', alpha=0.25) + + ax[1].legend(loc=2) + ax[1].set_ylim((900, 1020)) + ax[1].set_ylabel("Pressure (hPa)") + ax[2].plot(df['gate'], df['count_sum']/1000) + ax[2].fill_between(df['gate'], df['count_q90'], df['count_q10'], alpha=0.25) + ax[2].set_xlim((0,48)) + ax[2].set_xticks(np.arange(0,48,2)) + ax[2].set_yticks(np.arange(0,2.1,.2)) + ax[2].set_xticklabels(df['label'][::2], rotation='vertical') + ax[2].set_ylabel("Mean proportion of landfall") + plt.savefig(os.path.join(datapath, "landfall_intensity.png"), bbox_inches='tight') + +def plotIntensityDistribution(df, plotpath): + """ + Plot a distribution of intenstity at landfall + + :param df: `DataFrame` containing the gate data + :param str plotpath: Path to where the plots will be saved/ + """ + LOGGER.info("Plotting landfall intensity distribution") + width=0.4 + fig, ax = plt.subplots(1,1, figsize=(12,6), sharex=True) + cat12 = np.add(df['cat1_nanmean'], df['cat2_nanmean']).tolist() + cat123 = np.add(cat12, df['cat3_nanmean']).tolist() + cat1234 = np.add(cat123, df['cat4_nanmean']).tolist() + ax.bar(df['gate'], df['cat1_nanmean'], color='b', label="Cat 1") + ax.bar(df['gate'], df['cat2_nanmean'], bottom=df['cat1_nanmean'], color='g', label='Cat 2') + ax.bar(df['gate'], df['cat3_nanmean'], bottom=cat12, color='y', label='Cat 3') + ax.bar(df['gate'], df['cat4_nanmean'], bottom=cat123, color='orange', label='Cat 4') + ax.bar(df['gate'], df['cat5_nanmean'], bottom=cat1234, color='r', label='Cat 5') + + ax.legend() + ax.set_ylabel("Number of TCs") + ax.set_xlim((0,48)) + ax.set_xticks(np.arange(0,48,2)) + ax.set_yticks(np.arange(0,1.1,.2)) + ax.set_xticklabels(gatedata['label'][::2], rotation='vertical') + ax.set_ylabel("Mean rate of landfall") + plt.savefig(os.path.join(plotpath, "mean_landfall_rate_intensity.png"), bbox_inches='tight') + + +if __name__ == "__main__": + + logFormat = "%(asctime)s: %(funcName)s: %(message)s" + logging.basicConfig(level='INFO', + format=logFormat, + filename='landfallIntensity.log', + filemode='w', + datefmt="%Y-%m-%d %H:%M:%S") + console = logging.StreamHandler(sys.stdout) + console.setLevel(getattr(logging, 'INFO')) + formatter = logging.Formatter(logFormat, + datefmt='%H:%M:%S', ) + console.setFormatter(formatter) + LOGGER.addHandler(console) + LOGGER.info(f"Started {sys.argv[0]} (pid {os.getpid()})") + gatefile = "/g/data/w85/QFES_SWHA/hazard/input/gates.shp" + gates = loadGateFile(gatefile) + datapath = "/g/data/w85/QFES_SWHA/hazard/output/GROUP2_RCP45_1981-2010/tracks" + plotpath = "/g/data/w85/QFES_SWHA/hazard/output/GROUP2_RCP45_1981-2010/plots" + gatedata = loadLandfallRates(datapath, gates) + plotIntensityDistribution(gatedata, plotpath) + plotLandfallIntensity(gatedata, plotpath) + gatedata_nogeom = pd.DataFrame(gatedata.drop(columns='geometry')) + gatedata_nogeom.to_csv(os.path.join(plotpath, "simulated_landfall_rates.csv"), index=False) + gatedata.to_file(os.path.join(plotpath, "simulated_landfall_rates.shp")) diff --git a/Utilities/tracks2shp.py b/Utilities/tracks2shp.py index 16a3ecc9..ff981e28 100644 --- a/Utilities/tracks2shp.py +++ b/Utilities/tracks2shp.py @@ -21,10 +21,10 @@ OBSFIELD_NAMES = ('Indicator', 'TCID', 'Year', 'Month', 'Day', 'Hour', 'Minute', 'TElapsed', 'Longitude', 'Latitude', 'Speed', 'Bearing', 'Pcentre', - 'MaxWind', 'rMax', 'Penv') -OBSFIELD_TYPES = ('N',)*16 -OBSFIELD_WIDTH = (1, 6, 4, 2, 2, 2, 2, 6, 7, 7, 6, 6, 7, 6, 6, 7) -OBSFIELD_PREC = (0, 0, 0, 0, 0, 0, 0, 2, 2, 2, 1, 1, 1, 1, 1, 1) + 'MaxWind', 'rMax', 'Penv', 'Category') +OBSFIELD_TYPES = ('N',)*17 +OBSFIELD_WIDTH = (1, 6, 4, 2, 2, 2, 2, 6, 7, 7, 6, 6, 7, 6, 6, 7, 1) +OBSFIELD_PREC = (0, 0, 0, 0, 0, 0, 0, 2, 2, 2, 1, 1, 1, 1, 1, 1, 0) OBSFIELDS = [[n, t, w, p] for n, t, w, p in zip(OBSFIELD_NAMES, OBSFIELD_TYPES, @@ -32,10 +32,11 @@ OBSFIELD_PREC)] TCRM_FIELD_NAMES = ('CycloneNumber', 'TimeElapsed', 'Longitude', 'Latitude', - 'Speed', 'Bearing', 'CentralPressure', 'EnvPressure', 'rMax') -TCRM_FIELD_TYPES = ('N',) * 9 -TCRM_FIELD_WIDTH = (2, 6, 9, 9, 8, 8, 8, 8, 8) -TCRM_FIELD_PREC = (0, 2, 4, 4, 4, 4, 3, 3, 4) + 'Speed', 'Bearing', 'CentralPressure', 'EnvPressure', + 'rMax','Category') +TCRM_FIELD_TYPES = ('N',) * 10 +TCRM_FIELD_WIDTH = (2, 6, 9, 9, 8, 8, 8, 8, 8, 1) +TCRM_FIELD_PREC = (0, 2, 4, 4, 4, 4, 3, 3, 4, 0) TCRM_FIELDS = [[n, t, w, p] for n, t, w, p in zip(TCRM_FIELD_NAMES, TCRM_FIELD_TYPES, @@ -79,6 +80,40 @@ def recdropfields(rec, names): return newrec +def add_field(a, descr): + """ + Add a field to the description of a track. + """ + if a.dtype.fields is None: + raise ValueError("`A' must be a structured numpy array") + b = np.empty(a.shape, dtype=a.dtype.descr + descr) + for name in a.dtype.names: + b[name] = a[name] + return b + + +def add_category(tracks): + """ + Add a category field (for central pressure) to the tracks. + """ + for track in tracks: + track.data = add_field(track.data, [('Category', int)]) + + for rec in track.data: + if rec["CentralPressure"] < 930: + rec["Category"] = 5 + elif rec["CentralPressure"] < 955: + rec["Category"] = 4 + elif rec["CentralPressure"] < 970: + rec["Category"] = 3 + elif rec["CentralPressure"] < 985: + rec["Category"] = 2 + elif rec["CentralPressure"] < 999: + rec["Category"] = 1 + else: + rec["Category"] = 0 + + def tracks2point(tracks, outputFile, netcdf_format=False): """ Writes tracks to a shapefile as a collection of point features. @@ -151,6 +186,7 @@ def tracks2line(tracks, outputFile, dissolve=False, netcdf_format=False): for track in tracks: track.data = recdropfields(track.data, ['Datetime']) + if dissolve: if len(track.data) > 1: dlon = np.diff(track.Longitude) @@ -302,6 +338,7 @@ def tracks2line(tracks, outputFile, dissolve=False, netcdf_format=False): else: raise ValueError("format of {} is not recognizable".format(track_file)) + add_category(tracks) tracks2point(tracks, pt_output_file, netcdf_format=netcdf_format) tracks2line(tracks, line_output_file, netcdf_format=netcdf_format) tracks2line(tracks, dissolve_output_file, dissolve=True, netcdf_format=netcdf_format) diff --git a/setup.py b/setup.py index 88b0ca54..acad2904 100644 --- a/setup.py +++ b/setup.py @@ -7,7 +7,7 @@ setup( name = "TCRM", - version = '3.0rc3', + version = '3.1.0', packages=find_packages(), scripts=['tcrm.py', 'tcevent.py'], include_package_data=True, @@ -20,32 +20,39 @@ install_requires = [ 'numpy', 'scipy', - 'pandas', - 'cartopy', - 'affine', 'matplotlib', 'basemap', + 'shapely', + 'nose', 'netcdf4', 'cftime', - 'configparser', - 'gdal', + 'coverage', + 'coveralls', 'pycurl', 'pyproj', 'seaborn', - 'shapely', 'simplejson', + 'sqlite', 'statsmodels', + 'libgdal' + 'gdal', + 'configparser', + 'cartopy', + 'affine', + 'pandas', 'tqdm', 'xarray', - 'nose', - 'coverage', - 'coveralls'], + 'pthread-stubs', + 'imageio', + 'mpi4py', + 'boto3', + 'botocore'], # metadata: author = "Craig Arthur", author_email = "craig.arthur@ga.gov.au", description = "Tropical Cyclone Risk Model", keywords = "Tropical cyclone risk hazard", - url = "https://geoscienceastralia.github.io/tcrm", + url = "https://geoscienceaustralia.github.io/tcrm", ) diff --git a/tcrmenv.yml b/tcrmenv.yml index e1829615..84d02fa7 100644 --- a/tcrmenv.yml +++ b/tcrmenv.yml @@ -30,3 +30,5 @@ dependencies: - pthread-stubs - imageio - mpi4py + - boto3 + - botocore diff --git a/wind/__init__.py b/wind/__init__.py index a3540943..11b80a6b 100644 --- a/wind/__init__.py +++ b/wind/__init__.py @@ -52,6 +52,7 @@ from Utilities.track import ncReadTrackData, Track from ProcessMultipliers import processMultipliers as pM + class WindfieldAroundTrack(object): """ The windfield around the tropical cyclone track. @@ -86,7 +87,7 @@ class WindfieldAroundTrack(object): :type gustFactor: float :param gustFactor: Conversion from 1-min mean to 0.2-sec gust wind speed, - for off-sea conditions. See WMO TD1555 (2010) for + for off-sea conditions. See WMO TD1555 (2010) for details. :type gridLimit: :class:`dict` @@ -153,7 +154,8 @@ def pressureProfile(self, i, R): from PressureInterface.pressureProfile import PrsProfile as PressureProfile p = PressureProfile(R, convert(self.track.EnvPressure[i], 'hPa', 'Pa'), - convert(self.track.CentralPressure[i], 'hPa', 'Pa'), + convert( + self.track.CentralPressure[i], 'hPa', 'Pa'), self.track.rMax[i], self.track.Latitude[i], self.track.Longitude[i], @@ -183,7 +185,7 @@ def localWindField(self, i): thetaFm = bearing2theta(self.track.Bearing[i] * np.pi/180.) thetaMax = self.thetaMax - #FIXME: temporary way to do this + # FIXME: temporary way to do this cls = windmodels.profile(self.profileType) params = windmodels.profileParams(self.profileType) values = [getattr(self, p) for p in params if hasattr(self, p)] @@ -193,7 +195,7 @@ def localWindField(self, i): P = self.pressureProfile(i, R) - #FIXME: temporary way to do this + # FIXME: temporary way to do this cls = windmodels.field(self.windFieldType) params = windmodels.fieldParams(self.windFieldType) values = [getattr(self, p) for p in params if hasattr(self, p)] @@ -257,28 +259,32 @@ def regionalExtremes(self, gridLimit, timeStepCallback=None): # We only consider the times when the TC track falls in the region timesInRegion = np.where((xMin <= self.track.Longitude) & - (self.track.Longitude <= xMax) & - (yMin <= self.track.Latitude) & - (self.track.Latitude <= yMax))[0] + (self.track.Longitude <= xMax) & + (yMin <= self.track.Latitude) & + (self.track.Latitude <= yMax))[0] nsteps = len(self.track.TimeElapsed) for i in tqdm.tqdm(timesInRegion, disable=None): - log.debug("Calculating wind field at timestep {0} of {1}".format(i, nsteps)) + log.debug(("Calculating wind field at timestep " + "{0} of {1}".format(i, nsteps))) # Map the local grid to the regional grid # Set up max/min over the whole domain - jmin, jmax = 0, int((maxLat - minLat + 2. * gridMargin) \ - / gridStep) + 1 - imin, imax = 0, int((maxLon - minLon + 2. * gridMargin) \ - / gridStep) + 1 + jmin, jmax = 0, int((maxLat - minLat + 2. * gridMargin) / + gridStep) + 1 + imin, imax = 0, int((maxLon - minLon + 2. * gridMargin) / + gridStep) + 1 - # If domain is set in gridLimit in config, use only those limits instead + # If domain is set in gridLimit in config, use only those + # limits instead if self.domain == 'bounded': jmin = int((latCDegree[i] - minLat - gridMargin) / gridStep) - jmax = int((latCDegree[i] - minLat + gridMargin) / gridStep) + 1 + jmax = int((latCDegree[i] - minLat + + gridMargin) / gridStep) + 1 imin = int((lonCDegree[i] - minLon - gridMargin) / gridStep) - imax = int((lonCDegree[i] - minLon + gridMargin) / gridStep) + 1 + imax = int((lonCDegree[i] - minLon + + gridMargin) / gridStep) + 1 # Calculate the local wind speeds and pressure at time i Ux, Vy, P = self.localWindField(i) @@ -384,7 +390,8 @@ def setGridLimit(self, track): """ - track_limits = {'xMin':9999, 'xMax':-9999, 'yMin':9999, 'yMax':-9999} + track_limits = {'xMin': 9999, 'xMax': - + 9999, 'yMin': 9999, 'yMax': -9999} track_limits['xMin'] = min(track_limits['xMin'], track.Longitude.min()) track_limits['xMax'] = max(track_limits['xMax'], track.Longitude.max()) track_limits['yMin'] = min(track_limits['yMin'], track.Latitude.min()) @@ -395,7 +402,6 @@ def setGridLimit(self, track): self.gridLimit['yMin'] = np.floor(track_limits['yMin']) self.gridLimit['yMax'] = np.ceil(track_limits['yMax']) - def calculateExtremesFromTrack(self, track, callback=None): """ Calculate the wind extremes given a single tropical cyclone track. @@ -423,21 +429,20 @@ def calculateExtremesFromTrack(self, track, callback=None): resolution=self.resolution, gridLimit=self.gridLimit, domain=self.domain) - + if self.config.getboolean('Timeseries', 'Windfields', fallback=False): - from . import writer - output = pjoin(self.windfieldPath, + from . import writer + output = pjoin(self.windfieldPath, 'evolution.{0:03d}-{1:05d}.nc'.format( - *track.trackId)) - callback = writer.WriteFoliationCallback(output, - self.gridLimit, - self.resolution, - self.margin, + *track.trackId)) + callback = writer.WriteFoliationCallback(output, + self.gridLimit, + self.resolution, + self.margin, wraps=callback) return track, wt.regionalExtremes(self.gridLimit, callback) - def calculateExtremesFromTrackfile(self, trackfile, callback=None): """ Calculate the wind extremes from a `trackfile` that might contain a @@ -493,23 +498,23 @@ def dumpGustsFromTracks(self, trackiter, windfieldPath, """ if timeStepCallback: results = map(self.calculateExtremesFromTrack, - trackiter, - itertools.repeat(timeStepCallback)) + trackiter, + itertools.repeat(timeStepCallback)) else: results = map(self.calculateExtremesFromTrack, - trackiter) + trackiter) for track, result in results: - log.debug("Saving data for track {0:03d}-{1:05d}"\ + log.debug("Saving data for track {0:03d}-{1:05d}" .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:03d}-{1:05d}.nc'. format(*track.trackId)) plotfile = pjoin(windfieldPath, - 'gust.{0:03d}-{1:05d}.png'.\ + 'gust.{0:03d}-{1:05d}.png'. format(*track.trackId)) self.saveGustToFile(track.trackfile, (np.flipud(lat), lon, @@ -518,10 +523,9 @@ def dumpGustsFromTracks(self, trackiter, windfieldPath, np.flipud(Vy), np.flipud(P)), dumpfile) - + if self.multipliers is not None: self.calcLocalWindfield(results) - #self.plotGustToFile((lat, lon, gust, Vx, Vy, P), plotfile) def plotGustToFile(self, result, filename): """ @@ -703,15 +707,15 @@ def calcLocalWindfield(self, results): """ Calculate a local wind field using the regional windfield data - :param results: collection of :tuple: track and wind field data + :param results: collection of :tuple: track and wind field data """ # Load a multiplier file to determine the projection: - # m4_max_file = pjoin(self.multipliers, 'm4_max.img') + # m4_max_file = pjoin(self.multipliers, 'm4_max.img') log.info("Using M4 data from {0}".format(self.multipliers)) for track, result in results: - log.debug("Doing Multiplier for track {0:03d}-{1:05d}"\ + log.debug("Doing Multiplier for track {0:03d}-{1:05d}" .format(*track.trackId)) # gust, bearing, Vx, Vy, P, lon, lat = result = N # windfield_path = None @@ -721,9 +725,10 @@ def calcLocalWindfield(self, results): lon = lon - delta / 2. lat = lat - delta / 2. - pM.processMult(gust, Vx, Vy, lon, lat,self.windfieldPath, + pM.processMult(gust, Vx, Vy, lon, lat, self.windfieldPath, self.multipliers) + def loadTracksFromFiles(trackfiles): """ Generator that yields :class:`Track` objects from a list of track @@ -828,14 +833,12 @@ def run(configFile, callback=None): trackPath = pjoin(outputPath, 'tracks') gridLimit = None - if config.has_option('Region','gridLimit'): + if config.has_option('Region', 'gridLimit'): gridLimit = config.geteval('Region', 'gridLimit') if config.has_option('WindfieldInterface', 'gridLimit'): gridLimit = config.geteval('WindfieldInterface', 'gridLimit') - - #if callback: - # raise NotImplementedError + if config.getboolean('Timeseries', 'Extract', fallback=False): from Utilities.timeseries import Timeseries ts = Timeseries(configFile) diff --git a/wind/vmax.py b/wind/vmax.py index 3de42529..e09412c9 100644 --- a/wind/vmax.py +++ b/wind/vmax.py @@ -71,7 +71,8 @@ def vmax(pCentre, pEnv, type="holland", beta=1.3, rho=1.15): (Atkinson & Holliday) was determined using surface wind observations so should be used with caution at the gradient level. - :raises ValueError: if environmental pressure is lower than central pressure + :raises ValueError: if environmental pressure is lower than + central pressure Note: The pressure should ideally be passed in units of Pa, but the function will accept hPa and automatically convert to Pa. @@ -115,6 +116,7 @@ def vmax(pCentre, pEnv, type="holland", beta=1.3, rho=1.15): raise NotImplementedError("Vmax type " + type + " not implemented") return vMax + def pDiff(vMax, pEnv, vMaxType="holland", beta=1.3, rho=1.15): """ Inverse functions to calculate central pressure from vMax diff --git a/wind/windmodels.py b/wind/windmodels.py index cbbc46c1..b463fccc 100644 --- a/wind/windmodels.py +++ b/wind/windmodels.py @@ -43,6 +43,7 @@ log = logging.getLogger(__name__) log.addHandler(logging.NullHandler()) + class WindSpeedModel(object): """ @@ -275,8 +276,8 @@ def vorticity(self, R): """ rr = R * 1000. - Z = (np.sign(self.f) * 2 * self.vMax * - self.rMax / (self.rMax ** 2 + rr ** 2) + + Z = (np.sign(self.f) * 2 * self.vMax * + self.rMax / (self.rMax ** 2 + rr ** 2) + np.sign(self.f) * 2 * self.vMax * self.rMax * (self.rMax ** 2 - rr ** 2) / (self.rMax ** 2 + rr ** 2) ** 2) @@ -324,9 +325,9 @@ def secondDerivative(self): d2Vm = ((beta * dP * (-4 * beta ** 3 * dP / rho - (-2 + beta ** 2) * E * (np.abs(f) * rMax) ** 2)) / - (E * rho * sqrt((4 * beta * dP) / (E * rho) - + (f * rMax) ** 2) * (4 * beta * dP * rMax ** 2 / rho - + E * (f * rMax ** 2) ** 2))) + (E * rho * sqrt((4 * beta * dP) / (E * rho) + + (f * rMax) ** 2) * (4 * beta * dP * rMax ** 2 / rho + + E * (f * rMax ** 2) ** 2))) try: assert d2Vm < 0.0 @@ -349,9 +350,10 @@ def firstDerivative(self): E = exp(1) - dVm = (-np.abs(f)/2 + (E*(f**2)*rMax*np.sqrt((4*beta*dP/rho)/E + \ - (f*rMax)**2))/ \ - (2*(4*beta*dP/rho + E*(f*rMax)**2))) + dVm = (-np.abs(f)/2 + (E * (f**2) * rMax * + np.sqrt((4 * beta * dP / rho) / E + + (f * rMax) ** 2)) / + (2 * (4 * beta * dP / rho + E * (f * rMax)**2))) return dVm def velocity(self, R): @@ -369,8 +371,8 @@ def velocity(self, R): d2Vm = self.secondDerivative() dVm = self.firstDerivative() - aa = ((d2Vm / 2. - (dVm - self.vMax /self.rMax) - / self.rMax) / self.rMax) + 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 @@ -388,7 +390,7 @@ def velocity(self, R): def vorticity(self, R): """ Calculate the vorticity associated with the (gradient level) - vortex. + vortex. :param R: :class:`numpy.ndarray` of distance of grid from the TC centre (metres). @@ -402,11 +404,14 @@ def vorticity(self, R): 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 + 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() @@ -448,8 +453,8 @@ def __init__(self, lat, lon, eP, cP, rMax, windSpeedModel=WilloughbyWindSpeed): HollandWindProfile.__init__(self, lat, lon, eP, cP, rMax, 1.0, windSpeedModel) - self.beta = (1.0036 + 0.0173 * self.vMax - 0.0313 * np.log(rMax/1000.) - + 0.0087 * np.abs(lat)) + self.beta = (1.0036 + 0.0173 * self.vMax - 0.0313 * + np.log(rMax/1000.) + 0.0087 * np.abs(lat)) self.speed = HollandWindSpeed(self) @@ -511,10 +516,10 @@ def vorticity(self, R): """ V = self.velocity(R) - Z = V/R - self.alpha*self.vMax*self.rMax/(R**(self.alpha+1)) + Z = V/R - self.alpha * self.vMax * self.rMax / (R**(self.alpha + 1)) icore = np.where(R <= self.rMax) Z[icore] = np.sign(self.f) * ((self.vMax * (R[icore] / - self.rMax) + self.vMax / self.rMax)) + self.rMax) + self.vMax / self.rMax)) return Z @@ -611,33 +616,33 @@ def secondDerivative(self): (8 * (4 * beta1 * dp1 / (rho * E) + (4 * beta2 * dp2 / rho) * nu * exp(-nu) + - (rMax1 * f) ** 2) ** 1.5) - * (-(4 * (beta1 ** 2) * dp1 / (rho * rMax1 * E)) + + (rMax1 * f) ** 2) ** 1.5) * + (-(4 * (beta1 ** 2) * dp1 / (rho * rMax1 * E)) + (4 * (beta1 ** 2) * dp1 / (rho * rMax1 * E)) - (4 * (beta2 ** 2) * dp2 / rho) * - (nu / rMax1) * exp(-nu) - + (4 * (beta2 ** 2) * dp2 / rho) * - ((nu ** 2) / rMax1) * exp(-nu) - + 2 * rMax1 * f ** 2) ** 2 - + 1 / (4 * sqrt((4 * beta1 * dp1 / (rho * E)) + - (4 * beta2 * dp2 / rho) * nu * 2 + - exp(-nu) + (rMax1 * f) ** 2)) - * ((4 * (beta1 ** 3) * dp1 / (rho * (rMax1 ** 2) * E)) - + (4 * (beta1 ** 2) * dp1 / (rho * (rMax1 ** 2) * E)) - - (12 * (beta1 ** 3) * dp1 / (rho * (rMax1 ** 2) * E)) - - (4 * (beta1 ** 2) * dp1 / (rho * (rMax1 ** 2) * E)) - + (4 * (beta1 ** 3) * dp1 / (rho * (rMax1 ** 2) * E)) - + (4 * (beta2 ** 3) * dp2 / rho) * - (nu / (rMax1 ** 2)) * exp(-nu) - + (4 * (beta2 ** 2) * dp2 / rho) * - (nu / (rMax1 ** 2)) * exp(-nu) - - (12 * (beta2 ** 3) * dp2 / rho) * - (nu ** 2) / (rMax1 ** 2) * exp(-nu) - - (4 * (beta2 ** 2) * dp2 / rho) * - (nu ** 2) / (rMax1 ** 2) * exp(-nu) - + (4 * (beta2 ** 3) * dp2 / rho) * - (nu ** 3) / (rMax1 ** 2) * exp(-nu) - + 2 * f ** 2)) + (nu / rMax1) * exp(-nu) + + (4 * (beta2 ** 2) * dp2 / rho) * + ((nu ** 2) / rMax1) * exp(-nu) + + 2 * rMax1 * f ** 2) ** 2 + + 1 / (4 * sqrt((4 * beta1 * dp1 / (rho * E)) + + (4 * beta2 * dp2 / rho) * nu * 2 + + exp(-nu) + (rMax1 * f) ** 2)) * + ((4 * (beta1 ** 3) * dp1 / (rho * (rMax1 ** 2) * E)) + + (4 * (beta1 ** 2) * dp1 / (rho * (rMax1 ** 2) * E)) - + (12 * (beta1 ** 3) * dp1 / (rho * (rMax1 ** 2) * E)) - + (4 * (beta1 ** 2) * dp1 / (rho * (rMax1 ** 2) * E)) + + (4 * (beta1 ** 3) * dp1 / (rho * (rMax1 ** 2) * E)) + + (4 * (beta2 ** 3) * dp2 / rho) * + (nu / (rMax1 ** 2)) * exp(-nu) + + (4 * (beta2 ** 2) * dp2 / rho) * + (nu / (rMax1 ** 2)) * exp(-nu) - + (12 * (beta2 ** 3) * dp2 / rho) * + (nu ** 2) / (rMax1 ** 2) * exp(-nu) - + (4 * (beta2 ** 2) * dp2 / rho) * + (nu ** 2) / (rMax1 ** 2) * exp(-nu) + + (4 * (beta2 ** 3) * dp2 / rho) * + (nu ** 3) / (rMax1 ** 2) * exp(-nu) + + 2 * f ** 2)) assert d2Vm < 0.0 @@ -999,7 +1004,7 @@ def field(self, R, lam, vFm, thetaFm, thetaMax=0.): thetaMaxAbsolute = np.array(thetaFm) + thetaMax * np.pi/180. phi = inflow - lam - asym = (0.5 * (1. + np.cos(thetaMaxAbsolute - lam)) * + asym = (0.5 * (1. + np.cos(thetaMaxAbsolute - lam)) * vFm * (V / np.abs(V).max())) Vsf = V + asym @@ -1054,52 +1059,54 @@ def field(self, R, lam, vFm, thetaFm, thetaMax=0.): else: Umod = vFm Vt = Umod * np.ones(V.shape) - + core = np.where(R > 2. * self.rMax) - Vt[core] = Umod * np.exp(-((R[core] / (2.*self.rMax)) - 1.) ** 2. ) + Vt[core] = Umod * np.exp(-((R[core] / (2.*self.rMax)) - 1.) ** 2.) - al = ((2. * V / R ) + self.f) / (2. * K) + al = ((2. * V / R) + self.f) / (2. * K) be = (self.f + Z) / (2. * K) gam = (V / (2. * K * R)) - + albe = np.sqrt(al / be) ind = np.where(np.abs(gam) > np.sqrt(al * be)) chi = np.abs((Cd / K) * V / np.sqrt(np.sqrt(al * be))) eta = np.abs((Cd / K) * V / np.sqrt(np.sqrt(al * be) + np.abs(gam))) - psi = np.abs((Cd / K) * V / np.sqrt(np.abs(np.sqrt(al * be) - - np.abs(gam)))) + psi = np.abs((Cd / K) * V / np.sqrt(np.abs(np.sqrt(al * be) - + np.abs(gam)))) i = complex(0., 1.) - A0 = -(chi * (1 + i * (1 + chi))*V) / (2*chi**2 + 3*chi +2) - + A0 = -(chi * (1 + i * (1 + chi)) * V) / (2 * chi**2 + 3 * chi + 2) + # Symmetric surface wind component - u0s = np.sign(self.f) * albe * A0.real - v0s = A0.imag + u0s = A0.real * albe * np.sign(self.f) + v0s = A0.imag Am = -(psi * (1 + 2 * albe + (1 + i) * (1 + albe) * eta) * Vt) / \ - (albe * ((2 + 2 * i) * (1 + eta * psi) + 3 * psi + 3* i * eta)) + (albe * ((2 + 2 * i) * (1 + eta * psi) + 3 * psi + 3 * i * eta)) AmIII = -(psi * (1 + 2 * albe + (1 + i) * (1 + albe) * eta) * Vt) / \ - (albe * ((2 - 2 * i + 3 * (eta + psi) + (2 + 2 * i)*eta * psi))) + (albe * ((2 - 2 * i + 3 * (eta + psi) + (2 + 2 * i) * + eta * psi))) Am[ind] = AmIII[ind] # First asymmetric surface component - ums = albe * (Am * np.exp(-i * lam * np.sign(self.f))).real - vms = np.sign(self.f) * (Am * np.exp(-i * lam * np.sign(self.f))).imag - - Ap = -(eta * (1 - 2 * albe + (1 + i)*(1 - albe) * psi) * Vt) / \ - (albe * ((2 + 2*i)*(1 + eta * psi) + 3*eta + 3*i*psi)) - ApIII = -(eta * (1 - 2 * albe + (1 - i)*(1 - albe)*psi)*Vt) / \ - (albe * (2 + 2 * i + 3 * (eta + psi) + (2 -2 * i)*eta*psi)) + 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) + + Ap = -(eta * (1 - 2 * albe + (1 + i) * (1 - albe) * psi) * Vt) / \ + (albe * ((2 + 2 * i) * (1 + eta * psi) + 3 * eta + 3 * i * psi)) + ApIII = -(eta * (1 - 2 * albe + (1 - i) * (1 - albe)*psi) * Vt) / \ + (albe * (2 + 2 * i + 3 * (eta + psi) + (2 - 2 * i) * + eta * psi)) Ap[ind] = ApIII[ind] # Second asymmetric surface component - ups = albe * (Ap * np.exp(i * lam * np.sign(self.f))).real - vps = np.sign(self.f) * (Ap * np.exp(i * lam * np.sign(self.f))).imag + 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) # Total surface wind in (moving coordinate system) - us = u0s + ups + ums - vs = V + v0s + vps + vms + us = u0s + ups + ums + vs = v0s + vps + vms + V usf = us + Vt * np.cos(lam - thetaFm) vsf = vs - Vt * np.sin(lam - thetaFm) diff --git a/wind/writer.py b/wind/writer.py index c999f77a..6204c8e6 100644 --- a/wind/writer.py +++ b/wind/writer.py @@ -7,6 +7,7 @@ import affine import numpy as np + class WriteFoliationCallback(object): """ Incremental (by time step) recording of wind fields to NetCDF @@ -16,59 +17,77 @@ class WriteFoliationCallback(object): >>> with tempfile.NamedTemporaryFile(delete=False) as f: ... res = 0.5 ... lat, lon = np.r_[-10:-5+.1:res], np.r_[110:115+.1:res] - ... gridLimit = dict(xMin=min(lon), yMin=min(lat), yMax=max(lat), xMax=max(lon)) + ... gridLimit = dict(xMin=min(lon), yMin=min(lat), + ... yMax=max(lat), xMax=max(lon)) ... callback = WriteFoliationCallback(f.name, gridLimit, res) - ... callback(datetime.datetime(2018,3,17), *[np.ones((len(lat),len(lon)))]*4, lon=lon, lat=lat) - ... callback(datetime.datetime(2018,3,18), *[np.full((3,3), fill_value=100)]*4, lon=lon[2:5], lat=lat[2:5]) - ... callback(datetime.datetime(2018,3,20,12,0), *[np.ones((3,6))]*4, lon=lon[0:6], lat=lat[2:5]) + ... callback(datetime.datetime(2018,3,17), + ... *[np.ones((len(lat),len(lon)))]*4, lon=lon, lat=lat) + ... callback(datetime.datetime(2018,3,18), *[np.full((3,3), + ... fill_value=100)]*4, lon=lon[2:5], lat=lat[2:5]) + ... callback(datetime.datetime(2018,3,20,12,0), *[np.ones((3,6))]*4, + ... lon=lon[0:6], lat=lat[2:5]) ... callback.ds.close() ... x = xarray.open_dataset(f.name) >>> len(x.time) 3 >>> x.pressure.sum(dim=['lat','lon']).values.tolist() [121.0, 900.0, 18.0] - >>> (x.gust_speed.isel(time=2, lon=slice(0,6), lat=slice(2,5)).values == 1).all() + >>> (x.gust_speed.isel(time=2, lon=slice(0,6), + ... lat=slice(2,5)).values == 1).all() True + + Uses low-level NetCDF bindings to support windowed progressive writes. """ - # Uses low-level NetCDF bindings to support windowed progressive writes. - def __init__(self, filename, gridLimit, resolution, margin=0, maxchunk=256, wraps=None): + def __init__(self, filename, gridLimit, resolution, margin=0, + maxchunk=256, wraps=None): """ :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) + :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)) + + logging.debug( + "Preparing to record windfield evolution to {}".format(filename)) self.callback = wraps def series(start, stop, inc=resolution): return np.linspace(start, stop, int(round((stop-start)/inc)) + 1) + lat = series(gridLimit['yMin'] - margin, gridLimit['yMax'] + margin) lon = series(gridLimit['xMin'] - margin, gridLimit['xMax'] + margin) - self.affine = ~(affine.Affine.translation(xoff=gridLimit['xMin']-margin, - yoff=gridLimit['yMin']-margin) - * affine.Affine.scale(resolution)) + self.affine = ~(affine.Affine.translation( + xoff=gridLimit['xMin']-margin, + yoff=gridLimit['yMin']-margin) * + affine.Affine.scale(resolution)) self.ds = root = netCDF4.Dataset(filename, mode='w') root.description = "Simulated Windfield Timeseries" # declare shapes - root.createDimension('time', size=None) # growable + root.createDimension('time', size=None) root.createDimension('lat', size=len(lat)) root.createDimension('lon', size=len(lon)) # coordinate variables: - self.time = root.createVariable('time', datatype='f8', dimensions=('time',)) - self.lat = root.createVariable('lat', datatype='f4', dimensions=('lat',)) - self.lon = root.createVariable('lon', datatype='f4', dimensions=('lon',)) + self.time = root.createVariable('time', + datatype='f8', + dimensions=('time',)) + self.lat = root.createVariable('lat', + datatype='f4', + dimensions=('lat',)) + self.lon = root.createVariable('lon', + datatype='f4', + dimensions=('lon',)) self.lat[:] = lat self.lon[:] = lon @@ -76,8 +95,8 @@ def series(start, stop, inc=resolution): # data variables: etc = dict(datatype='f4', dimensions=('time', 'lat', 'lon'), - #fill_value=0, # np.NaN ? - chunksizes=(1, min(len(lat), maxchunk), min(len(lon), maxchunk)), + chunksizes=(1, min(len(lat), maxchunk), + min(len(lon), maxchunk)), zlib=True) self.speed = root.createVariable('gust_speed', fill_value=0, **etc) self.Ux = root.createVariable('velocity_east', fill_value=0, **etc) @@ -89,9 +108,9 @@ def __call__(self, time, gust, Ux, Uy, P, lon, lat): if self.callback: self.callback(time, gust, Ux, Uy, P, lon, lat) - t = len(self.time) # current time index + t = len(self.time) - if not t: # then initialise + if not t: self.time.units = "days since " + time.date().isoformat() # convert window extent to slice indices @@ -106,4 +125,5 @@ def __call__(self, time, gust, Ux, Uy, P, lon, lat): self.Uy[t, i, j] = Uy self.P[t, i, j] = P - self.ds.sync() # save (flush) layer to file + # save (flush) layer to file + self.ds.sync()