Commit 3a8230dd authored by Carl Schreck's avatar Carl Schreck
Browse files

Automated Nightly Commit - Sat Jun 8 00:00:21 EDT 2019

parent 3eb0face
""" Library of Carl's python functions. """
__author__ = "Carl Schreck"
__email__ = "cjschrec@ncsu.edu"
__copyright__ = "Copyright 2019, North Carolina State University"
__license__ = "BSD-3.0"
import datetime
def tstamp(message='',maxcols=55):
""" Print a string with a timestamp. """
print( str(message).ljust(maxcols),
datetime.datetime.today().strftime(' | %F %T'))
""" Library of functions for accessing IBTrACS data """
__author__ = "Carl Schreck"
__email__ = "cjschrec@ncsu.edu"
__copyright__ = "Copyright 2019, North Carolina State University"
__license__ = "BSD-3.0"
import xarray as xr
import numpy as np
import sys
#from shapely.geometry import Point, Polygon
sys.path.append('/home/carl/lib/python')
import cjs
class Ibtracs:
''' This class uses IBTrACS data. '''
def __init__(self, inpath='since1980'):
epochs = [ 'last3years', 'ALL', 'since1980' ]
if inpath in epochs:
self.inpath = ( '~/data/ibtracs/v04r00/IBTrACS.' + inpath
+ '.v04r00.nc' )
else:
self.inpath = inpath
self.ds = xr.open_dataset(self.inpath)
self.ds['lon'] = self.ds.lon.where(self.ds.lon > 0, self.ds.lon+360)
def find_tropical(self, thresh=35, keep_subtropical=True):
if keep_subtropical:
bad_natures = [b'DS', b'ET']
else:
bad_natures = [b'SS', b'DS', b'ET']
is_tropical = ~self.ds.nature.isin(bad_natures)
if thresh < 0:
retval = is_tropical
else:
is_strong_enough = self.ds.usa_wind >= thresh
retval = is_tropical & is_strong_enough
return(retval)
def find_first_ind(self, thresh=35, keep_subtropical=True):
is_main = self.ds.track_type != b'spur'
is_tropical = self.find_tropical(thresh, keep_subtropical)
if thresh >= 0:
is_strong = self.ds.usa_wind > thresh
else:
lmi = np.amax(self.ds.usa_wind, axis=1)
is_strong = self.ds.usa_wind == lmi
nstrong = np.sum(is_strong, axis=1)
is_eligible = is_main & is_tropical & is_strong
neligible = np.sum(is_eligible, axis=1)
ind_eligible = np.argmax(is_eligible,1)
retval = False & is_eligible
for i,j in enumerate(ind_eligible):
if neligible[i] > 0:
retval[i,j] = True
retval['oned'] = np.any(retval, axis=1)
# print(np.c_[self.ds.sid.values[retval1d],self.ds.usa_wind.values[retval],lmi[retval1d]])
return(retval)
def check_basin(self, target_basin='NA'):
target_basin = target_basin.upper()
# Handle special cases
if target_basin in ['ALL', 'GL']:
retval = self.ds.lat.notnull()
elif target_basin is 'NH':
retval = self.ds.lat >= 0
elif target_basin is 'SH':
retval = self.ds.lat < 0
elif target_basin is 'WMOSP':
retval = (self.ds.lat < 0) & (self.ds.lon >= 160)
elif target_basin is 'WMOSI':
retval = (self.ds.lat < 0) & (self.ds.lon < 90)
elif target_basin is 'WMOAUS':
retval = ((self.ds.lat < 0) & (self.ds.lon >= 90)
& (self.ds.lon < 160))
elif target_basin is 'WMOCP':
retval = ((self.ds.lat >= 0) & (self.ds.lat < 40)
& (self.ds.lon >= 180) & (self.ds.lon < 220))
elif target_basin is 'NA_MDR':
retval = ((self.ds.lat >= 10) & (self.ds.lat < 20)
& (self.ds.lon >= 275) & (self.ds.lon < 340))
elif target_basin is 'US_KOSSIN':
retval = ((self.ds.lat >= 25) & (self.ds.lat < 50)
& (self.ds.lon >= 260) & (self.ds.lon < 295))
elif target_basin is 'US_GC':
retval = ((self.ds.lat >= 25) & (self.ds.lat < 35)
& (self.ds.lon >= 260) & (self.ds.lon < 275))
elif target_basin is 'US_FL':
retval = ((self.ds.lat >= 20) & (self.ds.lat < 30)
& (self.ds.lon >= 275) & (self.ds.lon < 285))
#US_EC will require point in polygon
#manual function: http://www.ariel.com.au/a/python-point-int-poly.html
#https://automating-gis-processes.github.io/CSC18/lessons/L4/point-in-polygon.html
else
retval = target_basin == self.ds.basin
return(retval)
'''
NAME
Custom Colormaps for Matplotlib
PURPOSE
This program shows how to implement make_cmap which is a function that
generates a colorbar. If you want to look at different color schemes,
check out https://kuler.adobe.com/create.
PROGRAMMER(S)
Chris Slocum (original); Larry Gloeckler (revised 20141010)
REVISION HISTORY
20130411 -- Initial version created
20140313 -- Small changes made and code posted online
20140320 -- Added the ability to set the position of each color
20141010 -- Added a block that converts color map arrays (e.g., color maps
generated in MATLAB) into lists of tuples
'''
def mplCmap(colors, position=None, bit=False):
'''
mplCmap takes a list of tuples which contain RGB values. The RGB
values may either be in 8-bit [0 to 255] (in which bit must be set to
True when called) or arithmetic [0 to 1] (default). mplCmap returns
a cmap with equally spaced colors.
Arrange your tuples so that the first color is the lowest value for the
colorbar and the last is the highest.
position contains values from 0 to 1 to dictate the location of each color.
'''
import sys
import matplotlib as mpl
import numpy as np
# Convert array to list of tuples (if applicable)
if isinstance(colors, tuple):
colors = list(map(tuple, colors))
# Evenly space numbers from 0 to 255
bit_rgb = np.linspace(0,1,256)
if position == None:
position = np.linspace(0,1,len(colors))
else:
if len(position) != len(colors):
sys.exit("position length must be the same as colors")
elif position[0] != 0 or position[-1] != 1:
sys.exit("position must start with 0 and end with 1")
# Convert bit values to arithmetic values
if bit:
for i in range(len(colors)):
colors[i] = (bit_rgb[colors[i][0]],
bit_rgb[colors[i][1]],
bit_rgb[colors[i][2]])
# Initialize RGB dictionary
cdict = {'red':[], 'green':[], 'blue':[]}
# Assign values to red, green, and blue positions in cdict
for pos, color in zip(position, colors):
cdict['red'].append((pos, color[0], color[0]))
cdict['green'].append((pos, color[1], color[1]))
cdict['blue'].append((pos, color[2], color[2]))
cmap = mpl.colors.LinearSegmentedColormap('my_colormap',cdict,256)
return cmap
\ No newline at end of file
"""
myMethods.py
Contains several customized frequently-run Python methods
"""
import numpy as np
import numpy.linalg as lin
from matplotlib import mlab
import h5py as h5
import scipy.io as sio
import scipy.stats as stat
from scipy.ndimage.filters import minimum_filter, maximum_filter
import datetime as dt
#from spharm import gaussian_lats_wts, Spharmt,regrid
import matplotlib.pyplot as plt
from matplotlib.patches import Polygon
# A method to read data from a text or similar file
def readData(filename, startLine=0):
# Read in file
fileobj = open(filename, 'r')
outputstr = fileobj.readlines()
fileobj.close()
# Split strings into lists
for i in xrange(len(outputstr)):
outputstr[i] = outputstr[i].split(' ')
# Output to array with type double precision
# outputarray = outputstr
outputarray = np.array(outputstr[startLine:len(outputstr)])
outputarray = outputarray.astype('d')
return outputarray
# A method to load Matlab data
def loadMat(filename, var):
"""
Uses two different approaches to load HDF (.mat) files.
This ensures that HDF files of different versions can
be loaded using this method.
"""
##Initialize dictionary for temporary variable storage
#mdict = {}
try:
fn = sio.loadmat(filename)
except NotImplementedError:
fn = h5.File(filename, 'r')
data = fn.get(var)
data = np.array(data, dtype='d')
return data
# A method to covert regrid gaussian data
def gauss2grid(data,nlats_reg,nlons_reg,nlats_gau,nlons_gau):
""" This function uses spharm to regrid data from gaussian grid to lat/lon
grid. Input data must be 2-d (Lat x Lon).
data: data to convert (lat x lon)
nlats_reg: number of lats
nlons_reg: number of lons
nlats_gau: number of gaussian lats
nlons_gau: number of gaussian lons
"""
gaulats, wts = gaussian_lats_wts(nlats_gau)
# create Spharmt instances for regular and gaussian grids.
reggrid = Spharmt(nlons_reg,nlats_reg,gridtype='regular')
gaugrid = Spharmt(nlons_gau,nlats_gau,gridtype='gaussian')
# regrid from gaussian to regular grid.
data_regrid = regrid(gaugrid,reggrid,data)
return data_regrid
# A method to return a date list
def dateRange(start, end):
"""
dateRange returns a dictionary of arrays with keys 'day', 'month',
and 'year'. To access elements in these arrays, use array_name['key'].
Inputs:
:start: A tuple defining a start date (tuple; e.g., (1979,1,1))
:end: A tuple defining an end date (tuple; e.g., (2010,12,31))
"""
# Check for class instance of input dates
if isinstance(start[0], float):
startTemp = np.zeros_like(start, dtype='int')
for i in xrange(len(startTemp)):
startTemp[i] = start[i]
start = tuple(startTemp)
if isinstance(end[0], float):
endTemp = np.zeros_like(end, dtype='int')
for i in xrange(len(endTemp)):
endTemp[i] = end[i]
end = tuple(endTemp)
# Initialize dictionary
dateList = {}
# Fill dictionary with dates
start = dt.datetime(*start[:6]) # * is an argument-unpacking operator
end = dt.datetime(*end[:6])
n = (end + dt.timedelta(days=1) - start).days
dateList['year'] = np.array([(start + dt.timedelta(days=i)).year for i in xrange(n)])
dateList['month'] = np.array([(start + dt.timedelta(days=i)).month for i in xrange(n)])
dateList['day'] = np.array([(start + dt.timedelta(days=i)).day for i in xrange(n)])
return dateList
# A method to trim arrays
def compressArray(inArray, arrayStart, arrayEnd, trimStart, trimEnd, ax=0):
"""
compressArray returns a compressed array corresponding the the input array.
The input array is compressed along the nth axis (default is 0th).
Inputs:
:array: The array to be compressed (array)
:arrayStart: The start date of the uncompressed array (tuple; e.g., (1974,6,1))
:arrayEnd: The end date of the uncompressed array (tuple; e.g., (2013,12,31))
:trimStart: The start date of the compressed array (tuple; e.g., (1979,1,1))
:trimEnd: The end date of the compressed array (tuple; e.g., (2010,12,31))
:ax: The axis over which the compression occurs (int; default=0)
"""
# Generate date list
dateList = dateRange(arrayStart, arrayEnd)
# Generate index vector with length equal to the length of the date list
arrayIndex = np.array(range(0, len(dateList['year'])))
# Lower and upper bounds of compression vector
lowerBound = np.where(np.logical_and(np.logical_and(dateList['year'] == trimStart[0],
dateList['month'] == trimStart[1]), dateList['day'] == trimStart[2]))
upperBound = np.where(np.logical_and(np.logical_and(dateList['year'] == trimEnd[0],
dateList['month'] == trimEnd[1]), dateList['day'] == trimEnd[2]))
# Condition for which to compress the array
condition = np.squeeze(np.logical_and(arrayIndex >= lowerBound,
arrayIndex <= upperBound))
# Compress the array
outArray = np.compress(condition, inArray, ax)
return outArray
# A method to create seasonal cycles of data
def createCycles(days, harmonics):
# Check for NaNs --- to come in future version
# condition = np.isnan(data)
t = np.arange(days)
cycle = np.zeros((days, harmonics*2))
cycle = np.array(cycle, dtype='d')
j = 0
for i in xrange(1, harmonics+1):
cycle[:,j] = np.sin(i*2*np.pi*t/365.25)
cycle[:,j+1] = np.cos(i*2*np.pi*t/365.25)
j = j + 2
return cycle
# A method to create data anomalies
def createAnomalies(data, harmonics=4, removetrend=False, standardized=False):
"""
Creates anomalies over a specificed climatology by removing N harmonics of
the seasonal cycle, and regresses the full dataset onto the seasonal cycle.
Inputs:
:data: The data used to generate anomalies
:harmonics: The number of harmonics of the seasonal cycle to be removed. Default is 4.
:removetrend: Remove the long-term trend? Can be 1 or 0, or True or False. Default is False.
:standardized: Standardize the anomalies? Can be 1 or 0, or True or False. Default is False.
"""
# Preallocate anomaly array
#anomaly = N.zeros_like(data)
[le, wi, ht] = np.shape(data)
# Create seasonal cycle
cycle = createCycles(le, harmonics)
# Create array X (column of ones) and insert into cycle array
X = [1]*le
X = np.insert(cycle, 0, X, axis=1)
Y = np.zeros_like(data)
# Create anomalies
for i in xrange(ht):
"""
In Matlab, if A is an m-by-n matrix with m~=n and B is a column vector
with m components, or a matrix with several such columns, the X = A\B is
the solution in the least squares sense to the under- or overdetermined
system of equations AX = B. Thus, numpy.linalg.lstsq is the numpy equivalent
to the left divide operator in Matlab.
"""
# Equivalent to A = X \ squeeze(data(:,:,i))
A = lin.inv(np.dot(X.conj().T, X))
A = np.dot(A, X.conj().T)
A = np.dot(A, data[:,:,i])
Y[:,:,i] = np.dot(X,A)
# Fill the anomaly array
anomaly = data - Y
# Remove linear trend
if removetrend:
for lat in xrange(wi):
for lon in xrange(ht):
anomaly[:,lat,lon] = mlab.detrend_linear(anomaly[:,lat,lon])
# Create standardized anomalies
if standardized:
anomRep = anomaly # Reassign anomalies to anomRep
anomSq = np.power(anomaly, 2) # Squared anomalies
X = [1]*le
X = np.insert(cycle, 0, X, axis=1)
for i in xrange(ht):
A = lin.inv(np.dot(X.conj().T, X)) # Solve for coeffs.
A = np.dot(A, X.conj().T)
A = np.dot(A, anomSq[:,:,i])
Y[:,:,i] = np.dot(X,A) # Solve for variance of seasonal cycle
# Fill anomaly matrix
anomaly[:,:,i] = anomRep[:,:,i]/np.sqrt(Y) # Standardize anomalies
return anomaly
# A method to find extrema (highs/lows) on a plan view map
def extrema(data, mode='wrap', window=10):
# Find the indices of local extrema in the input array
mn = minimum_filter(data, size=window, mode=mode)
mx = maximum_filter(data, size=window, mode=mode)
# If pixel is equal to local max, then data == mx
# If pixel is equal to local min, then data == mn
return np.nonzero(data == mn), np.nonzero(data == mx)
# A method to compute divergence of the horizontal wind on spherical coordinates
def divergence(lats, lons, u, v):
# Constants
Re = 6.37E6
nTime = np.shape(u)[0]
latGrid = np.arange(np.squeeze(len(lats)))
lonGrid = np.arange(np.squeeze(len(lons)))
# Format arrays for derivatives
lam = np.radians(lons) # Lon vector in radians
phi = np.radians(lats) # Lat vector in radians
L, P = np.meshgrid(lam, phi) # Initialize radian grid
P2 = np.kron(np.ones((nTime, 1, 1)), P) # Equivalent to 'repmat' in Matlab
# Compute grid resolution and convert to radians
latRes = np.abs(lats[1] - lats[0])
lonRes = np.abs(lons[1] - lons[0])
dLam = np.radians(lonRes)
dPhi = np.radians(latRes)
# Compute divergence of the horizontal wind using centered differences
d = np.empty_like(P2) * np.nan
d[:, 1:-1, 1:-1] = 1 / Re * (np.cos(P2[:, 1:-1, 1:-1])) * \
(u[:, latGrid[0]+1:latGrid[-1], lonGrid[0]+2:lonGrid[-1]+1] - \
u[:, latGrid[0]+1:latGrid[-1], lonGrid[0]:lonGrid[-1]-1]) / \
(2 * dLam) + 1 / Re * (v[:, latGrid[0]:latGrid[-1]-1, lonGrid[0]+1:lonGrid[-1]] - \
v[:, latGrid[0]+2:latGrid[-1]+1, lonGrid[0]+1:lonGrid[-1]]) / (2 * dPhi) - \
v[:, latGrid[0]+1:latGrid[-1], lonGrid[0]+1:lonGrid[-1]] * np.tan(P2[:, 1:-1, 1:-1]) / Re
return d
# A method to iterate over randomDraw N times with a 1-D array
def bootstrap(indexArray, N=1):
"""
This method draws random elements equal to the number of elements in
indexArray and replaces them in indexArray before each subsequent
draw. This process is repeated N times (1 by default), and output is stored
in a len(indexArray)-N dimensional array.
"""
# Call randomDraw and iterate over method N times
if N == 1:
randArray = np.random.choice(indexArray, len(indexArray))
else:
randArray = np.random.choice(indexArray, (len(indexArray), N))
return randArray
# A method to test statistical significance using a Monte Carlo test
# Currently under construction
def monteCarlo(data, compTime, randTime, nSample, alpha):
"""
Returns results of a Monte Carlo (bootstrap resampling) test.
"""
if data.ndim == 3:
nLat = np.shape(data)[1]
nLon = np.shape(data)[2]
map = 1
elif data.ndim == 2:
nLat = 1
nLon = nLon = np.shape(data)[1]
map = 0
nTime = len(np.squeeze(randTime))
nSize = len(np.squeeze(compTime))
# Initialize arrays
result = np.zeros((nLat, nLon))
timeRand = np.zeros((nSize, nSample))
timeRep = np.zeros((nSize, nSample))
for r in np.arange(nSample):
timeRand[:,r] = np.random.choice(randTime, nSize)
timeRep[:,r] = np.random.choice(compTime, nSize)
for i in np.arange(nLat):
if map:
X = np.squeeze(data[:,i,:])
elif not map:
X = data
Null = np.nanmean(np.reshape(X[timeRand.astype('i'),:], (nSize, nSample, nLon)), 0)
Test = np.nanmean(np.reshape(X[timeRep.astype('i'),:], (nSize, nSample, nLon)), 0)
Null = np.sort(np.squeeze(Null), 0)
Test = np.sort(np.squeeze(Test), 0)
del X
left = np.nanmin(Test, 0) <= np.nanmin(Null, 0)
left = left.astype('i')
xPass = Null[np.floor(nSample * alpha), left]
nPass = np.sum(Test[:,left] < np.ones((nSample, 1)) * xPass, 0)
result[i, left] = nPass >= np.round(nSample * (1 - alpha))
xPass = Null[np.ceil(nSample * (1-alpha)), left.any()]
nPass = np.Sum(Test[:, left.any()] > np.ones((nSample, 1)) * xPass, 0)
result[i, left.any()] = nPass >= np.round(nSample * (1 - alpha))
return result
## A method for testing field significance
def fieldSig(field, stagnantArray, variedArray, N=1000):
"""
fieldSig tests field significance of two different fields. The null
hypothesis is that the two fields are as spatially correlated as two
randomly generated fields. Output takes the form of a 1-D array of length N,
where N is the number of randomly generated composites used to test field
significance.
Inputs:
:stagnantArray: The index array that will remain unchanged through the test
:variedArray: The index array that will be subject to random draws
:field: The data field used to test field significance (e.g., 300-hPa Z)
:N: The number of random composites generated to test field significance (default is 1000)
"""
# Add sqrt of cosine of latitude standardization
# Apply bootstrap to index vector to generate sets of composites
randArray = bootstrap(variedArray, N)
# Initialize correlation array
r = np.zeros((N, 1))
# Composite random and stagnant fields
stagFieldComp = np.nanmean(field[stagnantArray,:,:], axis=0, dtype='d')
fieldComp1 = np.ravel(stagFieldComp) # Ravel into 1-D array
for n in xrange(N):
randFieldComp = np.nanmean(field[randArray[:,n],:,:], axis=0, dtype='d')
fieldComp2 = np.ravel(randFieldComp) # Ravel into 1-D array
# Compute the population correlation between the two fields
# Analogous to using np.corrcoef
c = np.cov(fieldComp1, fieldComp2)
d = np.diag(c)
rTemp = c/np.sqrt(np.multiply.outer(d, d))
r[n] = rTemp[0,1]
return sorted(r)