Source code for oceans.RPSstuff

from datetime import datetime

import numpy as np
import numpy.ma as ma

earth_radius = 6371.0e3


[docs] def h2hms(hours): """ Converts hours to hours, minutes, and seconds. Examples -------- >>> h2hms(12.51) (np.float64(12.0), np.float64(30.0), np.float64(36.0)) """ hour = np.floor(hours) mins = np.remainder(hours, 1.0) * 60.0 mn = np.floor(mins) secs = np.round(np.remainder(mins, 1.0) * 60.0) return hour, mn, secs
[docs] def hms2h(h, m=None, s=None): """ Converts hours, minutes, and seconds to hours. Examples -------- >>> hms2h(12.0, 30.0, 36.0) 12.51 >>> # Or, >>> hms2h(123036) np.float64(12.51) """ if not m and not s: hms = h h = np.floor(hms / 10000) ms = hms - h * 10000 m = np.floor(ms / 100) s = ms - m * 100 hours = h + m / 60 + s / 3600 else: hours = h + (m + s / 60) / 60 return hours
[docs] def ms2hms(millisecs): """ Converts milliseconds to integer hour, minute, seconds. Examples -------- >>> ms2hms(1e3 * 60) (np.float64(0.0), np.float64(1.0), np.float64(0.0)) """ sec = np.round(millisecs / 1000) hour = np.floor(sec / 3600) mn = np.floor(np.remainder(sec, 3600) / 60) sec = np.round(np.remainder(sec, 60)) return hour, mn, sec
[docs] def julian(y, m=0, d=0, h=0, mi=0, s=0, noon=False): """ Converts Gregorian calendar dates to Julian dates USAGE: [j]=julian(y,m,d,h) DESCRIPTION: Converts Gregorian dates to decimal Julian days using the astronomical conversion, but with time zero starting at midnight instead of noon. In this convention, Julian day 2440000 begins at 0000 hours, May 23, 1968. The decimal Julian day, with double precision, yields an accuracy of decimal days of about 0.1 milliseconds. If you want Julian days to start and end at noon set `noon` to True. INPUT: y : year (e.g., 1979) component m : month (1-12) component d : day (1-31) component of Gregorian date hour : hours (0-23) min : minutes (0-59) sec : decimal seconds or h : decimal hours (assumed 0 if absent) OUTPUT: j : decimal Julian day number last revised 1/3/96 by Rich Signell (rsignell@usgs.gov) Examples -------- >>> julian(1968, 5, 23, 0) array([2440000.]) """ y, m, d, h, mi, s = list(map(np.atleast_1d, (y, m, d, h, mi, s))) h = hms2h(h, mi, s) mo = m + 9 yr = y - 1 i = m > 2 mo[i] = m[i] - 3 yr[i] = y[i] c = np.floor(yr / 100.0) yr = yr - c * 100 j = ( np.floor((146097 * c) / 4.0) + np.floor((1461 * yr) / 4.0) + np.floor((153 * mo + 2) / 5.0) + d + 1721119 ) if noon: j = j + (h - 12) / 24 else: j = j + h / 24 return j
[docs] def jdrps2jdmat(jd): """ Convert Signell's Julian days to Matlab's Serial day matlab's serial date = 1 at 0000 UTC, 1-Jan-0000. Examples -------- >>> jdrps2jdmat(2440000) array([718941.]) """ return jd - julian(0000, 1, 1, 0, 0, 0) + 1
[docs] def jdmat2jdrps(jdmat): """ Convert Matlab's Serial Day to Signell's Julian days matlab's serial date = 1 at 0000 UTC, 1-Jan-0000. Examples -------- >>> jdmat2jdrps(718941) array([2440000.]) """ return jdmat + julian(0000, 1, 1, 0, 0, 0) - 1
[docs] def gregorian(jd, noon=False): """ Converts decimal Julian days to Gregorian dates using the astronomical conversion, but with time zero starting at midnight instead of noon. In this convention, Julian day 2440000 begins at 0000 hours, May 23, 1968. The Julian day does not have to be an integer, and with Matlab's double precision, the accuracy of decimal days is about 0.1 milliseconds. INPUT: jd = decimal Julian days OUTPUT: gtime = six column Gregorian time matrix, where each row is: [yyyy mo da hr mi sec]. yyyy = year (e.g., 1979) mo = month (1-12) da = day (1-31) hr = hour (0-23) mi = minute (0-59) sec = decimal seconds example: [1990 12 12 0 0 0] is midnight on Dec 12, 1990. Examples -------- >>> gregorian(2440000) array([[1968., 5., 23., 0., 0., 0.]]) AUTHOR: Rich Signell (rsignell@usgs.gov) """ # Add 0.2 milliseconds before Gregorian calculation to prevent # roundoff error resulting from math operations on time # from occasionally representing midnight as # (for example) [1990 11 30 23 59 59.99...] instead of [1990 12 1 0 0 0]); # If adding a 0.2 ms to time (each time you go back and forth between # Julian and Gregorian) bothers you more than the inconvenient # representation of Gregorian time at midnight you can comment this # line out... jd = np.atleast_1d(jd) jd = jd + 2.0e-9 if noon: h = np.remainder(jd, 1) * 24 + 12 i = h >= 24 jd[i] = jd[i] + 1 h[i] = h[i] - 24 secs = np.remainder(jd, 1) * 24 * 3600 j = np.floor(jd) - 1721119 ini = 4 * j - 1 y = np.floor(ini / 146097) j = ini - 146097 * y ini = np.floor(j / 4) ini = 4 * ini + 3 j = np.floor(ini / 1461) d = np.floor(((ini - 1461 * j) + 4) / 4) ini = 5 * d - 3 m = np.floor(ini / 153) d = np.floor(((ini - 153 * m) + 5) / 5) y = y * 100 + j mo = m - 9 yr = y + 1 i = m < 10 mo[i] = m[i] + 3 yr[i] = y[i] hr, mi, sc = s2hms(secs) return np.c_[yr, mo, d, hr, mi, sc]
[docs] def s2hms(secs): """ Converts seconds to integer hour,minute,seconds Usage: hour, min, sec = s2hms(secs) Examples -------- >>> s2hms(3600 + 60 + 1) (np.float64(1.0), np.float64(1.0), np.int64(1)) """ hr = np.floor(secs / 3600) mi = np.floor(np.remainder(secs, 3600) / 60) sc = np.round(np.remainder(secs, 60)) return hr, mi, sc
# FIXME: STOPPED HERE
[docs] def ss2(jd): """ Return Gregorian start and stop dates of Julian day variable Usage: start, stop = ss2(jd) """ start = gregorian(jd[0]) stop = gregorian(jd[-1]) return start, stop
[docs] def angled(h): """ ANGLED: Returns the phase angles in degrees of a matrix with complex elements. Usage: deg = angled(h) h = complex matrix deg = angle in math convention (degrees counterclockwise from "east") """ pd = np.angle(h, deg=True) return pd
[docs] def ij2ind(a, i, j): m, n = a.shape return m * i - j + 1 # TODO: Check this +1
[docs] def ind2ij(a, ind): """ ind2ij returns i, j indices of array. """ m, n = a.shape j = np.ceil(ind / m) i = np.remainder(ind, m) i[i == 0] = m return i, j
[docs] def rms(u): """ Compute root mean square for each column of matrix u. """ # TODO: use an axis arg. if u.ndim > 1: m, n = u.shape else: m = u.size return np.sqrt(np.sum(u**2) / m)
[docs] def z0toCn(z0, H): """ Convert roughness height z0 to Chezy "C" and Manning's "n" which is a function of the water depth Inputs: z0 = roughness height (meters) H = water depth (meters) (can be vector) Outputs: C = Chezy "C" (non-dimensional) n = Manning's "n" (non-dimensional) Examples -------- >>> # finds vectors C and n corresponding to a z0=0.003 >>> # and a range of water depths from 2--200 meters. >>> C, n = z0toCn(0.003, np.arange(2, 200)) """ k_s = 30 * z0 C = 18 * np.log10(12 * H / k_s) n = (H ** (1.0 / 6.0)) / C return C, n
[docs] def z0tocd(z0, zr): """ Calculates CD at a given ZR corresponding to Z0. """ cd = (0.4 * np.ones(z0.size) / np.log(zr / z0)) ** 2 return cd
[docs] def short_calc(amin, amax): rang = 32767 - (-32767) add_offset = (amax + amin) * 0.5 scale_factor = (amax - amin) / rang return add_offset, scale_factor
[docs] def gsum(x, **kw): """ Just like sum, except that it skips over bad points. """ xnew = ma.masked_invalid(x) return np.sum(xnew, **kw)
[docs] def gmean(x, **kw): """ Just like mean, except that it skips over bad points. """ xnew = ma.masked_invalid(x) return np.mean(xnew, **kw)
[docs] def gmedian(x, **kw): """ Just like median, except that it skips over bad points. """ xnew = ma.masked_invalid(x) return np.median(xnew, **kw)
[docs] def gmin(x, **kw): """ Just like min, except that it skips over bad points. """ xnew = ma.masked_invalid(x) return np.min(xnew, **kw)
[docs] def gmax(x, **kw): """ Just like max, except that it skips over bad points. """ xnew = ma.masked_invalid(x) return np.max(xnew, **kw)
[docs] def gstd(x, **kw): """ Just like std, except that it skips over bad points. """ xnew = ma.masked_invalid(x) return np.std(xnew, **kw)
[docs] def near(x, x0, n=1): """ Given an 1D array `x` and a scalar `x0`, returns the `n` indices of the element of `x` closest to `x0`. """ distance = np.abs(x - x0) index = np.argsort(distance) return index[:n], distance[index[:n]]
[docs] def swantime(a): """ Converts SWAN default time format to datetime object. """ if isinstance(a, str): a = float(a) a = np.asanyarray(a) year = np.floor(a / 1e4) a = a - year * 1e4 mon = np.floor(a / 1e2) a = a - mon * 1e2 day = np.floor(a) a = a - day hour = np.floor(a * 1e2) a = a - hour / 1e2 mn = np.floor(a * 1e4) a = a - mn / 1e4 sec = np.floor(a * 1e6) return datetime(year, mon, day, hour, mn, sec)
[docs] def shift(a, b, n): """ a and b are vectors n is number of points of a to cut off anew and bnew will be the same length. """ # la, lb = a.size, lb = b.size anew = a[list(range(0 + n, len(a))), :] if len(anew) > len(b): anew = anew[list(range(0, len(b))), :] bnew = b else: bnew = b[list(range(0, len(anew))), :] return anew, bnew
[docs] def lagcor(a, b, n): """ Finds lagged correlations between two series. a and b are two column vectors n is range of lags cor is correlation as fn of lag. """ cor = [] for k in range(0, n + 1): d1, d2 = shift(a, b, k) ind = ~np.isnan(d1 + d2) c = np.corrcoef(d1[ind], d2[ind]) if len(c) > 1: cor.append(c[0, 1]) return cor
[docs] def coast2bln(coast, bln_file): """ Converts a matlab coast (two column array w/ nan for line breaks) into a Surfer blanking file. Where coast is a two column vector and bln_file is the output file name. Needs `fixcoast`. """ c2 = fixcoast(coast) ind = np.where(np.isnan(c2[:, 0]))[0] n = len(ind) - 1 bln = c2.copy() for k in range(0, n - 1): kk = list(range(ind[k] + 1, ind[k + 1])) NP = int(len(kk)) bln[ind[k], 0] = NP bln[ind[k], 1] = 1 bln = bln[:-1] np.savetxt(bln_file, bln, fmt="%g")
[docs] def fixcoast(coast): """ FIXCOAST Makes sure coastlines meet Signell's conventions. Fixes coastline is in the format we want. Assumes that lon/lat are in the first two columns of the matrix coast, and that coastline segments are separated by rows of NaNs (or -99999s). This routine ensures that only 1 row of NaNs separates each segment, and makes sure that the first and last rows contain NaNs. """ ind = coast == -99999.0 coast[ind] = np.nan ind = np.where(np.isnan(coast[:, 0]))[0] dind = np.diff(ind) idup = np.where(dind == 1)[0] coast = np.delete(coast, ind[idup], axis=0) if not np.isnan(coast[0, 0]): coast = np.insert(coast, 0, np.nan, axis=0) if not np.isnan(coast[-1, -1]): coast = np.append(coast, np.c_[np.nan, np.nan], axis=0) return coast