diff --git a/Evaluate/plotStatisics.py b/Evaluate/plotStatisics.py new file mode 100644 index 00000000..1a91e578 --- /dev/null +++ b/Evaluate/plotStatisics.py @@ -0,0 +1,373 @@ +#!/usr/bin/env python +# coding: utf-8 + +# # Plotting up the statistical output +# +# To gain an understanding of the meaning of the statistics, and the spatial distribution of these statistics, we plot up several of the metrics on maps. This allows you to explore how these variables vary over the simulation domain, and interpret the meaning of these variables. +# +# This only gives one perspective on the choice of variables used in TCRM to simulate TC track behaviour. There are a number of other graphical products that further help with understanding the way TCRM works. + +import logging +import io + +from Utilities.config import ConfigParser +from Utilities.nctools import ncLoadFile, ncGetDims, ncGetVar +from Utilities.grid import SampleGrid +from PlotInterface.maps import ArrayMapFigure, saveFigure +import numpy as np + +from os.path import join as pjoin +import seaborn as sns +sns.set_context('paper') + + +# Use the normal approach of setting up a configuration string, +# and set it up to read data from the output folders containing +# the data we've generated in previous notebooks. + +configstr = """ +[DataProcess] +InputFile=C:/WorkSpace/data/Allstorms.ibtracs_wmo.v03r10.csv +Source=IBTRACS +StartSeason=1981 +FilterSeasons=False + +[Region] +; Domain for windfield and hazard calculation +gridLimit={'xMin':90.,'xMax':180.,'yMin':-35.0,'yMax':-5.0} +gridSpace={'x':1.0,'y':1.0} +gridInc={'x':1.0,'y':0.5} + +[TrackGenerator] +NumSimulations=5000 +YearsPerSimulation=10 +SeasonSeed=68876543 +TrackSeed=334825 +TimeStep=1.0 + +[Input] +landmask = C:/WorkSpace/tcrm/input/landmask.nc +mslpfile = C:/WorkSpace/tcrm/MSLP/slp.day.ltm.nc +datasets = IBTRACS,LTMSLP + +[Output] +Path=C:/WorkSpace/data/tcha/ + +[Hazard] +Years=2,5,10,20,25,50,100,200,250,500,1000,2000,2500,5000 +MinimumRecords=10 +CalculateCI=False + +[Logging] +LogFile=C:/WorkSpace/data/tcha/log/tcha.log +LogLevel=INFO +Verbose=False + +[IBTRACS] +; Input data file settings +url = ftp://eclipse.ncdc.noaa.gov/pub/ibtracs/v03r10/wmo/csv/Allstorms.ibtracs_wmo.v03r10.csv.gz +path = C:/WorkSpace/data/ +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 = C:/WorkSpace/data/MSLP +filename = slp.day.ltm.nc +""" + +config = ConfigParser() +config.readfp(io.StringIO(configstr)) + +# We load a landmask dataset to allow us to determine +# when TCs are over water or over land + +landmask_file = config.get('Input', 'Landmask') +landmask = SampleGrid(landmask_file) + + +# One thing that we realise in the statistics is that TC behaviour +# changes when a TC makes landfall. And because the number of +# observations in any given cell may not be sufficient to calculate +# reliable statistics, TCRM automatically increases the region that +# is sampled. + +# This causes problems in regions close to land. If the model is +# determining statistics for one of these cells close to the coast +# (but offshore), and the expanded region starts capturing observations +# from over land, then the statistics are not truly representative of +# behaviour of offshore TCs. This is important for parameters related +# to intensity - TCs can often continue intensifying right up to landfall. +# But if we are sampling statistics of TCs overland, then we dilute the +# intensity statistics with observations from overland which (nearly without +# exception) act to weaken the storm. + +# To avoid any potential issues like this, TCRM flags over-land +# observations, and the statistics are calculated separately for +# over-water and over-land observations. When it comes to viewing +# the results, we need to read in the over-water and over-land +# observations separately, then combine them into a single grid +# for plotting. This next function does exactly that. + +def getData(ncobj, varname, ij): + var = ncGetVar(ncobj, varname)[:] + lvar = ncGetVar(ncobj, 'l'+varname)[:] + data = var + data[ij] = lvar[ij] + return data + + +# Let's breifly look at the statistical data generated. Here, we +# load the 'pressure_rate_stats.nc' file, which contains the statistics +# for the pressure rate of change (hPa/hr) of obsersved TCs. + +outputPath = config.get("Output", "Path") +processPath = pjoin(outputPath, "process") +plotsPath = pjoin(outputPath, "plots", "stats") +fname = pjoin(processPath, "pressure_rate_stats.nc") +ncobj = ncLoadFile(fname) + +lon = ncGetDims(ncobj, 'lon') +lat = ncGetDims(ncobj, 'lat') +ncobj.close() +xgrid, ygrid = np.meshgrid(lon, lat) + +ls = np.zeros(np.shape(xgrid)) + +for i in range(len(lon)): + for j in range(len(lat)): + if landmask.sampleGrid(lon[i], lat[j]) > 0.0: + ls[j, i] = 1 + +ij = np.where(ls==1) + +# Set the map keyword arguments that will help draw the basemap. + +map_kwargs = dict(llcrnrlon=90., llcrnrlat=-30, + urcrnrlon=180., urcrnrlat=-5., + resolution='h', projection='merc') + + +ncobj = ncLoadFile(pjoin(processPath, "pressure_stats.nc")) +lon = ncGetDims(ncobj, 'lon') +lat = ncGetDims(ncobj, 'lat') +ardata = getData(ncobj, 'alpha', ij) +mudata = getData(ncobj, 'mu', ij) +mindata = getData(ncobj, 'min', ij) +sigdata = getData(ncobj, 'sig', ij) +ncobj.close() +fig = ArrayMapFigure() +fig.add(mudata, xgrid, ygrid, 'Mean pressure ', [950, 1000], + 'Pressure (hPa)', map_kwargs) +fig.add(mindata, xgrid, ygrid, 'Minimum pressure', [900, 1000], + 'Pressure (hPa)', map_kwargs) +fig.add(sigdata, xgrid, ygrid, 'Pressure standard deviation', [0, 50], + 'Std dev.', map_kwargs) +fig.add(ardata, xgrid, ygrid, 'Pressure AR(1)', [-1, 1], + 'AR(1)', map_kwargs) +fig.plot() +saveFigure(fig, pjoin(plotsPath, "pressure_stats.png")) + + +ncobj = ncLoadFile(pjoin(processPath, "pressure_stats.nc")) +lon = ncGetDims(ncobj, 'lon') +lat = ncGetDims(ncobj, 'lat') +ardata = getData(ncobj, 'alpha', ij) +mudata = getData(ncobj, 'mu', ij) +mindata = getData(ncobj, 'min', ij) +sigdata = getData(ncobj, 'sig', ij) +ncobj.close() +fig = ArrayMapFigure() +fig.add(mudata, xgrid, ygrid, 'Mean pressure ', [950, 1000], + 'Pressure (hPa)', map_kwargs) +fig.plot() +fig.axes[0].set_ylim((-35, -5)) +saveFigure(fig, pjoin(plotsPath, "pressure_mean.png")) + +fig = ArrayMapFigure() +fig.add(mindata, xgrid, ygrid, 'Minimum pressure', [900, 1000], + 'Pressure (hPa)', map_kwargs) +fig.plot() +fig.axes[0].set_ylim((-35, -5)) +saveFigure(fig, pjoin(plotsPath, "pressure_min.png")) + +fig = ArrayMapFigure() +fig.add(sigdata, xgrid, ygrid, 'Pressure standard deviation', [0, 50], + r'$\sigma$ (hPa)', map_kwargs) +fig.plot() +fig.axes[0].set_ylim((-35, -5)) +saveFigure(fig, pjoin(plotsPath, "pressure_std.png")) + +fig = ArrayMapFigure() +fig.add(ardata, xgrid, ygrid, 'Pressure AR(1)', [-1, 1], + 'AR(1)', map_kwargs) +fig.plot() +fig.axes[0].set_ylim((-35, -5)) +saveFigure(fig, pjoin(plotsPath, "pressure_ar1.png")) + +# The first figure shows the mean pressure (in hPa) in each cell across the simulation domain. It's important to recognise that in areas with few observations, the mean may not be well-defined. TCRM makes an effort to address this, by expanding the sampling region if there are insufficient observations in an individual cell. +# +# Next is the minimum pressure. This is where we can see some significant variability, because the minimum pressure in any given cell depends on the vagaries of historical TC behaviour, and we are looking at an extreme value in the distribution of pressure values in each cell. +# +# The third figure is the standard deviation of the pressure observations. Areas where there's a lot of variation are likely near coastlines and at the outer limits of the simulation region, where TC intensity is more variable. +# +# The final figure shows the lag-1 autocorrelation of pressure observations. For this variable, the values are all very close to 1, because the pressure values do not change dramatically from one observation to the next. For example, the largest change in pressure over a 6-hour period is only about *XX* hPa - compare that to the pressure observation and the change is (in a numerical sense) quite small. But this is not the case for all the model variables, as we will see below. + +# Now lets examine the pressure rate statistics. This is measuring the intensification rate of TCs in each cell, rather than the intensity (above). +# +# The first figure is the mean pressure rate of change. Positive values indicate an increase in pressure, negative values a decrease in pressure. This shows us the (historically) favoured regions for intensification and weakening of TCs. +# +# In figure 4 is the lag-1 autocorrelation. Compare the values plotted here with the corresponding figure for the pressure variable (above). We saw previously that the pressure rate of change is largely uncorrelated when looking at the autocorrelation. But what spatial patterns emerge in the pressure rate statistics? + + + +ncobj = ncLoadFile(pjoin(processPath, "pressure_rate_stats.nc")) +lon = ncGetDims(ncobj, 'lon') +lat = ncGetDims(ncobj, 'lat') +ardata = getData(ncobj, 'alpha', ij) +mudata = getData(ncobj, 'mu', ij) +mindata = getData(ncobj, 'min', ij) +sigdata = getData(ncobj, 'sig', ij) +ncobj.close() + +fig = ArrayMapFigure() +fig.add(mudata, xgrid, ygrid, 'Mean pressure rate', [-1, 1], + 'Pressure rate (hPa/hr)', map_kwargs) +fig.add(mindata, xgrid, ygrid, 'Minimum pressure rate', [-10, 10], + 'Pressure rate (hPa/hr)', map_kwargs) +fig.add(sigdata, xgrid, ygrid, 'Pressure rate standard deviation', [0, 2.5], + r'$\sigma$ (hPa/hr)', map_kwargs) +fig.add(ardata, xgrid, ygrid, 'Pressure rate AR(1)', [-1, 1], + 'AR(1)', map_kwargs) +fig.plot() +saveFigure(fig, pjoin(plotsPath, "pressure_rate_stats.png")) + +ncobj = ncLoadFile(pjoin(processPath, "pressure_rate_stats.nc")) +lon = ncGetDims(ncobj, 'lon') +lat = ncGetDims(ncobj, 'lat') +ardata = getData(ncobj, 'alpha', ij) +mudata = getData(ncobj, 'mu', ij) +mindata = getData(ncobj, 'min', ij) +sigdata = getData(ncobj, 'sig', ij) +ncobj.close() + +fig = ArrayMapFigure() +fig.add(mudata, xgrid, ygrid, 'Mean pressure rate', [-1, 1], + 'Pressure rate (hPa/hr)', map_kwargs) +fig.plot() +fig.axes[0].set_ylim((-35, -5)) +saveFigure(fig, pjoin(plotsPath, "pressure_rate_mean.png")) + +fig = ArrayMapFigure() +fig.add(mindata, xgrid, ygrid, 'Minimum pressure rate', [-10, 10], + 'Pressure rate (hPa/hr)', map_kwargs) +fig.plot() +fig.axes[0].set_ylim((-35, -5)) +saveFigure(fig, pjoin(plotsPath, "pressure_rate_min.png")) +fig = ArrayMapFigure() +fig.add(sigdata, xgrid, ygrid, 'Pressure rate standard deviation', [0, 2.5], + r'$\sigma$ (hPa/hr)', map_kwargs) +fig.plot() +fig.axes[0].set_ylim((-35, -5)) +saveFigure(fig, pjoin(plotsPath, "pressure_rate_std.png")) + +fig = ArrayMapFigure() +fig.add(ardata, xgrid, ygrid, 'Pressure rate AR(1)', [-1, 1], + 'AR(1)', map_kwargs) +fig.plot() +fig.axes[0].set_ylim((-35, -5)) +saveFigure(fig, pjoin(plotsPath, "pressure_rate_ar1.png")) + + + + +ncobj = ncLoadFile(pjoin(processPath, "speed_rate_stats.nc")) +lon = ncGetDims(ncobj, 'lon') +lat = ncGetDims(ncobj, 'lat') +ardata = getData(ncobj, 'alpha', ij) +mudata = getData(ncobj, 'mu', ij) +mindata = getData(ncobj, 'min', ij) +sigdata = getData(ncobj, 'sig', ij) +ncobj.close() + + +fig = ArrayMapFigure() +fig.add(mudata, xgrid, ygrid, 'Mean speed rate', [-1, 1], 'Speed rate (m/s/hr)', map_kwargs) +fig.add(mindata, xgrid, ygrid, 'Minimum speed rate', [-10, 10], 'Speed rate (m/s/hr)', map_kwargs) +fig.add(sigdata, xgrid, ygrid, 'Speed standard deviation', [0, 5], r'$\sigma$ (m/s/hr)', map_kwargs) +fig.add(ardata, xgrid, ygrid, 'Speed rate AR(1)', [-1, 1], 'AR(1)', map_kwargs) +fig.plot() +saveFigure(fig, pjoin(plotsPath, "speed_rate_stats.png")) + + +fig = ArrayMapFigure() +fig.add(mudata, xgrid, ygrid, 'Mean speed', [0, 25], 'Speed (m/s)', map_kwargs) +fig.plot() +fig.axes[0].set_ylim((-35, -5)) +saveFigure(fig, pjoin(plotsPath, "speed_rate_mean.png")) + +fig = ArrayMapFigure() +fig.add(mindata, xgrid, ygrid, 'Minimum speed', [0, 25], 'Speed (m/s)', map_kwargs) +fig.plot() +fig.axes[0].set_ylim((-35, -5)) +saveFigure(fig, pjoin(plotsPath, "speed_rate_min.png")) + +fig = ArrayMapFigure() +fig.add(sigdata, xgrid, ygrid, 'Speed standard deviation', [0, 20], r'$\sigma$ (m/s/hr)', map_kwargs) +fig.plot() +fig.axes[0].set_ylim((-35, -5)) +saveFigure(fig, pjoin(plotsPath, "speed_rate_std.png")) + +fig = ArrayMapFigure() +fig.add(ardata, xgrid, ygrid, 'Speed AR(1)', [-1, 1], 'AR(1)', map_kwargs) +fig.plot() +fig.axes[0].set_ylim((-35, -5)) +saveFigure(fig, pjoin(plotsPath, "speed_rate_ar1.png")) + + +ncobj = ncLoadFile(pjoin(processPath, "bearing_rate_stats.nc")) +lon = ncGetDims(ncobj, 'lon') +lat = ncGetDims(ncobj, 'lat') +ardata = getData(ncobj, 'alpha', ij) +mudata = getData(ncobj, 'mu', ij) +mindata = getData(ncobj, 'min', ij) +sigdata = getData(ncobj, 'sig', ij) +ncobj.close() + +fig = ArrayMapFigure() +fig.add(mudata, xgrid, ygrid, 'Mean bearing rate', [-10, 20], + 'Bearing rate (degrees/hr)', map_kwargs) +fig.add(mindata, xgrid, ygrid, 'Minimum bearing rate', [-100., 10], + 'Bearing rate (degrees/hr)', map_kwargs) +fig.add(sigdata, xgrid, ygrid, 'Bearing rate standard deviation', [0, 1], + r'$\sigma$ (degrees/hr)', map_kwargs) +fig.add(ardata, xgrid, ygrid, 'Bearing rate AR(1)', [-1, 1], + 'AR(1)', map_kwargs) +fig.plot() +saveFigure(fig, pjoin(plotsPath, "bearing_rate_stats.png")) + +ncobj = ncLoadFile(pjoin(processPath, "bearing_stats.nc")) +lon = ncGetDims(ncobj, 'lon') +lat = ncGetDims(ncobj, 'lat') +ardata = getData(ncobj, 'alpha', ij) +mudata = getData(ncobj, 'mu', ij) +mindata = getData(ncobj, 'min', ij) +sigdata = getData(ncobj, 'sig', ij) + +ncobj.close() +fig = ArrayMapFigure() +fig.add(mudata*180./np.pi, xgrid, ygrid, 'Mean bearing', [0, 360.], 'Bearing (degrees)', map_kwargs) +fig.add(mindata*180, xgrid, ygrid, 'Minimum bearing', [0, 180], 'Bearing (degrees)', map_kwargs) +fig.add(sigdata, xgrid, ygrid, 'Bearing standard deviation', [0, 1], 'Std dev.', map_kwargs) +fig.add(ardata, xgrid, ygrid, 'Bearing AR(1)', [-1, 1], 'AR(1)', map_kwargs) +fig.plot() + +bearing = mudata +saveFigure(fig, pjoin(plotsPath, "bearing_stats.png")) diff --git a/PlotInterface/maps.py b/PlotInterface/maps.py index 97ef598b..4d3e59ff 100644 --- a/PlotInterface/maps.py +++ b/PlotInterface/maps.py @@ -193,8 +193,9 @@ def addGraticule(self, axes, mapobj): meridians = np.arange(xmin // dl * dl, xmax + dl, dl) parallels = np.arange(ymin // dl * dl, ymax + dl, dl) - mapobj.gridlines(xlocs=meridians, ylocs=parallels, - draw_labels=True) + gl = mapobj.gridlines(xlocs=meridians, ylocs=parallels, + draw_labels=True) + gl.xlabels_top = False def addCoastline(self, mapobj): """ @@ -438,8 +439,8 @@ def subplot(self, axes, subfigure): vmax = datarange[1] CS = mapobj.pcolormesh(mx, my, data, vmin=vmin, vmax=vmax, cmap=cmap) - #CB = self.colorbar(CS, location='bottom', pad=0.1, - # fig=self, ax=axes) + CB = self.colorbar(CS, orientation='horizontal', pad=0.1, + ax=axes) CB.set_label(cbarlab) axes.set_title(title) self.addGraticule(axes, mapobj) diff --git a/Utilities/loadData.py b/Utilities/loadData.py index bc1981ae..b5b3a7cc 100644 --- a/Utilities/loadData.py +++ b/Utilities/loadData.py @@ -171,7 +171,7 @@ def getSpeedBearing(index, lon, lat, deltatime, ieast=1, -1 = positive longiture westwards. :param missingValue: Replace questionable values with `missingValue`. - :type missingValue: int or float, default = `sys.maxint` + :type missingValue: int or float, default = `sys.maxsize` :returns: speed and bearing : :class:`numpy.ndarray` @@ -729,7 +729,7 @@ def filterPressure(pressure, inputPressureUnits='hPa', :param missingValue: replace all null values in the input data with this value. :type pressure: :class:`numpy.ndarray` - :type missingValue: int or float (default ``sys.maxint``) + :type missingValue: int or float (default ``sys.maxsize``) :returns: :class:`numpy.ndarray` with only valid pressure values. @@ -752,7 +752,7 @@ def getMinPressure(track, missingValue=sys.maxsize): :param track: A :class:`Track` instance :param missingValue: Replace missing values with this value - (default ``sys.maxint``). + (default ``sys.maxsize``). :returns: :class:`Track.trackMinPressure` attribute updated @@ -771,7 +771,7 @@ def getMaxWind(track, missingValue=sys.maxsize): :param track: A :class:`Track` instance :param missingValue: replace all null values in the input data with this value. - :type missingValue: int or float (default ``sys.maxint``) + :type missingValue: int or float (default ``sys.maxsize``) :returns: :class:`Track.trackMaxWind` attribute updated with calculated wind speed updated. diff --git a/Utilities/track.py b/Utilities/track.py index d1b741f9..1e6d94e1 100644 --- a/Utilities/track.py +++ b/Utilities/track.py @@ -101,6 +101,9 @@ def __getattr__(self, key): return self.data[key] + def __repr__(self): + return "".format(", ".join(self.data.dtype.names)) + def inRegion(self, gridLimit): """ Check if the tropical cyclone track starts within a region. diff --git a/Utilities/tracks2shp.py b/Utilities/tracks2shp.py index a0f0adc8..16a3ecc9 100644 --- a/Utilities/tracks2shp.py +++ b/Utilities/tracks2shp.py @@ -31,6 +31,17 @@ OBSFIELD_WIDTH, 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) + +TCRM_FIELDS = [[n, t, w, p] for n, t, w, p in zip(TCRM_FIELD_NAMES, + TCRM_FIELD_TYPES, + TCRM_FIELD_WIDTH, + TCRM_FIELD_PREC)] + # For storing events as a single polyline: EVENTFIELD_NAMES = ('TCID', 'Year', 'Month', 'Day', 'Hour', 'Minute', 'Age', 'MinCP', 'MaxWind' ) @@ -68,7 +79,7 @@ def recdropfields(rec, names): return newrec -def tracks2point(tracks, outputFile): +def tracks2point(tracks, outputFile, netcdf_format=False): """ Writes tracks to a shapefile as a collection of point features. @@ -76,6 +87,7 @@ def tracks2point(tracks, outputFile): :param tracks: :class:`Track` features to store in a shape file :param str outputFile: Path to output file destination + :param bool netcdf_format: Whether tracks are in TCRM format :raises: :mod:`shapefile.ShapefileException` if there is an error when attempting to save the file. @@ -83,7 +95,10 @@ def tracks2point(tracks, outputFile): """ LOG.info("Writing point shape file: {0}".format(outputFile)) sf = shapefile.Writer(shapefile.POINT) - sf.fields = OBSFIELDS + if netcdf_format: + sf.fields = TCRM_FIELDS + else: + sf.fields = OBSFIELDS LOG.debug("Processing {0} tracks".format(len(tracks))) @@ -101,7 +116,7 @@ def tracks2point(tracks, outputFile): return -def tracks2line(tracks, outputFile, dissolve=False): +def tracks2line(tracks, outputFile, dissolve=False, netcdf_format=False): """ Writes tracks to a shapefile as a collection of line features @@ -118,12 +133,16 @@ def tracks2line(tracks, outputFile, dissolve=False): :type dissolve: boolean :param dissolve: Store track features or track segments. + :param bool netcdf_format: Whether tracks are in TCRM format + :raises: :mod:`shapefile.ShapefileException` if there is an error when attempting to save the file. """ LOG.info("Writing line shape file: {0}".format(outputFile)) sf = shapefile.Writer(shapefile.POLYLINE) - if dissolve: + if netcdf_format: + sf.fields = TCRM_FIELDS + elif dissolve: sf.fields = EVENTFIELDS else: sf.fields = OBSFIELDS @@ -157,19 +176,22 @@ def tracks2line(tracks, outputFile, dissolve=False): sf.line([lines]) - minPressure = track.trackMinPressure - maxWind = track.trackMaxWind + if netcdf_format: + sf.record(*track.data[0]) + else: + minPressure = track.trackMinPressure + maxWind = track.trackMaxWind - age = track.TimeElapsed.max() + age = track.TimeElapsed.max() - startYear = track.Year[0] - startMonth = track.Month[0] - startDay = track.Day[0] - startHour = track.Hour[0] - startMin = track.Minute[0] - record = [track.CycloneNumber[0], startYear, startMonth, startDay, - startHour, startMin, age, minPressure, maxWind] - sf.record(*record) + startYear = track.Year[0] + startMonth = track.Month[0] + startDay = track.Day[0] + startHour = track.Hour[0] + startMin = track.Minute[0] + record = [track.CycloneNumber[0], startYear, startMonth, startDay, + startHour, startMin, age, minPressure, maxWind] + sf.record(*record) else: if len(track.data) == 1: @@ -266,16 +288,22 @@ def tracks2line(tracks, outputFile, dissolve=False): line_output_file = filename + '_line.shp' dissolve_output_file = filename + '_dissolve.shp' - if track_file.endswith("nc"): + if track_file.endswith(".nc"): + from Utilities.track import ncReadTrackData tracks = ncReadTrackData(track_file) - else: + netcdf_format = True + + elif track_file.endswith(".csv"): tracks = loadTrackFile(config_file, track_file, source, calculateWindSpeed=True) + netcdf_format = False + else: + raise ValueError("format of {} is not recognizable".format(track_file)) - tracks2point(tracks, pt_output_file) - tracks2line(tracks, line_output_file) - tracks2line(tracks, dissolve_output_file, dissolve=True) + 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) LOG.info("Completed tracks2shp") diff --git a/docs/install.rst b/docs/install.rst index d3c9ca70..28908cac 100644 --- a/docs/install.rst +++ b/docs/install.rst @@ -3,15 +3,14 @@ Installation ============ -Installing TCRM is intended to be a simple process, requiring only a -small amount of compilation and basic understanding of command line -operations. TCRM has been installed and (lightly) tested on a range of -unix-based systems, Windows and Mac OS/X systems. +Installing TCRM is intended to be a simple process, requiring only basic +understanding of command line operations. TCRM has been installed and (lightly) +tested on a range of unix-based systems, Windows and Mac OS/X systems. .. _downloading: Downloading ----------- +----------- The TCRM code can be downloaded from the `Geoscience Australia GitHub page `_. @@ -28,7 +27,7 @@ TCRM. .. _environment: Setting the environment ----------------------- +----------------------- To enable TCRM to run correctly, you may need to change some environment settings. The important variable to set is the @@ -40,6 +39,7 @@ A complete discussion on environment variables in Python is given in the `Python documentation `_. + Windows ~~~~~~~ @@ -52,7 +52,7 @@ variables on different Windows systems. BASH shell ~~~~~~~~~~ -:: +.. code-block:: bash export PYTHONPATH=$PYTHONPATH:/path/to/tcrm:/path/to/tcrm/Utilities @@ -60,18 +60,15 @@ BASH shell CSH/TCSH shell ~~~~~~~~~~~~~~ -:: +.. code-block:: tcsh setenv PYTHONPATH $PYTHONPATH:/path/to/tcrm:/path/to/tcrm/Utilities - - - .. _dependencies: Dependencies ------------ +------------ TCRM relies on a number of additional libraries that are not part of the standard library. There are several ways to obtain the required @@ -79,7 +76,7 @@ libraries -- using Python's recommended tool `pip `_, installing a distribution such as `Python(x,y) package `_ (for Windows environments) or `Anaconda -`_ (cross-platform), or +`_ (cross-platform), or installing the libraries from source or binary installers (pre-compiled binary Windows installer versions for all the libraries (both 32-bit and 64-bit) can be obtained `here @@ -88,7 +85,7 @@ installing the libraries from source or binary installers For detailed instructions on installation of these dependencies, please see the documentation for each individual library. -* `Python `_ - v2.7 preferred +* `Python `_ - v3.5 or later * `Numpy `_ - v1.6 or later * `Scipy `_ - v0.12 or later * `Matplotlib `_ v1.2 or later. @@ -108,9 +105,11 @@ Using pip If you have `pip `_ installed, the required modules can be installed using the following command, -executed in the main TCRM directory:: +executed in the main TCRM directory + +.. code-block:: bash - pip -v install -r requirements.txt + pip -v install -r requirements.txt This will automatically build the required libraries (listed in the ``requirements.txt`` file) and any dependencies. ``pip`` must be on @@ -118,63 +117,27 @@ the ``$PATH`` for this to work. .. _compilation: -Compiling the extensions ------------------------- - -The model requires a number of C extensions to be compiled before -execution. These can be built using Python's inbuilt :mod:`distutils` -module. - - -Unix -~~~~ -From the base directory, execute the build process:: - python installer/setup.py build_ext -i - -Ubuntu -~~~~~~ -The github branch issue_25 (created from the v2.0 branch) had an environment created by `installing miniconda -`_ and executing the following commands:: - - ~/miniconda2/bin/conda create --name tcrm - ~/miniconda2/bin/source activate tcrm - ~/miniconda2/bin/conda install numpy - ~/miniconda2/bin/conda install scipy - ~/miniconda2/bin/conda install matplotlib - ~/miniconda2/bin/conda install basemap - ~/miniconda2/bin/conda install netcdf4 - ~/miniconda2/bin/conda install shapely - ~/miniconda2/bin/conda install Tornado - ~/miniconda2/bin/conda install statsmodel - ~/miniconda2/bin/conda install seaborn - ~/miniconda2/bin/pip --proxy=http://localhost:3128 install simplejson - - -The following libraries were needed to compile the C extensions, and run the unit tests:: +Using Anaconda +~~~~~~~~~~~~~~ - sudo apt install libgl1-mesa-glx - sudo apt-get install python-numpy-dev +To install ``tcrm``, make a new environment: -The C extensions were compiled from the trcm directory with:: +.. code-block:: bash - (tcrm) user@server:~/tcrm$ python intaller/setup.py build_ext -i + conda env create -f tcrmenv.yml -An error occurred where the include file seems to have changed paths. It may be a one off, -or it may reoccur in another version of Linux. The error was in KPDF.c and the change was to -comment out one line and replace it with another.:: +After creating the environment the user needs to move to that environment using the command - #include "numpy/arrayobject.h" - /* #include "arrayobject.h" */ +.. code-block:: bash -A requiremements file was created in the root directory called ``linux_v20.yml`` and should (it hasn't been tested) -replace the ``conda install`` commands above. The command to use this file is:: + conda activate tcrm - conda env create -f linux_v20.yml +The bash promt will look like -Activating the environment would be:: +.. code-block:: - source activate linux_v20 + (tcrm) user@server:~/tcrm$ Windows @@ -199,7 +162,7 @@ installations. .. _testing: Testing the installation ------------------------ +------------------------ The model code includes a suite of unit tests that ensure elements of the code base will work as expected, even if a user makes @@ -207,7 +170,9 @@ modificaitons to the code. The test suite can be run from the main directory. On Windows, run the ``run_test_all.cmd`` script from the main TCRM directory. On Unix, use -the command:: +the command + +.. code-block:: bash python ./tests/run.py @@ -254,7 +219,9 @@ system. Test the installation ~~~~~~~~~~~~~~~~~~~~~ -Run this command :: +Run this command + +.. code-block:: bash docker run olivierdalang/tcrm nosetests --exe @@ -269,7 +236,9 @@ To run TCRM though Docker, you need to mount a folders containing your inputs and the output folder in the container. This can be done like this (assuming you have a my_conf.ini file in -a folder) :: +a folder) + +.. code-block:: bash docker run -v /path_to/my_data_folder:/home/src/mount -v /path_to/my_output_folder:/home/src/output olivierdalang/tcrm python tcevent.py -v -c mount/my_conf.ini @@ -283,13 +252,17 @@ take some time. Developement ~~~~~~~~~~~~ -You can also use Docker when developping TCRM by mounting the source:: +You can also use Docker when developping TCRM by mounting the source + +.. code-block:: bash git checkout https://github.com/GeoscienceAustralia/tcrm.git cd tcrm docker run -v ${PWD}:/home/src olivierdalang/tcrm python tcevent.py -c example/yasi.ini -If you wish to make changes to the builds steps or dependencies, you need to rebuild the image locally :: +If you wish to make changes to the builds steps or dependencies, you need to rebuild the image locally + +.. code-block:: bash docker build -t olivierdalang/tcrm . @@ -297,10 +270,12 @@ Releases ~~~~~~~~ For users to be able to use the docker image out of the box without having to rebuild it locally, -the image must be pushed to the docker hub repository like this :: +the image must be pushed to the docker hub repository like this + +.. code-block:: bash docker build -t olivierdalang/tcrm . docker login docker push olivierdalang/tcrm -This can be setup to be done automatically after pushes through docker hub. \ No newline at end of file +This can be setup to be done automatically after pushes through docker hub. diff --git a/setup.py b/setup.py index 5543b6e6..88b0ca54 100644 --- a/setup.py +++ b/setup.py @@ -21,7 +21,7 @@ 'numpy', 'scipy', 'pandas', - 'cartopy' + 'cartopy', 'affine', 'matplotlib', 'basemap', @@ -34,7 +34,6 @@ 'seaborn', 'shapely', 'simplejson', - 'sqlite', 'statsmodels', 'tqdm', 'xarray', diff --git a/tests/test_data/windFieldTestData.pkl b/tests/test_data/windFieldTestData.pkl index af0fab75..e831c6f5 100644 Binary files a/tests/test_data/windFieldTestData.pkl and b/tests/test_data/windFieldTestData.pkl differ diff --git a/wind/__init__.py b/wind/__init__.py index fa63900c..a3540943 100644 --- a/wind/__init__.py +++ b/wind/__init__.py @@ -7,7 +7,7 @@ primary vortex of the simulated TC, and bounday layer models that define the asymmetry induced by surface friction and forward motion of the TC over the earth's surface. The final output from the module is a -netCDF file containing the maximum surface gust wind speed (a 10-minute +netCDF file containing the maximum surface gust wind speed (a 10-minute mean wind speed, at 10 metres above ground level), along with the components (eastward and westward) that generated the wind gust and the minimum mean sea level pressure over the lifetime of the event. If multiple @@ -28,17 +28,18 @@ """ -import numpy as np import logging as log import itertools import math import os import sys -import tqdm -from . import windmodels from os.path import join as pjoin from collections import defaultdict +import numpy as np +import tqdm +from . import windmodels + from PlotInterface.maps import saveWindfieldMap from Utilities.files import flModDate, flProgramVersion @@ -48,9 +49,7 @@ from Utilities.parallel import attemptParallel import Utilities.nctools as nctools - from Utilities.track import ncReadTrackData, Track - from ProcessMultipliers import processMultipliers as pM class WindfieldAroundTrack(object): @@ -521,7 +520,7 @@ def dumpGustsFromTracks(self, trackiter, windfieldPath, dumpfile) if self.multipliers is not None: - self.calcLocalWindfield(results_mult) + self.calcLocalWindfield(results) #self.plotGustToFile((lat, lon, gust, Vx, Vy, P), plotfile) def plotGustToFile(self, result, filename): @@ -702,6 +701,9 @@ def dumpGustsFromTrackfiles(self, trackfiles, windfieldPath, def calcLocalWindfield(self, results): """ + Calculate a local wind field using the regional windfield data + + :param results: collection of :tuple: track and wind field data """ # Load a multiplier file to determine the projection: @@ -721,7 +723,7 @@ def calcLocalWindfield(self, results): 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 @@ -777,7 +779,7 @@ def loadTracksFromPath(path): """ files = os.listdir(path) trackfiles = [pjoin(path, f) for f in files if f.startswith('tracks')] - msg = 'Processing %d track files in %s' % (len(trackfiles), path) + msg = 'Processing {0} track files in {1}'.format(len(trackfiles), path) log.info(msg) return loadTracksFromFiles(sorted(trackfiles)) @@ -805,7 +807,7 @@ def run(configFile, callback=None): """ log.info('Loading wind field calculation settings') - + # Get configuration config = ConfigParser() @@ -840,11 +842,11 @@ def run(configFile, callback=None): timestepCallback = ts.extract else: timestepCallback = None - + multipliers = None - if config.has_option('Input','Multipliers'): + if config.has_option('Input', 'Multipliers'): multipliers = config.get('Input', 'Multipliers') - + thetaMax = math.radians(thetaMax) # Attempt to start the track generator in parallel @@ -868,7 +870,7 @@ def run(configFile, callback=None): multipliers=multipliers, windfieldPath=windfieldPath) - log.info('Dumping gusts to %s' % windfieldPath) + log.info('Dumping gusts to {0}'.format(windfieldPath)) # Get the trackfile names and count @@ -876,7 +878,7 @@ def run(configFile, callback=None): trackfiles = [pjoin(trackPath, f) for f in files if f.startswith('tracks')] nfiles = len(trackfiles) - log.info('Processing %d track files in %s' % (nfiles, trackPath)) + log.info('Processing {0} track files in {1}'.format(nfiles, trackPath)) # Do the work diff --git a/wind/windmodels.py b/wind/windmodels.py index 88158030..bde28052 100644 --- a/wind/windmodels.py +++ b/wind/windmodels.py @@ -275,8 +275,9 @@ def vorticity(self, R): """ rr = R * 1000. - Z = (np.sign(self.f) * 2 * self.vMax * self.rMax / (self.rMax - ** 2 + rr ** 2) + np.sign(self.f) * 2 * self.vMax * + 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) return Z @@ -330,7 +331,8 @@ def secondDerivative(self): try: assert d2Vm < 0.0 except AssertionError: - log.critical("Pressure deficit: {0:.2f} hPa, RMW: {1:%2f} km".format(dP/100., rMax/1000.)) + log.critical(("Pressure deficit: {0:.2f} hPa," + " RMW: {1:%2f} km".format(dP/100., rMax/1000.))) raise return d2Vm @@ -349,7 +351,7 @@ def firstDerivative(self): 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))) + (2*(4*beta*dP/rho + E*(f*rMax)**2))) return dVm def velocity(self, R): @@ -374,9 +376,9 @@ def velocity(self, R): delta = (self.rMax / R) ** self.beta edelta = np.exp(-delta) - V = (np.sqrt((self.dP * self.beta / self.rho) - * delta * edelta + (R * self.f / 2.) ** 2) - R * - np.abs(self.f) / 2.) + V = (np.sqrt((self.dP * self.beta / self.rho) * + delta * edelta + (R * self.f / 2.) ** 2) - + R * np.abs(self.f) / 2.) icore = np.where(R <= self.rMax) V[icore] = (R[icore] * (R[icore] * (R[icore] * aa + bb) + cc)) @@ -387,7 +389,7 @@ def vorticity(self, R): """ Calculate the vorticity associated with the (gradient level) vortex. - + :param R: :class:`numpy.ndarray` of distance of grid from the TC centre (metres). @@ -708,9 +710,6 @@ def vorticity(self, R): """ - rm = self.rMax * 1000. - rm2 = self.rMax2 * 1000. - rr = R*1000. # Scale dp2 if dP is less than 1500 Pa: if self.dP < 1500.: dp2 = (self.dP / 1500.) * (800. + (self.dP - 800.) / 2000.) @@ -1051,7 +1050,7 @@ def field(self, R, lam, vFm, thetaFm, thetaMax=0.): Cd = 0.002 # Constant drag coefficient Vm = np.abs(V).max() if (vFm > 0) and (Vm/vFm < 5.): - Umod = vFm * np.abs(1. - (vFm/Vm)/vFm) + Umod = vFm * np.abs(1.25*(1. - (vFm/Vm))) else: Umod = vFm Vt = Umod * np.ones(V.shape)