import numpy as np
import datetime
import xarray as xr
[docs]def display_message():
print('########################################################################')
print('This is a close source toolbox\n To get more information please email:\n\
r.zyngfogel@calypso.science or\n\
b.beamsley@metocean.co.nz')
print('########################################################################')
import numpy as np
[docs]def lanc(numwt, haf):
"""
Generates a numwt + 1 + numwt lanczos cosine low pass filter with -6dB
(1/4 power, 1/2 amplitude) point at haf
Function from Oceans Toolbox
Parameters
----------
numwt : int
number of points
haf : float
frequency (in 'cpi' of -6dB point, 'cpi' is cycles per interval.
For hourly data cpi is cph,
Examples
--------
>>> from oceans.filters import lanc
>>> import matplotlib.pyplot as plt
>>> t = np.arange(500) # Time in hours.
>>> h = 2.5 * np.sin(2 * np.pi * t / 12.42)
>>> h += 1.5 * np.sin(2 * np.pi * t / 12.0)
>>> h += 0.3 * np.random.randn(len(t))
>>> wt = lanc(96+1+96, 1./40)
>>> low = np.convolve(wt, h, mode='same')
>>> high = h - low
>>> fig, (ax0, ax1) = plt.subplots(nrows=2)
>>> _ = ax0.plot(high, label='high')
>>> _ = ax1.plot(low, label='low')
>>> _ = ax0.legend(numpoints=1)
>>> _ = ax1.legend(numpoints=1)
"""
summ = 0
numwt += 1
wt = np.zeros(numwt)
# Filter weights.
ii = np.arange(numwt)
wt = 0.5 * (1.0 + np.cos(np.pi * ii * 1.0 / numwt))
ii = np.arange(1, numwt)
xx = np.pi * 2 * haf * ii
wt[1 : numwt + 1] = wt[1 : numwt + 1] * np.sin(xx) / xx
summ = wt[1 : numwt + 1].sum()
xx = wt.sum() + summ
wt /= xx
return np.r_[wt[::-1], wt[1 : numwt + 1]]
[docs]def qkhfs( w, h ):
"""
Quick iterative calculation of kh in gravity-wave dispersion relationship
kh = qkhfs(w, h )
Input
w - angular wave frequency = 2*pi/T where T = wave period [1/s]
h - water depth [m]
Returns
kh - wavenumber * depth [ ]
Orbital velocities from kh are accurate to 3e-12 !
RL Soulsby (2006) "Simplified calculation of wave orbital velocities"
HR Wallingford Report TR 155, February 2006
Eqns. 12a - 14
from https://github.com/csherwood-usgs/crspy/blob/master/crspy.py
"""
g = 9.81
x = w**2.0 *h/g
y = np.sqrt(x) * (x<1.) + x *(x>=1.)
# is this faster than a loop?
t = np.tanh( y )
y = y-( (y*t -x)/(t+y*(1.0-t**2.0)))
t = np.tanh( y )
y = y-( (y*t -x)/(t+y*(1.0-t**2.0)))
t = np.tanh( y )
y = y-( (y*t -x)/(t+y*(1.0-t**2.0)))
kh = y
return kh
[docs]def peaks(y):
""" return peaks and trough indx
from Nagi Hatoum peaks.m
copyright 2005"""
ds=np.diff(y)
ds = np.insert(ds, 0, ds[0], axis=0) #pad diff
fil=np.nonzero((ds[1:]==0))[0]+1 #find zeros
ds[fil]=ds[fil-1] #replace zeros
ds=np.sign(ds)
ds=np.diff(ds)
t=np.nonzero((ds>0))[0]
p=np.nonzero((ds<0))[0]
return p,t
[docs]def PolyArea(x,y):
bad=np.logical_or(np.isnan(x),np.isnan(y))
x=x[~bad]
y=y[~bad]
return 0.5*np.abs(np.dot(x,np.roll(y,1))-np.dot(y,np.roll(x,1)))
[docs]def do_occurence(dpm,min_occ):
dir_int=dir_interval(45,'centred')
dir_int_name=np.array(['N','NE','E','SE','S','SW','W','NW'])
occ=np.ones((len(dir_int)-1))
for j in range(0,len(dir_int)-1):
if dir_int[j+1] <= dir_int[j]:
D=((np.mod(dpm,360)>dir_int[j]) | (np.mod(dpm,360)<=dir_int[j+1])).nonzero()[0]
else:
D=((np.mod(dpm,360)>dir_int[j]) & (np.mod(dpm,360)<=dir_int[j+1])).nonzero()[0]
occ[j]=(len(D)/len(dpm[~np.isnan(dpm)]))*100;
Occ=dir_int_name[np.where(occ>=min_occ)]
return Occ
[docs]def dyadlength(x):
'''% dyadlength -- Find length and dyadic length of array
% Usage
% [n,J] = dyadlength(x)
% Inputs
% x array of length n = 2^J (hopefully)
% Outputs
% n length(x)
% J least power of two greater than n
%
% Side Effects
% A warning is issued if n is not a power of 2.
%
% See Also
% quadlength, dyad, dyad2ix
%'''
n = len(x)
J = np.ceil(np.log(n)/np.log(2));
if 2**J != n:
print('Warning in dyadlength: n != 2^J')
'''%
% Part of Wavelab Version 850
% Built Tue Jan 3 13:20:40 EST 2006
% This is Copyrighted Material
% For Copying permissions see COPYING.m
% Comments? e-mail wavelab@stat.stanford.edu '''
return n,J
[docs]def spdir2uv(spd,direc,origin='going to'):
ang_rot = 180 if origin=='coming from' else 0
direcR = np.deg2rad(direc + ang_rot)
u = spd * np.sin(direcR)
v = spd * np.cos(direcR)
return u,v
[docs]def uv2spdir(u,v,origin='going to'):
ang_rot = 180 if origin=='coming from' else 0
vetor = u + v * 1j
mag = np.abs(vetor)
direc = xr.ufuncs.angle(vetor, deg=True) + ang_rot
direc = np.mod(90 - direc, 360)
return mag,direc
[docs]def get_opt(var,label_name,default):
if hasattr(var,label_name):
label=getattr(var,label_name)
else:
label=''
if label=='':
label=default
return label
[docs]def get_number_of_loops(time_blocking):
if 'annual' in time_blocking.lower():
number_of_loops=1
identidifers=['Annual']
month_identidier=[np.arange(1,13,1)]
elif 'seasonal' in time_blocking.lower():
identidifers,month_identidier=get_seasons(time_blocking)
number_of_loops=4+1
month_identidier.append(np.arange(1,13,1))
identidifers.append('Annual')
elif 'monthly' in time_blocking.lower():
identidifers=[]
month_identidier=[]
for m in range(1,13):
identidifers.append(datetime.date(1900, m, 1).strftime('%B'))
month_identidier.append([m])
identidifers.append('Annual')
number_of_loops=12+1
month_identidier.append(np.arange(1,13,1))
return number_of_loops,identidifers,month_identidier
[docs]def get_seasons(typ):
seasons_name=['Summer','Autumn','Winter','Spring']
if 'south' in typ.lower():
seasons=[[12,1,2],[3,4,5],[6,7,8],[9,10,11]]
else:
seasons=[[6,7,8],[9,10,11],[12,1,2],[3,4,5]]
return seasons_name,seasons
[docs]def get_increment(data,s):
if len(s)>3:
interval=s
else:
if len(s)==3:
interval=np.arange(s[0],s[2]+s[1],s[1])
else:
upper_limit=np.ceil(np.max(data/s[1]))*s[1]
interval=np.arange(s[0],upper_limit+s[1],s[1])
return interval
[docs]def wavenuma(ang_freq, water_depth):
"""Chen and Thomson wavenumber approximation."""
k0h = 0.10194 * ang_freq * ang_freq * water_depth
D = [0, 0.6522, 0.4622, 0, 0.0864, 0.0675]
a = 1.0
for i in range(1, 6):
a += D[i] * k0h ** i
return (k0h * (1 + 1.0 / (k0h * a)) ** 0.5) / water_depth
[docs]def cart2sph(x,y,z):
azimuth = np.arctan2(y,x)
elevation = np.arctan2(z,np.sqrt(x**2 + y**2))
r = np.sqrt(x**2 + y**2 + z**2)
return azimuth, elevation, r
[docs]def sph2cart(azimuth,elevation,r):
x = r * np.cos(elevation) * np.cos(azimuth)
y = r * np.cos(elevation) * np.sin(azimuth)
z = r * np.sin(elevation)
return x, y, z
[docs]def cart2pol(x,y):
th = np.arctan2(y,x)
r = np.sqrt(x**2 + y**2)
return th,r
[docs]def dir_interval(dir_swath=45,mode='centred'):
#Calculates directional intervals based on dir_swath which must
#be a whole number divisor of 360.
#mode can be either "centred" or "not-centred"
#Default is centred
# here you specify the directional swath that you wish to analyse
if dir_swath==0:
dir_swath=360;
if 360 % dir_swath !=0:
print('The chosen interval is not a multiple of 360')
raise
interval=[]
if mode == 'centred' or mode=='centered':
for i in range(0,int(360//dir_swath)):
if i==0:
interval.append(360-(dir_swath/2))
elif i==1:
interval.append(dir_swath/2);
elif i>1:
interval.append(interval[i-1]+dir_swath)
interval.append(interval[i]+dir_swath-1e-10) # this is a fudge factor (-1e-10) but is a good soln.
else: # 'not-centred'
for i in range(0,int(360//dir_swath)):
if i==0:
interval.append(0)
else:
interval.append(interval[i-1]+dir_swath)
interval.append(interval[i]+dir_swath)
return interval
[docs]def degToCompass(num):
if isinstance(num,list) or isinstance(num,np.array):
if num[0]==0 and num[1]==360:
return 'Omni'
rads = np.deg2rad(num)
av_sin = np.mean(np.sin(rads))
av_cos = np.mean(np.cos(rads))
ang_rad = np.arctan2(av_sin,av_cos)
num = np.mod(np.round(np.rad2deg(ang_rad),1),360)
val=int((num/22.5)+.5)
arr=["N","NNE","NE","ENE","E","ESE", "SE", "SSE","S","SSW","SW","WSW","W","WNW","NW","NNW"]
return arr[(val % 16)]
if __name__ == '__main__':
#print(dir_interval(dir_swath=22.5,mode='centred'))
#print(dir_interval(dir_swath=45,mode='not-centred'))
print(degToCompass([350,3,10]))