diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..0f82a58 --- /dev/null +++ b/.gitignore @@ -0,0 +1,8 @@ +# pip installer +*.egg-info + +# vim +.*.swp + +# python +*.pyc diff --git a/.travis.yml b/.travis.yml new file mode 100644 index 0000000..70a172d --- /dev/null +++ b/.travis.yml @@ -0,0 +1,18 @@ +language: python +python: + - 2.7 + - 3.2 +matrix: + fast_finish: true + allow_failures: + - python: 3.2 +virtualenv: + system_site_packages: true +before_install: + - sudo apt-get install -qq python-numpy python-scipy python3-numpy python3-scipy +install: + - LAPACK=/usr/lib/liblapack.so ATLAS=/usr/lib/libatlas.so BLAS=/usr/lib/libblas.so pip install -r requirements_for_travis.txt +script: + - ./run_tests.sh +notifications: + email: false diff --git a/README.md b/README.md index b62cda0..6f52c86 100644 --- a/README.md +++ b/README.md @@ -1,5 +1,6 @@ pytides ======= +[![Build Status](https://travis-ci.org/sam-cox/pytides.svg)](https://travis-ci.org/sam-cox/pytides) ####About Pytides is small Python package for the analysis and prediction of tides. Pytides can be used to extrapolate the tidal behaviour at a given location from its previous behaviour. The method used is that of harmonic constituents, in particular as presented by P. Schureman in Special Publication 98. The fitting of amplitudes and phases is handled by Scipy's leastsq minimisation function. Pytides currently supports the constituents used by NOAA, with plans to add more constituent sets. It is therefore possible to use the amplitudes and phases published by NOAA directly, without the need to perform the analysis again (although there may be slight discrepancies for some constituents). diff --git a/pytides/astro.py b/pytides/astro.py index b9d8ecf..9c5b730 100644 --- a/pytides/astro.py +++ b/pytides/astro.py @@ -8,201 +8,221 @@ # to make a significant difference to the resulting accuracy of harmonic # analysis. -#Convert a sexagesimal angle into decimal degrees -def s2d(degrees, arcmins = 0, arcsecs = 0, mas = 0, muas = 0): - return ( - degrees - + (arcmins / 60.0) - + (arcsecs / (60.0*60.0)) - + (mas / (60.0*60.0*1e3)) - + (muas / (60.0*60.0*1e6)) - ) - -#Evaluate a polynomial at argument + +def s2d(degrees, arcmins=0, arcsecs=0, mas=0, muas=0): + # Convert a sexagesimal angle into decimal degrees + return ( + degrees + + (arcmins / 60.0) + + (arcsecs / (60.0*60.0)) + + (mas / (60.0*60.0*1e3)) + + (muas / (60.0*60.0*1e6)) + ) + + def polynomial(coefficients, argument): - return sum([c * (argument ** i) for i,c in enumerate(coefficients)]) + # Evaluate a polynomial at argument + return sum([c * (argument ** i) for i, c in enumerate(coefficients)]) + -#Evaluate the first derivative of a polynomial at argument def d_polynomial(coefficients, argument): - return sum([c * i * (argument ** (i-1)) for i,c in enumerate(coefficients)]) + # Evaluate the first derivative of a polynomial at argument + return sum([c * i * (argument ** (i - 1)) + for i, c in enumerate(coefficients)]) + -#Meeus formula 11.1 def T(t): - return (JD(t) - 2451545.0)/36525 + # Meeus formula 11.1 + return (JD(t) - 2451545.0)/36525 + -#Meeus formula 7.1 def JD(t): - Y, M = t.year, t.month - D = ( - t.day - + t.hour / (24.0) - + t.minute / (24.0*60.0) - + t.second / (24.0*60.0*60.0) - + t.microsecond / (24.0 * 60.0 * 60.0 * 1e6) - ) - if M <= 2: - Y = Y - 1 - M = M + 12 - A = np.floor(Y / 100.0) - B = 2 - A + np.floor(A / 4.0) - return np.floor(365.25*(Y+4716)) + np.floor(30.6001*(M+1)) + D + B - 1524.5 - -#Meeus formula 21.3 + # Meeus formula 7.1 + Y, M = t.year, t.month + D = ( + t.day + + t.hour / (24.0) + + t.minute / (24.0*60.0) + + t.second / (24.0*60.0*60.0) + + t.microsecond / (24.0 * 60.0 * 60.0 * 1e6) + ) + if M <= 2: + Y = Y - 1 + M = M + 12 + A = np.floor(Y / 100.0) + B = 2 - A + np.floor(A / 4.0) + return np.floor(365.25*(Y+4716)) + np.floor(30.6001*(M+1)) + D + B - 1524.5 + +# Meeus formula 21.3 terrestrial_obliquity_coefficients = ( - s2d(23,26,21.448), - -s2d(0,0,4680.93), - -s2d(0,0,1.55), - s2d(0,0,1999.25), - -s2d(0,0,51.38), - -s2d(0,0,249.67), - -s2d(0,0,39.05), - s2d(0,0,7.12), - s2d(0,0,27.87), - s2d(0,0,5.79), - s2d(0,0,2.45) + s2d(23, 26, 21.448), + -s2d(0, 0, 4680.93), + -s2d(0, 0, 1.55), + s2d(0, 0, 1999.25), + -s2d(0, 0, 51.38), + -s2d(0, 0, 249.67), + -s2d(0, 0, 39.05), + s2d(0, 0, 7.12), + s2d(0, 0, 27.87), + s2d(0, 0, 5.79), + s2d(0, 0, 2.45) ) -#Adjust these coefficients for parameter T rather than U +# Adjust these coefficients for parameter T rather than U terrestrial_obliquity_coefficients = [ - c * (1e-2) ** i for i,c in enumerate(terrestrial_obliquity_coefficients) + c * (1e-2) ** i for i, c in enumerate(terrestrial_obliquity_coefficients) ] -#Not entirely sure about this interpretation, but this is the difference -#between Meeus formulae 24.2 and 24.3 and seems to work +# Not entirely sure about this interpretation, but this is the difference +# between Meeus formulae 24.2 and 24.3 and seems to work solar_perigee_coefficients = ( - 280.46645 - 357.52910, - 36000.76932 - 35999.05030, - 0.0003032 + 0.0001559, - 0.00000048 + 280.46645 - 357.52910, + 36000.76932 - 35999.05030, + 0.0003032 + 0.0001559, + 0.00000048 ) -#Meeus formula 24.2 +# Meeus formula 24.2 solar_longitude_coefficients = ( - 280.46645, - 36000.76983, - 0.0003032 + 280.46645, + 36000.76983, + 0.0003032 ) -#This value is taken from JPL Horizon and is essentially constant +# This value is taken from JPL Horizon and is essentially constant lunar_inclination_coefficients = ( - 5.145, + 5.145, ) -#Meeus formula 45.1 +# Meeus formula 45.1 lunar_longitude_coefficients = ( - 218.3164591, - 481267.88134236, - -0.0013268, - 1/538841.0 - -1/65194000.0 + 218.3164591, + 481267.88134236, + -0.0013268, + 1 / 538841.0 + -1 / 65194000.0 ) -#Meeus formula 45.7 +# Meeus formula 45.7 lunar_node_coefficients = ( - 125.0445550, - -1934.1361849, - 0.0020762, - 1/467410.0, - -1/60616000.0 + 125.0445550, + -1934.1361849, + 0.0020762, + 1/467410.0, + -1/60616000.0 ) -#Meeus, unnumbered formula directly preceded by 45.7 +# Meeus, unnumbered formula directly preceded by 45.7 lunar_perigee_coefficients = ( - 83.3532430, - 4069.0137111, - -0.0103238, - -1/80053.0, - 1/18999000.0 + 83.3532430, + 4069.0137111, + -0.0103238, + -1/80053.0, + 1/18999000.0 ) -#Now follow some useful auxiliary values, we won't need their speed. -#See notes on Table 6 in Schureman for I, nu, xi, nu', 2nu'' +# Now follow some useful auxiliary values, we won't need their speed. +# See notes on Table 6 in Schureman for I, nu, xi, nu', 2nu'' + + def _I(N, i, omega): - N, i, omega = d2r * N, d2r*i, d2r*omega - cosI = np.cos(i)*np.cos(omega)-np.sin(i)*np.sin(omega)*np.cos(N) - return r2d*np.arccos(cosI) + N, i, omega = d2r * N, d2r*i, d2r*omega + cosI = np.cos(i) * np.cos(omega) - np.sin(i) * np.sin(omega) * np.cos(N) + return r2d * np.arccos(cosI) + def _xi(N, i, omega): - N, i, omega = d2r * N, d2r*i, d2r*omega - e1 = np.cos(0.5*(omega-i))/np.cos(0.5*(omega+i)) * np.tan(0.5*N) - e2 = np.sin(0.5*(omega-i))/np.sin(0.5*(omega+i)) * np.tan(0.5*N) - e1, e2 = np.arctan(e1), np.arctan(e2) - e1, e2 = e1 - 0.5*N, e2 - 0.5*N - return -(e1 + e2)*r2d + N, i, omega = d2r * N, d2r*i, d2r*omega + e1 = (np.cos(0.5 * (omega - i)) + / np.cos(0.5 * (omega + i)) * np.tan(0.5 * N)) + e2 = (np.sin(0.5 * (omega - i)) + / np.sin(0.5 * (omega + i)) * np.tan(0.5 * N)) + e1, e2 = np.arctan(e1), np.arctan(e2) + e1, e2 = e1 - 0.5*N, e2 - 0.5 * N + return -(e1 + e2)*r2d + def _nu(N, i, omega): - N, i, omega = d2r * N, d2r*i, d2r*omega - e1 = np.cos(0.5*(omega-i))/np.cos(0.5*(omega+i)) * np.tan(0.5*N) - e2 = np.sin(0.5*(omega-i))/np.sin(0.5*(omega+i)) * np.tan(0.5*N) - e1, e2 = np.arctan(e1), np.arctan(e2) - e1, e2 = e1 - 0.5*N, e2 - 0.5*N - return (e1 - e2)*r2d - -#Schureman equation 224 -#Can we be more precise than B "the solar coefficient" = 0.1681? + N, i, omega = d2r * N, d2r*i, d2r*omega + e1 = (np.cos(0.5 * (omega - i)) + / np.cos(0.5 * (omega + i)) * np.tan(0.5 * N)) + e2 = (np.sin(0.5 * (omega - i)) + / np.sin(0.5 * (omega + i)) + * np.tan(0.5 * N)) + e1, e2 = np.arctan(e1), np.arctan(e2) + e1, e2 = e1 - 0.5 * N, e2 - 0.5 * N + return (e1 - e2) * r2d + + +# Can we be more precise than B "the solar coefficient" = 0.1681? def _nup(N, i, omega): - I = d2r * _I(N, i, omega) - nu = d2r * _nu(N, i, omega) - return r2d * np.arctan(np.sin(2*I)*np.sin(nu)/(np.sin(2*I)*np.cos(nu)+0.3347)) + # Schureman equation 224 + I = d2r * _I(N, i, omega) + nu = d2r * _nu(N, i, omega) + return (r2d * np.arctan( + np.sin(2 * I) * np.sin(nu) / (np.sin(2 * I) * np.cos(nu) + 0.3347)) + ) + -#Schureman equation 232 def _nupp(N, i, omega): - I = d2r * _I(N, i, omega) - nu = d2r * _nu(N, i, omega) - tan2nupp = (np.sin(I)**2*np.sin(2*nu))/(np.sin(I)**2*np.cos(2*nu)+0.0727) - return r2d * 0.5 * np.arctan(tan2nupp) + # Schureman equation 232 + I = d2r * _I(N, i, omega) + nu = d2r * _nu(N, i, omega) + tan2nupp = (np.sin(I)**2*np.sin(2*nu))/(np.sin(I)**2*np.cos(2*nu)+0.0727) + return r2d * 0.5 * np.arctan(tan2nupp) AstronomicalParameter = namedtuple('AstronomicalParameter', ['value', 'speed']) + def astro(t): - a = {} - #We can use polynomial fits from Meeus to obtain good approximations to - #some astronomical values (and therefore speeds). - polynomials = { - 's': lunar_longitude_coefficients, - 'h': solar_longitude_coefficients, - 'p': lunar_perigee_coefficients, - 'N': lunar_node_coefficients, - 'pp': solar_perigee_coefficients, - '90': (90.0,), - 'omega': terrestrial_obliquity_coefficients, - 'i': lunar_inclination_coefficients - } - #Polynomials are in T, that is Julian Centuries; we want our speeds to be - #in the more convenient unit of degrees per hour. - dT_dHour = 1 / (24 * 365.25 * 100) - for name, coefficients in polynomials.items(): - a[name] = AstronomicalParameter( - np.mod(polynomial(coefficients, T(t)), 360.0), - d_polynomial(coefficients, T(t)) * dT_dHour - ) - - #Some other parameters defined by Schureman which are dependent on the - #parameters N, i, omega for use in node factor calculations. We don't need - #their speeds. - args = list(each.value for each in [a['N'], a['i'], a['omega']]) - for name, function in { - 'I': _I, - 'xi': _xi, - 'nu': _nu, - 'nup': _nup, - 'nupp': _nupp - }.items(): - a[name] = AstronomicalParameter(np.mod(function(*args), 360.0), None) - - #We don't work directly with the T (hours) parameter, instead our spanning - #set for equilibrium arguments #is given by T+h-s, s, h, p, N, pp, 90. - #This is in line with convention. - hour = AstronomicalParameter((JD(t) - np.floor(JD(t))) * 360.0, 15.0) - a['T+h-s'] = AstronomicalParameter( - hour.value + a['h'].value - a['s'].value, - hour.speed + a['h'].speed - a['s'].speed - ) - #It is convenient to calculate Schureman's P here since several node - #factors need it, although it could be argued that these - #(along with I, xi, nu etc) belong somewhere else. - a['P'] = AstronomicalParameter( - np.mod(a['p'].value -a['xi'].value,360.0), - None - ) - return a + a = {} + # We can use polynomial fits from Meeus to obtain good approximations to + # some astronomical values (and therefore speeds). + polynomials = { + 's': lunar_longitude_coefficients, + 'h': solar_longitude_coefficients, + 'p': lunar_perigee_coefficients, + 'N': lunar_node_coefficients, + 'pp': solar_perigee_coefficients, + '90': (90.0,), + 'omega': terrestrial_obliquity_coefficients, + 'i': lunar_inclination_coefficients + } + # Polynomials are in T, that is Julian Centuries; we want our speeds to be + # in the more convenient unit of degrees per hour. + dT_dHour = 1 / (24 * 365.25 * 100) + for name, coefficients in polynomials.items(): + a[name] = AstronomicalParameter( + np.mod(polynomial(coefficients, T(t)), 360.0), + d_polynomial(coefficients, T(t)) * dT_dHour + ) + + # Some other parameters defined by Schureman which are dependent on the + # parameters N, i, omega for use in node factor calculations. We don't need + # their speeds. + args = list(each.value for each in [a['N'], a['i'], a['omega']]) + for name, function in { + 'I': _I, + 'xi': _xi, + 'nu': _nu, + 'nup': _nup, + 'nupp': _nupp + }.items(): + a[name] = AstronomicalParameter(np.mod(function(*args), 360.0), None) + + # We don't work directly with the T (hours) parameter, instead our spanning + # set for equilibrium arguments # is given by T+h-s, s, h, p, N, pp, 90. + # This is in line with convention. + hour = AstronomicalParameter((JD(t) - np.floor(JD(t))) * 360.0, 15.0) + a['T+h-s'] = AstronomicalParameter( + hour.value + a['h'].value - a['s'].value, + hour.speed + a['h'].speed - a['s'].speed + ) + # It is convenient to calculate Schureman's P here since several node + # factors need it, although it could be argued that these + # (along with I, xi, nu etc) belong somewhere else. + a['P'] = AstronomicalParameter( + np.mod(a['p'].value - a['xi'].value, 360.0), + None + ) + return a diff --git a/pytides/constituent.py b/pytides/constituent.py index 3ff8197..6c6230a 100644 --- a/pytides/constituent.py +++ b/pytides/constituent.py @@ -4,157 +4,166 @@ import numpy as np import nodal_corrections as nc -class BaseConstituent(object): - xdo_int = { - 'A': 1, 'B': 2, 'C': 3, 'D': 4, 'E': 5, 'F': 6, 'G': 7, 'H': 8, 'I': 9, - 'J': 10, 'K': 11, 'L': 12, 'M': 13, 'N': 14, 'O': 15, 'P': 16, 'Q': 17, - 'R': -8, 'S': -7, 'T': -6, 'U': -5, 'V': -4, 'W': -3, 'X': -2, 'Y': -1, - 'Z': 0 - } - - int_xdo = {v:k for k, v in xdo_int.items()} - - def __init__(self, name, xdo='', coefficients=[], u=nc.u_zero, f=nc.f_unity): - if xdo == '': - self.coefficients = np.array(coefficients) - else: - self.coefficients = np.array(self.xdo_to_coefficients(xdo)) - self.name = name - self.u = u - self.f = f - - def xdo_to_coefficients(self, xdo): - return [self.xdo_int[l.upper()] for l in xdo if l in string.ascii_letters] - - def coefficients_to_xdo(self, coefficients): - return ''.join([self.int_xdo[c] for c in cooefficients]) - - def V(self, astro): - return np.dot(self.coefficients, self.astro_values(astro)) - - def xdo(self): - return self.coefficients_to_xdo(self.coefficients) - - def speed(self, a): - return np.dot(self.coefficients, self.astro_speeds(a)) - - def astro_xdo(self, a): - return [a['T+h-s'], a['s'], a['h'], a['p'], a['N'], a['pp'], a['90']] - - def astro_speeds(self, a): - return np.array([each.speed for each in self.astro_xdo(a)]) - - def astro_values(self, a): - return np.array([each.value for each in self.astro_xdo(a)]) - - #Consider two out of phase constituents which travel at the same speed to - #be identical - def __eq__(self, c): - return np.all(self.coefficients[:-1] == c.coefficients[:-1]) - def __hash__(self): - return hash(tuple(self.coefficients[:-1])) - -class CompoundConstituent(BaseConstituent): - - def __init__(self, members = [], **kwargs): - self.members = members - - if 'u' not in kwargs: - kwargs['u'] = self.u - if 'f' not in kwargs: - kwargs['f'] = self.f - - super(CompoundConstituent,self).__init__(**kwargs) - - self.coefficients = reduce(op.add,[c.coefficients * n for (c,n) in members]) - - def speed(self, a): - return reduce(op.add, [n * c.speed(a) for (c,n) in self.members]) - - def V(self, a): - return reduce(op.add, [n * c.V(a) for (c,n) in self.members]) - - def u(self, a): - return reduce(op.add, [n * c.u(a) for (c,n) in self.members]) +class BaseConstituent(object): + xdo_int = { + 'A': 1, 'B': 2, 'C': 3, 'D': 4, 'E': 5, 'F': 6, 'G': 7, 'H': 8, 'I': 9, + 'J': 10, 'K': 11, 'L': 12, 'M': 13, 'N': 14, 'O': 15, 'P': 16, 'Q': 17, + 'R': -8, 'S': -7, 'T': -6, 'U': -5, 'V': -4, 'W': -3, 'X': -2, 'Y': -1, + 'Z': 0 + } - def f(self, a): - return reduce(op.mul, [c.f(a) ** abs(n) for (c,n) in self.members]) + int_xdo = {v: k for k, v in xdo_int.items()} -###### Base Constituents -#Long Term -_Z0 = BaseConstituent(name = 'Z0', xdo = 'Z ZZZ ZZZ', u = nc.u_zero, f = nc.f_unity) -_Sa = BaseConstituent(name = 'Sa', xdo = 'Z ZAZ ZZZ', u = nc.u_zero, f = nc.f_unity) -_Ssa = BaseConstituent(name = 'Ssa', xdo = 'Z ZBZ ZZZ', u = nc.u_zero, f = nc.f_unity) -_Mm = BaseConstituent(name = 'Mm', xdo = 'Z AZY ZZZ', u = nc.u_zero, f = nc.f_Mm) -_Mf = BaseConstituent(name = 'Mf', xdo = 'Z BZZ ZZZ', u = nc.u_Mf, f = nc.f_Mf) + def __init__(self, name, xdo='', coefficients=[], + u=nc.u_zero, f=nc.f_unity): + if xdo == '': + self.coefficients = np.array(coefficients) + else: + self.coefficients = np.array(self.xdo_to_coefficients(xdo)) + self.name = name + self.u = u + self.f = f -#Diurnals -_Q1 = BaseConstituent(name = 'Q1', xdo = 'A XZA ZZA', u = nc.u_O1, f = nc.f_O1) -_O1 = BaseConstituent(name = 'O1', xdo = 'A YZZ ZZA', u = nc.u_O1, f = nc.f_O1) -_K1 = BaseConstituent(name = 'K1', xdo = 'A AZZ ZZY', u = nc.u_K1, f = nc.f_K1) -_J1 = BaseConstituent(name = 'J1', xdo = 'A BZY ZZY', u = nc.u_J1, f = nc.f_J1) + def xdo_to_coefficients(self, xdo): + return [self.xdo_int[l.upper()] + for l in xdo + if l in string.ascii_letters] -#M1 is a tricky business for reasons of convention, rather than theory. The -#reasons for this are best summarised by Schureman paragraphs 126, 127 and in -#the comments found in congen_input.txt of xtides, so I won't go over all this -#again here. + def coefficients_to_xdo(self, coefficients): + return ''.join([self.int_xdo[c] for c in cooefficients]) -_M1 = BaseConstituent(name = 'M1', xdo = 'A ZZZ ZZA', u = nc.u_M1, f = nc.f_M1) -_P1 = BaseConstituent(name = 'P1', xdo = 'A AXZ ZZA', u = nc.u_zero, f = nc.f_unity) -_S1 = BaseConstituent(name = 'S1', xdo = 'A AYZ ZZZ', u = nc.u_zero, f = nc.f_unity) -_OO1 = BaseConstituent(name = 'OO1', xdo = 'A CZZ ZZY', u = nc.u_OO1, f = nc.f_OO1) + def V(self, astro): + return np.dot(self.coefficients, self.astro_values(astro)) -#Semi-Diurnals -_2N2 = BaseConstituent(name = '2N2', xdo = 'B XZB ZZZ', u = nc.u_M2, f = nc.f_M2) -_N2 = BaseConstituent(name = 'N2', xdo = 'B YZA ZZZ', u = nc.u_M2, f = nc.f_M2) -_nu2 = BaseConstituent(name = 'nu2', xdo = 'B YBY ZZZ', u = nc.u_M2, f = nc.f_M2) -_M2 = BaseConstituent(name = 'M2', xdo = 'B ZZZ ZZZ', u = nc.u_M2, f = nc.f_M2) -_lambda2 = BaseConstituent(name = 'lambda2', xdo = 'B AXA ZZB', u = nc.u_M2, f = nc.f_M2) -_L2 = BaseConstituent(name = 'L2', xdo = 'B AZY ZZB', u = nc.u_L2, f = nc.f_L2) -_T2 = BaseConstituent(name = 'T2', xdo = 'B BWZ ZAZ', u = nc.u_zero, f = nc.f_unity) -_S2 = BaseConstituent(name = 'S2', xdo = 'B BXZ ZZZ', u = nc.u_zero, f = nc.f_unity) -_R2 = BaseConstituent(name = 'R2', xdo = 'B BYZ ZYB', u = nc.u_zero, f = nc.f_unity) -_K2 = BaseConstituent(name = 'K2', xdo = 'B BZZ ZZZ', u = nc.u_K2, f = nc.f_K2) + def xdo(self): + return self.coefficients_to_xdo(self.coefficients) -#Third-Diurnals -_M3 = BaseConstituent(name = 'M3', xdo = 'C ZZZ ZZZ', u = lambda a: nc.u_Modd(a,3), f = lambda a: nc.f_Modd(a,3)) + def speed(self, a): + return np.dot(self.coefficients, self.astro_speeds(a)) -###### Compound Constituents -#Long Term -_MSF = CompoundConstituent(name = 'MSF', members = [(_S2, 1), (_M2, -1)]) + def astro_xdo(self, a): + return [a['T+h-s'], a['s'], a['h'], a['p'], a['N'], a['pp'], a['90']] -#Diurnal -_2Q1 = CompoundConstituent(name = '2Q1', members = [(_N2, 1), (_J1, -1)]) -_rho1 = CompoundConstituent(name = 'rho1', members = [(_nu2, 1), (_K1, -1)]) + def astro_speeds(self, a): + return np.array([each.speed for each in self.astro_xdo(a)]) -#Semi-Diurnal + def astro_values(self, a): + return np.array([each.value for each in self.astro_xdo(a)]) -_mu2 = CompoundConstituent(name = 'mu2', members = [(_M2, 2), (_S2, -1)]) #2MS2 -_2SM2 = CompoundConstituent(name = '2SM2', members = [(_S2, 2), (_M2, -1)]) + # Consider two out of phase constituents which travel at the same speed to + # be identical + def __eq__(self, c): + return np.all(self.coefficients[:-1] == c.coefficients[:-1]) -#Third-Diurnal -_2MK3 = CompoundConstituent(name = '2MK3', members = [(_M2, 1), (_O1, 1)]) -_MK3 = CompoundConstituent(name = 'MK3', members = [(_M2, 1), (_K1, 1)]) + def __hash__(self): + return hash(tuple(self.coefficients[:-1])) -#Quarter-Diurnal -_MN4 = CompoundConstituent(name = 'MN4', members = [(_M2, 1), (_N2, 1)]) -_M4 = CompoundConstituent(name = 'M4', members = [(_M2, 2)]) -_MS4 = CompoundConstituent(name = 'MS4', members = [(_M2, 1), (_S2, 1)]) -_S4 = CompoundConstituent(name = 'S4', members = [(_S2, 2)]) -#Sixth-Diurnal -_M6 = CompoundConstituent(name = 'M6', members = [(_M2, 3)]) -_S6 = CompoundConstituent(name = 'S6', members = [(_S2, 3)]) +class CompoundConstituent(BaseConstituent): -#Eighth-Diurnals -_M8 = CompoundConstituent(name = 'M8', members = [(_M2, 4)]) + def __init__(self, members=[], **kwargs): + self.members = members + + if 'u' not in kwargs: + kwargs['u'] = self.u + if 'f' not in kwargs: + kwargs['f'] = self.f + + super(CompoundConstituent, self).__init__(**kwargs) + + self.coefficients = reduce( + op.add, + [c.coefficients * n for (c, n) in members]) + + def speed(self, a): + return reduce(op.add, [n * c.speed(a) for (c, n) in self.members]) + + def V(self, a): + return reduce(op.add, [n * c.V(a) for (c, n) in self.members]) + + def u(self, a): + return reduce(op.add, [n * c.u(a) for (c, n) in self.members]) + + def f(self, a): + return reduce(op.mul, [c.f(a) ** abs(n) for (c, n) in self.members]) + +# ## ## # Base Constituents +# Long Term +_Z0 = BaseConstituent(name='Z0', xdo='Z ZZZ ZZZ', u=nc.u_zero, f=nc.f_unity) +_Sa = BaseConstituent(name='Sa', xdo='Z ZAZ ZZZ', u=nc.u_zero, f=nc.f_unity) +_Ssa = BaseConstituent(name='Ssa', xdo='Z ZBZ ZZZ', u=nc.u_zero, f=nc.f_unity) +_Mm = BaseConstituent(name='Mm', xdo='Z AZY ZZZ', u=nc.u_zero, f=nc.f_Mm) +_Mf = BaseConstituent(name='Mf', xdo='Z BZZ ZZZ', u=nc.u_Mf, f=nc.f_Mf) + +# Diurnals +_Q1 = BaseConstituent(name='Q1', xdo='A XZA ZZA', u=nc.u_O1, f=nc.f_O1) +_O1 = BaseConstituent(name='O1', xdo='A YZZ ZZA', u=nc.u_O1, f=nc.f_O1) +_K1 = BaseConstituent(name='K1', xdo='A AZZ ZZY', u=nc.u_K1, f=nc.f_K1) +_J1 = BaseConstituent(name='J1', xdo='A BZY ZZY', u=nc.u_J1, f=nc.f_J1) + +# M1 is a tricky business for reasons of convention, rather than theory. The +# reasons for this are best summarised by Schureman paragraphs 126, 127 and in +# the comments found in congen_input.txt of xtides, so I won't go over all this +# again here. + +_M1 = BaseConstituent(name='M1', xdo='A ZZZ ZZA', u=nc.u_M1, f=nc.f_M1) +_P1 = BaseConstituent(name='P1', xdo='A AXZ ZZA', u=nc.u_zero, f=nc.f_unity) +_S1 = BaseConstituent(name='S1', xdo='A AYZ ZZZ', u=nc.u_zero, f=nc.f_unity) +_OO1 = BaseConstituent(name='OO1', xdo='A CZZ ZZY', u=nc.u_OO1, f=nc.f_OO1) + +# Semi-Diurnals +_2N2 = BaseConstituent(name='2N2', xdo='B XZB ZZZ', u=nc.u_M2, f=nc.f_M2) +_N2 = BaseConstituent(name='N2', xdo='B YZA ZZZ', u=nc.u_M2, f=nc.f_M2) +_nu2 = BaseConstituent(name='nu2', xdo='B YBY ZZZ', u=nc.u_M2, f=nc.f_M2) +_M2 = BaseConstituent(name='M2', xdo='B ZZZ ZZZ', u=nc.u_M2, f=nc.f_M2) +_lambda2 = BaseConstituent(name='lambda2', xdo='B AXA ZZB', + u=nc.u_M2, f=nc.f_M2) +_L2 = BaseConstituent(name='L2', xdo='B AZY ZZB', u=nc.u_L2, f=nc.f_L2) +_T2 = BaseConstituent(name='T2', xdo='B BWZ ZAZ', u=nc.u_zero, f=nc.f_unity) +_S2 = BaseConstituent(name='S2', xdo='B BXZ ZZZ', u=nc.u_zero, f=nc.f_unity) +_R2 = BaseConstituent(name='R2', xdo='B BYZ ZYB', u=nc.u_zero, f=nc.f_unity) +_K2 = BaseConstituent(name='K2', xdo='B BZZ ZZZ', u=nc.u_K2, f=nc.f_K2) + +# Third-Diurnals +_M3 = BaseConstituent(name='M3', xdo='C ZZZ ZZZ', + u=lambda a: nc.u_Modd(a, 3), + f=lambda a: nc.f_Modd(a, 3)) + +# ## ## # Compound Constituents +# Long Term +_MSF = CompoundConstituent(name='MSF', members=[(_S2, 1), (_M2, -1)]) + +# Diurnal +_2Q1 = CompoundConstituent(name='2Q1', members=[(_N2, 1), (_J1, -1)]) +_rho1 = CompoundConstituent(name='rho1', members=[(_nu2, 1), (_K1, -1)]) + +# Semi-Diurnal + +_mu2 = CompoundConstituent(name='mu2', members=[(_M2, 2), (_S2, -1)]) # 2MS2 +_2SM2 = CompoundConstituent(name='2SM2', members=[(_S2, 2), (_M2, -1)]) + +# Third-Diurnal +_2MK3 = CompoundConstituent(name='2MK3', members=[(_M2, 1), (_O1, 1)]) +_MK3 = CompoundConstituent(name='MK3', members=[(_M2, 1), (_K1, 1)]) + +# Quarter-Diurnal +_MN4 = CompoundConstituent(name='MN4', members=[(_M2, 1), (_N2, 1)]) +_M4 = CompoundConstituent(name='M4', members=[(_M2, 2)]) +_MS4 = CompoundConstituent(name='MS4', members=[(_M2, 1), (_S2, 1)]) +_S4 = CompoundConstituent(name='S4', members=[(_S2, 2)]) + +# Sixth-Diurnal +_M6 = CompoundConstituent(name='M6', members=[(_M2, 3)]) +_S6 = CompoundConstituent(name='S6', members=[(_S2, 3)]) + +# Eighth-Diurnals +_M8 = CompoundConstituent(name='M8', members=[(_M2, 4)]) noaa = [ - _M2, _S2, _N2, _K1, _M4, _O1, _M6, _MK3, _S4, _MN4, _nu2, _S6, _mu2, - _2N2, _OO1, _lambda2, _S1, _M1, _J1, _Mm, _Ssa, _Sa, _MSF, _Mf, - _rho1, _Q1, _T2, _R2, _2Q1, _P1, _2SM2, _M3, _L2, _2MK3, _K2, - _M8, _MS4 + _M2, _S2, _N2, _K1, _M4, _O1, _M6, _MK3, _S4, _MN4, _nu2, _S6, _mu2, + _2N2, _OO1, _lambda2, _S1, _M1, _J1, _Mm, _Ssa, _Sa, _MSF, _Mf, + _rho1, _Q1, _T2, _R2, _2Q1, _P1, _2SM2, _M3, _L2, _2MK3, _K2, + _M8, _MS4 ] - diff --git a/pytides/nodal_corrections.py b/pytides/nodal_corrections.py index 170f0d6..4114d81 100644 --- a/pytides/nodal_corrections.py +++ b/pytides/nodal_corrections.py @@ -2,140 +2,171 @@ import numpy as np d2r, r2d = np.pi/180.0, 180.0/np.pi -#The following functions take a dictionary of astronomical values (in degrees) -#and return dimensionless scale factors for constituent amplitudes. +# The following functions take a dictionary of astronomical values (in degrees) +# and return dimensionless scale factors for constituent amplitudes. + def f_unity(a): - return 1.0 + return 1.0 + -#Schureman equations 73, 65 def f_Mm(a): - omega = d2r*a['omega'].value - i = d2r*a['i'].value - I = d2r*a['I'].value - mean = (2/3.0 - np.sin(omega)**2)*(1 - 3/2.0 * np.sin(i)**2) - return (2/3.0 - np.sin(I)**2) / mean + # Schureman equations 73, 65 + omega = d2r*a['omega'].value + i = d2r*a['i'].value + I = d2r*a['I'].value + mean = (2/3.0 - np.sin(omega)**2)*(1 - 3/2.0 * np.sin(i)**2) + return (2/3.0 - np.sin(I)**2) / mean + -#Schureman equations 74, 66 def f_Mf(a): - omega = d2r*a['omega'].value - i = d2r*a['i'].value - I = d2r*a['I'].value - mean = np.sin(omega)**2 * np.cos(0.5*i)**4 - return np.sin(I)**2 / mean + # Schureman equations 74, 66 + omega = d2r*a['omega'].value + i = d2r*a['i'].value + I = d2r*a['I'].value + mean = np.sin(omega)**2 * np.cos(0.5*i)**4 + return np.sin(I)**2 / mean + -#Schureman equations 75, 67 def f_O1(a): - omega = d2r*a['omega'].value - i = d2r*a['i'].value - I = d2r*a['I'].value - mean = np.sin(omega) * np.cos(0.5*omega)**2 * np.cos(0.5*i)**4 - return (np.sin(I) * np.cos(0.5*I)**2) / mean + # Schureman equations 75, 67 + omega = d2r*a['omega'].value + i = d2r*a['i'].value + I = d2r*a['I'].value + mean = np.sin(omega) * np.cos(0.5*omega)**2 * np.cos(0.5*i)**4 + return (np.sin(I) * np.cos(0.5*I)**2) / mean + -#Schureman equations 76, 68 def f_J1(a): - omega = d2r*a['omega'].value - i = d2r*a['i'].value - I = d2r*a['I'].value - mean = np.sin(2*omega) * (1-3/2.0 * np.sin(i)**2) - return np.sin(2*I) / mean + # Schureman equations 76, 68 + omega = d2r*a['omega'].value + i = d2r*a['i'].value + I = d2r*a['I'].value + mean = np.sin(2*omega) * (1-3/2.0 * np.sin(i)**2) + return np.sin(2*I) / mean + -#Schureman equations 77, 69 def f_OO1(a): - omega = d2r*a['omega'].value - i = d2r*a['i'].value - I = d2r*a['I'].value - mean = np.sin(omega) * np.sin(0.5*omega)**2 * np.cos(0.5*i)**4 - return np.sin(I) * np.sin(0.5*I)**2 / mean + # Schureman equations 77, 69 + omega = d2r*a['omega'].value + i = d2r*a['i'].value + I = d2r*a['I'].value + mean = np.sin(omega) * np.sin(0.5*omega)**2 * np.cos(0.5*i)**4 + return np.sin(I) * np.sin(0.5*I)**2 / mean + -#Schureman equations 78, 70 def f_M2(a): - omega = d2r*a['omega'].value - i = d2r*a['i'].value - I = d2r*a['I'].value - mean = np.cos(0.5*omega)**4 * np.cos(0.5*i)**4 - return np.cos(0.5*I)**4 / mean - -#Schureman equations 227, 226, 68 -#Should probably eventually include the derivations of the magic numbers (0.5023 etc). + # Schureman equations 78, 70 + omega = d2r*a['omega'].value + i = d2r*a['i'].value + I = d2r*a['I'].value + mean = np.cos(0.5*omega)**4 * np.cos(0.5*i)**4 + return np.cos(0.5*I)**4 / mean + + def f_K1(a): - omega = d2r*a['omega'].value - i = d2r*a['i'].value - I = d2r*a['I'].value - nu = d2r*a['nu'].value - sin2Icosnu_mean = np.sin(2*omega) * (1-3/2.0 * np.sin(i)**2) - mean = 0.5023*sin2Icosnu_mean + 0.1681 - return (0.2523*np.sin(2*I)**2 + 0.1689*np.sin(2*I)*np.cos(nu)+0.0283)**(0.5) / mean - -#Schureman equations 215, 213, 204 -#It can be (and has been) confirmed that the exponent for R_a reads 1/2 via Schureman Table 7 + # Schureman equations 227, 226, 68 + # Should probably eventually include the derivations of the magic numbers + # (0.5023 etc). + omega = d2r*a['omega'].value + i = d2r*a['i'].value + I = d2r*a['I'].value + nu = d2r*a['nu'].value + sin2Icosnu_mean = np.sin(2*omega) * (1-3/2.0 * np.sin(i)**2) + mean = 0.5023*sin2Icosnu_mean + 0.1681 + return (0.2523 * np.sin(2 * I) ** 2 + 0.1689 * np.sin(2 * I) * np.cos(nu) + + 0.0283) ** (0.5) / mean + + def f_L2(a): - P = d2r*a['P'].value - I = d2r*a['I'].value - R_a_inv = (1 - 12*np.tan(0.5*I)**2 * np.cos(2*P)+36*np.tan(0.5*I)**4)**(0.5) - return f_M2(a) * R_a_inv + # Schureman equations 215, 213, 204 + # It can be (and has been) confirmed that the exponent for R_a reads 1/2 + # via Schureman Table 7 + P = d2r*a['P'].value + I = d2r*a['I'].value + R_a_inv = (1 - 12 * np.tan(0.5 * I) ** 2 * np.cos(2 * P) + 36 + * np.tan(0.5 * I) ** 4) ** (0.5) + return f_M2(a) * R_a_inv + -#Schureman equations 235, 234, 71 -#Again, magic numbers def f_K2(a): - omega = d2r*a['omega'].value - i = d2r*a['i'].value - I = d2r*a['I'].value - nu = d2r*a['nu'].value - sinsqIcos2nu_mean = np.sin(omega)**2 * (1-3/2.0 * np.sin(i)**2) - mean = 0.5023*sinsqIcos2nu_mean + 0.0365 - return (0.2533*np.sin(I)**4 + 0.0367*np.sin(I)**2 *np.cos(2*nu)+0.0013)**(0.5) / mean - -#Schureman equations 206, 207, 195 + # Schureman equations 235, 234, 71 + # Again, magic numbers + omega = d2r*a['omega'].value + i = d2r*a['i'].value + I = d2r*a['I'].value + nu = d2r*a['nu'].value + sinsqIcos2nu_mean = np.sin(omega)**2 * (1-3/2.0 * np.sin(i)**2) + mean = 0.5023 * sinsqIcos2nu_mean + 0.0365 + return (0.2533 * (np.sin(I) ** 4) + 0.0367 * np.sin(I) ** 2 + * np.cos(2 * nu) + 0.0013) ** (0.5) / mean + + def f_M1(a): - P = d2r*a['P'].value - I = d2r*a['I'].value - Q_a_inv = (0.25 + 1.5*np.cos(I)*np.cos(2*P)*np.cos(0.5*I)**(-0.5) + 2.25*np.cos(I)**2 * np.cos(0.5*I)**(-4))**(0.5) - return f_O1(a) * Q_a_inv + # Schureman equations 206, 207, 195 + P = d2r*a['P'].value + I = d2r*a['I'].value + Q_a_inv = (0.25 + 1.5 * np.cos(I) * np.cos(2 * P) + * np.cos(0.5 * I) ** (-0.5) + 2.25 * np.cos(I) ** 2 + * np.cos(0.5 * I) ** (-4)) ** (0.5) + return f_O1(a) * Q_a_inv + -#See e.g. Schureman equation 149 def f_Modd(a, n): - return f_M2(a) ** (n / 2.0) + # See e.g. Schureman equation 149 + return f_M2(a) ** (n / 2.0) + +# Node factors u, see Table 2 of Schureman. -#Node factors u, see Table 2 of Schureman. def u_zero(a): - return 0.0 + return 0.0 + def u_Mf(a): - return -2.0 * a['xi'].value + return -2.0 * a['xi'].value + def u_O1(a): - return 2.0 * a['xi'].value - a['nu'].value + return 2.0 * a['xi'].value - a['nu'].value + def u_J1(a): - return -a['nu'].value + return -a['nu'].value + def u_OO1(a): - return -2.0 * a['xi'].value - a['nu'].value + return -2.0 * a['xi'].value - a['nu'].value + def u_M2(a): - return 2.0 * a['xi'].value - 2.0 * a['nu'].value + return 2.0 * a['xi'].value - 2.0 * a['nu'].value + def u_K1(a): - return -a['nup'].value + return -a['nup'].value + -#Schureman 214 def u_L2(a): - I = d2r*a['I'].value - P = d2r*a['P'].value - R = r2d*np.arctan(np.sin(2*P)/(1/6.0 * np.tan(0.5*I) **(-2) -np.cos(2*P))) - return 2.0 * a['xi'].value - 2.0 * a['nu'].value - R + # Schureman 214 + I = d2r * a['I'].value + P = d2r * a['P'].value + R = r2d * np.arctan( + np.sin(2 * P) / (1 / 6.0 * np.tan(0.5 * I) ** (-2) - np.cos(2 * P))) + return 2.0 * a['xi'].value - 2.0 * a['nu'].value - R + def u_K2(a): - return -2.0 * a['nupp'].value + return -2.0 * a['nupp'].value + -#Schureman 202 def u_M1(a): - I = d2r*a['I'].value - P = d2r*a['P'].value - Q = r2d*np.arctan((5*np.cos(I)-1)/(7*np.cos(I)+1)*np.tan(P)) - return a['xi'].value - a['nu'].value + Q + # Schureman 202 + I = d2r*a['I'].value + P = d2r*a['P'].value + Q = r2d*np.arctan((5*np.cos(I)-1)/(7*np.cos(I)+1)*np.tan(P)) + return a['xi'].value - a['nu'].value + Q + def u_Modd(a, n): - return n/2.0 * u_M2(a) + return n/2.0 * u_M2(a) diff --git a/pytides/tide.py b/pytides/tide.py index bc7089c..34fb8a0 100644 --- a/pytides/tide.py +++ b/pytides/tide.py @@ -2,10 +2,10 @@ from collections import OrderedDict, Iterable from itertools import takewhile, count try: - from itertools import izip, ifilter -except ImportError: #Python3 - izip = zip - ifilter = filter + from itertools import izip, ifilter +except ImportError: # Python3 + izip = zip + ifilter = filter from datetime import datetime, timedelta import numpy as np from scipy.optimize import leastsq, fsolve @@ -14,393 +14,435 @@ d2r, r2d = np.pi/180.0, 180.0/np.pi + class Tide(object): - dtype = np.dtype([ - ('constituent', object), - ('amplitude', float), - ('phase', float)]) - - def __init__( - self, - constituents = None, - amplitudes = None, - phases = None, - model = None, - radians = False - ): - """ - Initialise a tidal model. Provide constituents, amplitudes and phases OR a model. - Arguments: - constituents -- list of constituents used in the model. - amplitudes -- list of amplitudes corresponding to constituents - phases -- list of phases corresponding to constituents - model -- an ndarray of type Tide.dtype representing the constituents, amplitudes and phases. - radians -- boolean representing whether phases are in radians (default False) - """ - if None not in [constituents, amplitudes, phases]: - if len(constituents) == len(amplitudes) == len(phases): - model = np.zeros(len(phases), dtype=Tide.dtype) - model['constituent'] = np.array(constituents) - model['amplitude'] = np.array(amplitudes) - model['phase'] = np.array(phases) - else: - raise ValueError("Constituents, amplitudes and phases should all be arrays of equal length.") - elif model is not None: - if not model.dtype == Tide.dtype: - raise ValueError("Model must be a numpy array with dtype == Tide.dtype") - else: - raise ValueError("Must be initialised with constituents, amplitudes and phases; or a model.") - if radians: - model['phase'] = r2d*model['phase'] - self.model = model[:] - self.normalize() - - def prepare(self, *args, **kwargs): - return Tide._prepare(self.model['constituent'], *args, **kwargs) - - @staticmethod - def _prepare(constituents, t0, t = None, radians = True): - """ - Return constituent speed and equilibrium argument at a given time, and constituent node factors at given times. - Arguments: - constituents -- list of constituents to prepare - t0 -- time at which to evaluate speed and equilibrium argument for each constituent - t -- list of times at which to evaluate node factors for each constituent (default: t0) - radians -- whether to return the angular arguments in radians or degrees (default: True) - """ - #The equilibrium argument is constant and taken at the beginning of the - #time series (t0). The speed of the equilibrium argument changes very - #slowly, so again we take it to be constant over any length of data. The - #node factors change more rapidly. - if isinstance(t0, Iterable): - t0 = t0[0] - if t is None: - t = [t0] - if not isinstance(t, Iterable): - t = [t] - a0 = astro(t0) - a = [astro(t_i) for t_i in t] - - #For convenience give u, V0 (but not speed!) in [0, 360) - V0 = np.array([c.V(a0) for c in constituents])[:, np.newaxis] - speed = np.array([c.speed(a0) for c in constituents])[:, np.newaxis] - u = [np.mod(np.array([c.u(a_i) for c in constituents])[:, np.newaxis], 360.0) - for a_i in a] - f = [np.mod(np.array([c.f(a_i) for c in constituents])[:, np.newaxis], 360.0) - for a_i in a] - - if radians: - speed = d2r*speed - V0 = d2r*V0 - u = [d2r*each for each in u] - return speed, u, f, V0 - - def at(self, t): - """ - Return the modelled tidal height at given times. - Arguments: - t -- array of times at which to evaluate the tidal height - """ - t0 = t[0] - hours = self._hours(t0, t) - partition = 240.0 - t = self._partition(hours, partition) - times = self._times(t0, [(i + 0.5)*partition for i in range(len(t))]) - speed, u, f, V0 = self.prepare(t0, times, radians = True) - H = self.model['amplitude'][:, np.newaxis] - p = d2r*self.model['phase'][:, np.newaxis] - - return np.concatenate([ - Tide._tidal_series(t_i, H, p, speed, u_i, f_i, V0) - for t_i, u_i, f_i in izip(t, u, f) - ]) - - def highs(self, *args): - """ - Generator yielding only the high tides. - Arguments: - see Tide.extrema() - """ - for t in ifilter(lambda e: e[2] == 'H', self.extrema(*args)): - yield t - - def lows(self, *args): - """ - Generator yielding only the low tides. - Arguments: - see Tide.extrema() - """ - for t in ifilter(lambda e: e[2] == 'L', self.extrema(*args)): - yield t - - def form_number(self): - """ - Returns the model's form number, a helpful heuristic for classifying tides. - """ - k1, o1, m2, s2 = ( - np.extract(self.model['constituent'] == c, self.model['amplitude']) - for c in [constituent._K1, constituent._O1, constituent._M2, constituent._S2] - ) - return (k1+o1)/(m2+s2) - - def classify(self): - """ - Classify the tide according to its form number - """ - form = self.form_number() - if 0 <= form <= 0.25: - return 'semidiurnal' - elif 0.25 < form <= 1.5: - return 'mixed (semidiurnal)' - elif 1.5 < form <= 3.0: - return 'mixed (diurnal)' - else: - return 'diurnal' - - def extrema(self, t0, t1 = None, partition = 2400.0): - """ - A generator for high and low tides. - Arguments: - t0 -- time after which extrema are sought - t1 -- optional time before which extrema are sought (if not given, the generator is infinite) - partition -- number of hours for which we consider the node factors to be constant (default: 2400.0) - """ - if t1: - #yield from in python 3.4 - for e in takewhile(lambda t: t[0] < t1, self.extrema(t0)): - yield e - else: - #We assume that extrema are separated by at least delta hours - delta = np.amin([ - 90.0 / c.speed(astro(t0)) for c in self.model['constituent'] - if not c.speed(astro(t0)) == 0 - ]) - #We search for stationary points from offset hours before t0 to - #ensure we find any which might occur very soon after t0. - offset = 24.0 - partitions = ( - Tide._times(t0, i*partition) for i in count()), (Tide._times(t0, i*partition) for i in count(1) - ) - - #We'll overestimate to be on the safe side; - #values outside (start,end) won't get yielded. - interval_count = int(np.ceil((partition + offset) / delta)) + 1 - amplitude = self.model['amplitude'][:, np.newaxis] - phase = d2r*self.model['phase'][:, np.newaxis] - - for start, end in izip(*partitions): - speed, [u], [f], V0 = self.prepare(start, Tide._times(start, 0.5*partition)) - #These derivatives don't include the time dependence of u or f, - #but these change slowly. - def d(t): - return np.sum(-speed*amplitude*f*np.sin(speed*t + (V0 + u) - phase), axis=0) - def d2(t): - return np.sum(-speed**2.0 * amplitude*f*np.cos(speed*t + (V0 + u) - phase), axis=0) - #We'll overestimate to be on the safe side; - #values outside (start,end) won't get yielded. - intervals = ( - delta * i -offset for i in range(interval_count)), (delta*(i+1) - offset for i in range(interval_count) - ) - for a, b in izip(*intervals): - if d(a)*d(b) < 0: - extrema = fsolve(d, (a + b) / 2.0, fprime = d2)[0] - time = Tide._times(start, extrema) - [height] = self.at([time]) - hilo = 'H' if d2(extrema) < 0 else 'L' - if start < time < end: - yield (time, height, hilo) - - @staticmethod - def _hours(t0, t): - """ - Return the hourly offset(s) of a (list of) time from a given time. - Arguments: - t0 -- time from which offsets are sought - t -- times to find hourly offsets from t0. - """ - if not isinstance(t, Iterable): - return Tide._hours(t0, [t])[0] - elif isinstance(t[0], datetime): - return np.array([(ti-t0).total_seconds() / 3600.0 for ti in t]) - else: - return t - - @staticmethod - def _partition(hours, partition = 3600.0): - """ - Partition a sorted list of numbers (or in this case hours). - Arguments: - hours -- sorted ndarray of hours. - partition -- maximum partition length (default: 3600.0) - """ - partition = float(partition) - relative = hours - hours[0] - total_partitions = np.ceil(relative[-1] / partition + 10*np.finfo(np.float).eps).astype('int') - return [hours[np.floor(np.divide(relative, partition)) == i] for i in range(total_partitions)] - - @staticmethod - def _times(t0, hours): - """ - Return a (list of) datetime(s) given an initial time and an (list of) hourly offset(s). - Arguments: - t0 -- initial time - hours -- hourly offsets from t0 - """ - if not isinstance(hours, Iterable): - return Tide._times(t0, [hours])[0] - elif not isinstance(hours[0], datetime): - return np.array([t0 + timedelta(hours=h) for h in hours]) - else: - return np.array(hours) - - @staticmethod - def _tidal_series(t, amplitude, phase, speed, u, f, V0): - return np.sum(amplitude*f*np.cos(speed*t + (V0 + u) - phase), axis=0) - - def normalize(self): - """ - Adapt self.model so that amplitudes are positive and phases are in [0,360) as per convention - """ - for i, (_, amplitude, phase) in enumerate(self.model): - if amplitude < 0: - self.model['amplitude'][i] = -amplitude - self.model['phase'][i] = phase + 180.0 - self.model['phase'][i] = np.mod(self.model['phase'][i], 360.0) - - @classmethod - def decompose( - cls, - heights, - t = None, - t0 = None, - interval = None, - constituents = constituent.noaa, - initial = None, - n_period = 2, - callback = None, - full_output = False - ): - """ - Return an instance of Tide which has been fitted to a series of tidal observations. - Arguments: - It is not necessary to provide t0 or interval if t is provided. - heights -- ndarray of tidal observation heights - t -- ndarray of tidal observation times - t0 -- datetime representing the time at which heights[0] was recorded - interval -- hourly interval between readings - constituents -- list of constituents to use in the fit (default: constituent.noaa) - initial -- optional Tide instance to use as first guess for least squares solver - n_period -- only include constituents which complete at least this many periods (default: 2) - callback -- optional function to be called at each iteration of the solver - full_output -- whether to return the output of scipy's leastsq solver (default: False) - """ - if t is not None: - if isinstance(t[0], datetime): - hours = Tide._hours(t[0], t) - t0 = t[0] - elif t0 is not None: - hours = t - else: - raise ValueError("t can be an array of datetimes, or an array " - "of hours since t0 in which case t0 must be " - "specified.") - elif None not in [t0, interval]: - hours = np.arange(len(heights)) * interval - else: - raise ValueError("Must provide t(datetimes), or t(hours) and " - "t0(datetime), or interval(hours) and t0(datetime) " - "so that each height can be identified with an " - "instant in time.") - - #Remove duplicate constituents (those which travel at exactly the same - #speed, irrespective of phase) - constituents = list(OrderedDict.fromkeys(constituents)) - - #No need for least squares to find the mean water level constituent z0, - #work relative to mean - constituents = [c for c in constituents if not c == constituent._Z0] - z0 = np.mean(heights) - heights = heights - z0 - - #Only analyse frequencies which complete at least n_period cycles over - #the data period. - constituents = [ - c for c in constituents - if 360.0 * n_period < hours[-1] * c.speed(astro(t0)) - ] - n = len(constituents) - - sort = np.argsort(hours) - hours = hours[sort] - heights = heights[sort] - - #We partition our time/height data into intervals over which we consider - #the values of u and f to assume a constant value (that is, their true - #value at the midpoint of the interval). Constituent - #speeds change much more slowly than the node factors, so we will - #consider these constant and equal to their speed at t0, regardless of - #the length of the time series. - - partition = 240.0 - - t = Tide._partition(hours, partition) - times = Tide._times(t0, [(i + 0.5)*partition for i in range(len(t))]) - - speed, u, f, V0 = Tide._prepare(constituents, t0, times, radians = True) - - #Residual to be minimised by variation of parameters (amplitudes, phases) - def residual(hp): - H, p = hp[:n, np.newaxis], hp[n:, np.newaxis] - s = np.concatenate([ - Tide._tidal_series(t_i, H, p, speed, u_i, f_i, V0) - for t_i, u_i, f_i in izip(t, u, f) - ]) - res = heights - s - if callback: - callback(res) - return res - - #Analytic Jacobian of the residual - this makes solving significantly - #faster than just using gradient approximation, especially with many - #measurements / constituents. - def D_residual(hp): - H, p = hp[:n, np.newaxis], hp[n:, np.newaxis] - ds_dH = np.concatenate([ - f_i*np.cos(speed*t_i+u_i+V0-p) - for t_i, u_i, f_i in izip(t, u, f)], - axis = 1) - - ds_dp = np.concatenate([ - H*f_i*np.sin(speed*t_i+u_i+V0-p) - for t_i, u_i, f_i in izip(t, u, f)], - axis = 1) - - return np.append(-ds_dH, -ds_dp, axis=0) - - #Initial guess for solver, haven't done any analysis on this since the - #solver seems to converge well regardless of the initial guess We do - #however scale the initial amplitude guess with some measure of the - #variation - amplitudes = np.ones(n) * (np.sqrt(np.dot(heights, heights)) / len(heights)) - phases = np.ones(n) - - if initial: - for (c0, amplitude, phase) in initial.model: - for i, c in enumerate(constituents): - if c0 == c: - amplitudes[i] = amplitude - phases[i] = d2r*phase - - initial = np.append(amplitudes, phases) - - lsq = leastsq(residual, initial, Dfun=D_residual, col_deriv=True, ftol=1e-7) - - model = np.zeros(1+n, dtype=cls.dtype) - model[0] = (constituent._Z0, z0, 0) - model[1:]['constituent'] = constituents[:] - model[1:]['amplitude'] = lsq[0][:n] - model[1:]['phase'] = lsq[0][n:] - - if full_output: - return cls(model = model, radians = True), lsq - return cls(model = model, radians = True) + dtype = np.dtype([ + ('constituent', object), + ('amplitude', float), + ('phase', float)]) + + def __init__( + self, + constituents=None, + amplitudes=None, + phases=None, + model=None, + radians=False): + """ + Initialise a tidal model. Provide constituents, amplitudes and phases + OR a model. + Arguments: + constituents -- list of constituents used in the model. + amplitudes -- list of amplitudes corresponding to constituents + phases -- list of phases corresponding to constituents + model -- an ndarray of type Tide.dtype representing the constituents, + amplitudes and phases. + radians -- boolean representing whether phases are in radians + (default False) + """ + if None not in [constituents, amplitudes, phases]: + if len(constituents) == len(amplitudes) == len(phases): + model = np.zeros(len(phases), dtype=Tide.dtype) + model['constituent'] = np.array(constituents) + model['amplitude'] = np.array(amplitudes) + model['phase'] = np.array(phases) + else: + raise ValueError("Constituents, amplitudes and phases should " + "all be arrays of equal length.") + elif model is not None: + if not model.dtype == Tide.dtype: + raise ValueError("Model must be a numpy array with " + "dtype == Tide.dtype") + else: + raise ValueError("Must be initialised with constituents, " + "amplitudes and phases; or a model.") + if radians: + model['phase'] = r2d * model['phase'] + self.model = model[:] + self.normalize() + + def prepare(self, *args, **kwargs): + return Tide._prepare(self.model['constituent'], *args, **kwargs) + + @staticmethod + def _prepare(constituents, t0, t=None, radians=True): + """ + Return constituent speed and equilibrium argument at a given time, and + constituent node factors at given times. + Arguments: + constituents -- list of constituents to prepare + t0 -- time at which to evaluate speed and equilibrium argument for each + constituent + t -- list of times at which to evaluate node factors for each + constituent (default: t0) + radians -- whether to return the angular arguments in radians or + degrees (default: True) + """ + # The equilibrium argument is constant and taken at the beginning of + # the time series (t0). The speed of the equilibrium argument changes + # very slowly, so again we take it to be constant over any length of + # data. The node factors change more rapidly. + if isinstance(t0, Iterable): + t0 = t0[0] + if t is None: + t = [t0] + if not isinstance(t, Iterable): + t = [t] + a0 = astro(t0) + a = [astro(t_i) for t_i in t] + + # For convenience give u, V0 (but not speed!) in [0, 360) + V0 = np.array([c.V(a0) for c in constituents])[:, np.newaxis] + speed = np.array([c.speed(a0) for c in constituents])[:, np.newaxis] + u = [np.mod(np.array([c.u(a_i) for c in constituents])[:, np.newaxis], + 360.0) for a_i in a] + f = [np.mod(np.array([c.f(a_i) for c in constituents])[:, np.newaxis], + 360.0) for a_i in a] + + if radians: + speed = d2r * speed + V0 = d2r * V0 + u = [d2r * each for each in u] + return speed, u, f, V0 + + def at(self, t): + """ + Return the modelled tidal height at given times. + Arguments: + t -- array of times at which to evaluate the tidal height + """ + t0 = t[0] + hours = self._hours(t0, t) + partition = 240.0 + t = self._partition(hours, partition) + times = self._times(t0, [(i + 0.5) * partition for i in range(len(t))]) + speed, u, f, V0 = self.prepare(t0, times, radians=True) + H = self.model['amplitude'][:, np.newaxis] + p = d2r * self.model['phase'][:, np.newaxis] + + return np.concatenate([ + Tide._tidal_series(t_i, H, p, speed, u_i, f_i, V0) + for t_i, u_i, f_i in izip(t, u, f) + ]) + + def highs(self, *args): + """ + Generator yielding only the high tides. + Arguments: + see Tide.extrema() + """ + for t in ifilter(lambda e: e[2] == 'H', self.extrema(*args)): + yield t + + def lows(self, *args): + """ + Generator yielding only the low tides. + Arguments: + see Tide.extrema() + """ + for t in ifilter(lambda e: e[2] == 'L', self.extrema(*args)): + yield t + + def form_number(self): + """ + Returns the model's form number, a helpful heuristic for classifying + tides. + """ + k1, o1, m2, s2 = ( + np.extract(self.model['constituent'] == c, self.model['amplitude']) + for c in [constituent._K1, constituent._O1, constituent._M2, + constituent._S2] + ) + return (k1 + o1) / (m2 + s2) + + def classify(self): + """ + Classify the tide according to its form number + """ + form = self.form_number() + if 0 <= form <= 0.25: + return 'semidiurnal' + elif 0.25 < form <= 1.5: + return 'mixed (semidiurnal)' + elif 1.5 < form <= 3.0: + return 'mixed (diurnal)' + else: + return 'diurnal' + + def extrema(self, t0, t1=None, partition=2400.0): + """ + A generator for high and low tides. + Arguments: + t0 -- time after which extrema are sought + t1 -- optional time before which extrema are sought (if not given, the + generator is infinite) + partition -- number of hours for which we consider the node factors to + be constant (default: 2400.0) + """ + if t1: + # yield from in python 3.4 + for e in takewhile(lambda t: t[0] < t1, self.extrema(t0)): + yield e + else: + # We assume that extrema are separated by at least delta hours + delta = np.amin([ + 90.0 / c.speed(astro(t0)) for c in self.model['constituent'] + if not c.speed(astro(t0)) == 0 + ]) + # We search for stationary points from offset hours before t0 to + # ensure we find any which might occur very soon after t0. + offset = 24.0 + partitions = ( + (Tide._times(t0, i * partition) for i in count()), + (Tide._times(t0, i * partition) for i in count(1)) + ) + + # We'll overestimate to be on the safe side; + # values outside (start,end) won't get yielded. + interval_count = int(np.ceil((partition + offset) / delta)) + 1 + amplitude = self.model['amplitude'][:, np.newaxis] + phase = d2r * self.model['phase'][:, np.newaxis] + + for start, end in izip(*partitions): + speed, [u], [f], V0 = self.prepare( + start, Tide._times(start, 0.5 * partition)) + # These derivatives don't include the time dependence of u or + # f, but these change slowly. + + def d(t): + return np.sum( + -speed * amplitude * f + * np.sin( + speed * t + (V0 + u) - phase), + axis=0) + + def d2(t): + return np.sum( + -speed ** 2.0 * amplitude * f + * np.cos( + speed * t + (V0 + u) - phase), + axis=0) + # We'll overestimate to be on the safe side; + # values outside (start,end) won't get yielded. + intervals = ( + (delta * i - offset for i in range(interval_count)), + (delta * (i + 1) - offset for i in range(interval_count)) + ) + for a, b in izip(*intervals): + if d(a) * d(b) < 0: + extrema = fsolve(d, (a + b) / 2.0, fprime=d2)[0] + time = Tide._times(start, extrema) + [height] = self.at([time]) + hilo = 'H' if d2(extrema) < 0 else 'L' + if start < time < end: + yield (time, height, hilo) + + @staticmethod + def _hours(t0, t): + """ + Return the hourly offset(s) of a (list of) time from a given time. + Arguments: + t0 -- time from which offsets are sought + t -- times to find hourly offsets from t0. + """ + if not isinstance(t, Iterable): + return Tide._hours(t0, [t])[0] + elif isinstance(t[0], datetime): + return np.array([(ti-t0).total_seconds() / 3600.0 for ti in t]) + else: + return t + + @staticmethod + def _partition(hours, partition=3600.0): + """ + Partition a sorted list of numbers (or in this case hours). + Arguments: + hours -- sorted ndarray of hours. + partition -- maximum partition length (default: 3600.0) + """ + partition = float(partition) + relative = hours - hours[0] + total_partitions = np.ceil( + relative[-1] / partition + 10 * np.finfo(np.float).eps + ).astype('int') + return [hours[np.floor(np.divide(relative, partition)) == i] + for i in range(total_partitions)] + + @staticmethod + def _times(t0, hours): + """ + Return a (list of) datetime(s) given an initial time and an (list of) + hourly offset(s). + Arguments: + t0 -- initial time + hours -- hourly offsets from t0 + """ + if not isinstance(hours, Iterable): + return Tide._times(t0, [hours])[0] + elif not isinstance(hours[0], datetime): + return np.array([t0 + timedelta(hours=h) for h in hours]) + else: + return np.array(hours) + + @staticmethod + def _tidal_series(t, amplitude, phase, speed, u, f, V0): + return np.sum(amplitude * f * np.cos(speed * t + (V0 + u) - phase), + axis=0) + + def normalize(self): + """ + Adapt self.model so that amplitudes are positive and phases are in + [0,360) as per convention + """ + for i, (_, amplitude, phase) in enumerate(self.model): + if amplitude < 0: + self.model['amplitude'][i] = -amplitude + self.model['phase'][i] = phase + 180.0 + self.model['phase'][i] = np.mod(self.model['phase'][i], 360.0) + + @classmethod + def decompose( + cls, + heights, + t=None, + t0=None, + interval=None, + constituents=constituent.noaa, + initial=None, + n_period=2, + callback=None, + full_output=False + ): + """ + Return an instance of Tide which has been fitted to a series of tidal + observations. + Arguments: + It is not necessary to provide t0 or interval if t is provided. + heights -- ndarray of tidal observation heights + t -- ndarray of tidal observation times + t0 -- datetime representing the time at which heights[0] was recorded + interval -- hourly interval between readings + constituents -- list of constituents to use in the fit + (default: constituent.noaa) + initial -- optional Tide instance to use as first guess for least + squares solver + n_period -- only include constituents which complete at least this many + periods (default: 2) + callback -- optional function to be called at each iteration of the + solver + full_output -- whether to return the output of scipy's leastsq solver + (default: False) + """ + if t is not None: + if isinstance(t[0], datetime): + hours = Tide._hours(t[0], t) + t0 = t[0] + elif t0 is not None: + hours = t + else: + raise ValueError("t can be an array of datetimes, or an array " + "of hours since t0 in which case t0 must be " + "specified.") + elif None not in [t0, interval]: + hours = np.arange(len(heights)) * interval + else: + raise ValueError("Must provide t(datetimes), or t(hours) and " + "t0(datetime), or interval(hours) and " + "t0(datetime) so that each height can be " + "identified with an instant in time.") + + # Remove duplicate constituents (those which travel at exactly the same + # speed, irrespective of phase) + constituents = list(OrderedDict.fromkeys(constituents)) + + # No need for least squares to find the mean water level constituent + # z0, work relative to mean + constituents = [c for c in constituents if not c == constituent._Z0] + z0 = np.mean(heights) + heights = heights - z0 + + # Only analyse frequencies which complete at least n_period cycles over + # the data period. + constituents = [ + c for c in constituents + if 360.0 * n_period < hours[-1] * c.speed(astro(t0)) + ] + n = len(constituents) + + sort = np.argsort(hours) + hours = hours[sort] + heights = heights[sort] + + # We partition our time/height data into intervals over which we + # consider the values of u and f to assume a constant value (that is, + # their true value at the midpoint of the interval). Constituent + # speeds change much more slowly than the node factors, so we will + # consider these constant and equal to their speed at t0, regardless of + # the length of the time series. + + partition = 240.0 + + t = Tide._partition(hours, partition) + times = Tide._times(t0, [(i + 0.5) * partition for i in range(len(t))]) + + speed, u, f, V0 = Tide._prepare(constituents, t0, times, radians=True) + + # Residual to be minimised by variation of parameters + # (amplitudes, phases) + def residual(hp): + H, p = hp[:n, np.newaxis], hp[n:, np.newaxis] + s = np.concatenate([ + Tide._tidal_series(t_i, H, p, speed, u_i, f_i, V0) + for t_i, u_i, f_i in izip(t, u, f) + ]) + res = heights - s + if callback: + callback(res) + return res + + # Analytic Jacobian of the residual - this makes solving significantly + # faster than just using gradient approximation, especially with many + # measurements / constituents. + def D_residual(hp): + H, p = hp[:n, np.newaxis], hp[n:, np.newaxis] + ds_dH = np.concatenate([ + f_i * np.cos(speed * t_i + u_i + V0 - p) + for t_i, u_i, f_i in izip(t, u, f)], + axis=1) + + ds_dp = np.concatenate([ + H * f_i * np.sin(speed * t_i + u_i + V0 - p) + for t_i, u_i, f_i in izip(t, u, f)], + axis=1) + + return np.append(-ds_dH, -ds_dp, axis=0) + + # Initial guess for solver, haven't done any analysis on this since the + # solver seems to converge well regardless of the initial guess We do + # however scale the initial amplitude guess with some measure of the + # variation + amplitudes = np.ones(n) * ( + np.sqrt(np.dot(heights, heights)) / len(heights)) + phases = np.ones(n) + + if initial: + for (c0, amplitude, phase) in initial.model: + for i, c in enumerate(constituents): + if c0 == c: + amplitudes[i] = amplitude + phases[i] = d2r * phase + + initial = np.append(amplitudes, phases) + + lsq = leastsq(residual, initial, Dfun=D_residual, col_deriv=True, + ftol=1e-7) + + model = np.zeros(1 + n, dtype=cls.dtype) + model[0] = (constituent._Z0, z0, 0) + model[1:]['constituent'] = constituents[:] + model[1:]['amplitude'] = lsq[0][:n] + model[1:]['phase'] = lsq[0][n:] + + if full_output: + return cls(model=model, radians=True), lsq + return cls(model=model, radians=True) diff --git a/requirements.txt b/requirements.txt new file mode 100644 index 0000000..ed284ab --- /dev/null +++ b/requirements.txt @@ -0,0 +1,2 @@ +-r requirements/_common.txt +-r requirements/_numpy_and_scipy.txt diff --git a/requirements/_common.txt b/requirements/_common.txt new file mode 100644 index 0000000..cf45f82 --- /dev/null +++ b/requirements/_common.txt @@ -0,0 +1 @@ +# any other requirements in here diff --git a/requirements/_numpy_and_scipy.txt b/requirements/_numpy_and_scipy.txt new file mode 100644 index 0000000..2092ca1 --- /dev/null +++ b/requirements/_numpy_and_scipy.txt @@ -0,0 +1,5 @@ +# These are broken out because Travis CI needs to use the system-wide packages +# rather than installing them into the virtualenv. + +scipy==0.14.0 +numpy==1.8.1 diff --git a/requirements_for_travis.txt b/requirements_for_travis.txt new file mode 100644 index 0000000..55f4716 --- /dev/null +++ b/requirements_for_travis.txt @@ -0,0 +1,5 @@ +-r requirements/_common.txt + +nose==1.3.3 +flake8==2.1.0 +pyflakes==0.8.1 diff --git a/run_tests.sh b/run_tests.sh new file mode 100755 index 0000000..569a283 --- /dev/null +++ b/run_tests.sh @@ -0,0 +1,13 @@ +#!/bin/bash -e + +function run_pep8_style_checks { + flake8 . +} + + +function run_python_unit_tests { + nosetests +} + +# run_pep8_style_checks # disabled until PEP8 changes are merged +run_python_unit_tests