diff --git a/seaexplorertools/BYQtools.py b/seaexplorertools/BYQtools.py deleted file mode 100644 index 2cad64e..0000000 --- a/seaexplorertools/BYQtools.py +++ /dev/null @@ -1,368 +0,0 @@ -import os, gzip -from scipy.interpolate import interp1d, interp2d -from scipy.optimize import fsolve, fmin -from scipy import signal -import pandas as pd -import numpy as np -import gsw -from datetime import datetime as dt - -### GENERAL SCRIPTS -def RunningMedian(x,N): - grid = np.ones((len(x)+2*N, 1 + 2*N ))*np.NaN - for istep in range(np.shape(grid)[1]): - grid[istep:len(x)+istep, istep] = x - return np.nanmedian(grid,axis=1)[N:-N] - -def RunningMax(x,N): - grid = np.ones((len(x)+2*N, 1 + 2*N ))*np.NaN - for istep in range(np.shape(grid)[1]): - grid[istep:len(x)+istep, istep] = x - return np.nanmax(grid,axis=1)[N:-N] - -def RunningMin(x,N): - grid = np.ones((len(x)+2*N, 1 + 2*N ))*np.NaN - for istep in range(np.shape(grid)[1]): - grid[istep:len(x)+istep, istep] = x - return np.nanmin(grid,axis=1)[N:-N] - -def RunningMean(x,N): - grid = np.ones((len(x)+2*N, 1 + 2*N ))*np.NaN - for istep in range(np.shape(grid)[1]): - grid[istep:len(x)+istep, istep] = x - return np.nanmean(grid,axis=1)[N:-N] - -def interp(x,y,xi): - _gg = np.isfinite(x+y) - return interp1d(x[_gg], y[_gg], bounds_error=False, fill_value=np.NaN)(xi) - -def rmsd(x): - return np.sqrt(np.nanmean(x**2)) - -def plog(msg): - print(str(dt.now().replace(microsecond=0))+' : '+msg) - return None - -### DATA DOWNLOADss -def download_ERA5(years=[2021], months=[1], days=[1], hours=[0], area=[0,0,0,0], variables=None, output_name=None): - """ - Returns xarray dataset. - - years = list of year numbers - months = list of month numbers - days = list of days, from 1 to 31 - hours = list of hours from 0 to 23 - area = list of [N,W,S,E] - variables = default collects heat flux variables, otherwise, specify. - output_name = file to save to if desired - """ - - # Area = [N,W,S,E] - - # First provided by Marcel du Plessis - # Loading ERA5 data directly from the Copernicus Data Service - # The goal of this notebook is to be able to access and analysis ERA5 data using cloud computing. - # This has the obvious advantage of not having to download and store the data on your local computer, - # which can quicly add up to terrabytes if you're looking for long term data. - - # I am following an example from https://towardsdatascience.com/read-era5-directly-into-memory-with-python-511a2740bba0 - - # Variables on the single levels reanalysis can be found here: - # https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-single-levels?tab=overview - - # Let's see how we go: - - if variables is None: - variables = [ - '2m_temperature', - '10m_u_component_of_wind', - '10m_v_component_of_wind', - 'sea_surface_temperature', - 'skin_temperature', - 'surface_pressure', - 'surface_latent_heat_flux', - 'surface_sensible_heat_flux', - 'surface_net_solar_radiation', - 'surface_net_thermal_radiation', - 'total_precipitation', - 'evaporation'] - - import cdsapi - import xarray as xr - from urllib.request import urlopen # start the client - - import certifi - import urllib3 - http = urllib3.PoolManager( - cert_reqs='CERT_REQUIRED', - ca_certs=certifi.where() - ) - - cds = cdsapi.Client() # dataset you want to read - - dataset = "reanalysis-era5-single-levels" # flag to download data - download_flag = True # api parameters - params = { - 'variable': variables, - 'product_type': 'reanalysis', - 'year': [str(x) for x in years], - 'month': [str(format(x, '02')) for x in months], - 'day': [str(format(x, '02')) for x in days], - 'time': [str(format(x, '02'))+':00' for x in hours], - 'area': area, - 'format': 'netcdf' - } - - fl = cds.retrieve(dataset, params) # download the file - - if output_name is not None: - if output_name[-3:] != '.nc': - output_name = output_name+'.nc' - - fl.download(str(output_name)) - - with urlopen(fl.location) as f: - ds = xr.open_dataset(f.read()) - - return ds - -### OPTICAL FUNCTIONS -def austinPetzold_1986(wavelength, k490): - wave = np.array([350, 360, 370, 380, 390, 400, - 410, 420, 430, 440, 450, 460, 470, 480, 490, 500, - 510, 520, 530, 540, 550, 560, 570, 580, 590, 600, - 610, 620, 630, 640, 650, 660, 670, 680, 690, 700 ]) - M = np.array([2.1442, 2.0504, 1.9610, 1.8772, 1.8009, 1.7383, - 1.7591, 1.6974, 1.6108, 1.5169, 1.4158, 1.3077, 1.1982, 1.0955, 1.0000, 0.9118, - 0.8310, 0.7578, 0.6924, 0.6350, 0.5860, 0.5457, 0.5146, 0.4935, 0.4840, 0.4903, - 0.5090, 0.5380, 0.6231, 0.7001, 0.7300, 0.7301, 0.7008, 0.6245, 0.4901, 0.2891 ]) - Kdw = np.array([0.0510, 0.0405, 0.0331, 0.0278, 0.0242, 0.0217, - 0.0200, 0.0189, 0.0182, 0.0178, 0.0176, 0.0176, 0.0179, 0.0193, 0.0224, 0.0280, - 0.0369, 0.0498, 0.0526, 0.0577, 0.0640, 0.0723, 0.0842, 0.1065, 0.1578, 0.2409, - 0.2892, 0.3124, 0.3296, 0.3290, 0.3559, 0.4105, 0.4278, 0.4521, 0.5116, 0.6514 ]) - - for i in np.arange(1,36): - if wave[i] >= wavelength: - l1 = wave[i] - k1 = Kdw[i] - m1 = M[i] - l0 = wave[i-1] - k0 = Kdw[i-1] - m0 = M[i-1] - break - - num = wavelength - l0 - den = l1 - l0 - frac = num / den - - kdiff = k1 - k0 - Kdw_l = k0 + frac*kdiff - - mdiff = m1 - m0 - M_l = m0 + frac*mdiff - - ref = np.argmin(np.abs(wave-490)) - Kd = (M_l/M[ref]) * (k490 - Kdw[ref]) + Kdw_l - return Kd - -def parFraction(wavelength): - """ - Returns PAR fraction assuming boring, standard, generic incoming solar radiation and flat response from PAR sensor. - Units in "per nm". - """ - wavelengths = np.array([400, 412, 443, 490, 510, 555, 625, 670, 700]); - parFraction = np.array([0.0029, 0.0032, 0.0035, 0.0037, 0.0037, 0.0036, 0.0032, 0.0030, 0.0024]); - return interp1d(wavelengths,parFraction, bounds_error=True, fill_value=np.NaN)(wavelength) - -def morelCoefficients(wavelength): - """ - input: wavelength - output: Kw, eChl, XChl - - link: 'http://onlinelibrary.wiley.com/doi/10.1029/2000JC000319/epdf' - title: 'Bio-optical properties of oceanic waters: A reappraisal' - author: 'André Morel, Stéphane Maritorena' - doi: '10.1029/2000JC000319' - algorithm: 'Kd(lambda) = Kw(lambda) + XChl(lambda) * [Chl]^(eChl(lambda))' - """ - wavelengths = np.array([350,355,360,365,370,375,380,385,390,395, - 400,405,410,415,420,425,430,435,440,445, - 450,455,460,465,470,475,480,485,490,495, - 500,505,510,515,520,525,530,535,540,545, - 550,555,560,565,570,575,580,585,590,595, - 600,605,610,615,620,625,630,635,640,645, - 650,655,660,665,670,675,680,685,690,695,700]) - - interp = lambda arr: interp1d(wavelengths, arr, bounds_error=True, fill_value=np.NaN)(wavelength) - - Kw = np.array([0.0271,0.0238,0.0216,0.0188,0.0177,0.01595,0.0151,0.01376,0.01271,0.01208, - 0.01042,0.0089,0.00812,0.00765,0.00758,0.00768,0.0077,0.00792,0.00885,0.0099, - 0.01148,0.01182,0.01188,0.01211,0.01251,0.0132,0.01444,0.01526,0.0166,0.01885, - 0.02188,0.02701,0.03385,0.0409,0.04214,0.04287,0.04454,0.0463,0.04846,0.05212, - 0.05746,0.06053,0.0628,0.06507,0.07034,0.07801,0.09038,0.11076,0.13584,0.16792, - 0.2231,0.25838,0.26506,0.26843,0.27612,0.284,0.29218,0.30176,0.31134,0.32553, - 0.34052,0.3715,0.41048,0.42947,0.43946,0.44844,0.46543,0.48642,0.5164,0.55939,0.62438]) - - eChl = np.array([0.778,0.767,0.756,0.737,0.72,0.7,0.685,0.673,0.67,0.66,0.64358,0.64776, - 0.65175,0.65555,0.65917,0.66259,0.66583,0.66889,0.67175,0.67443,0.67692, - 0.67923,0.68134,0.68327,0.68501,0.68657,0.68794,0.68903,0.68955,0.68947, - 0.6888,0.68753,0.68567,0.6832,0.68015,0.67649,0.67224,0.66739,0.66195,0.65591, - 0.64927,0.64204,0.64,0.63,0.623,0.615,0.61,0.614,0.618,0.622,0.626,0.63, - 0.634,0.638,0.642,0.647,0.653,0.658,0.663,0.667,0.672,0.677,0.682,0.687,0.695, - 0.697,0.693,0.665,0.64,0.62,0.6]) - - XChl = np.array([0.153,0.149,0.144,0.14,0.136,0.131,0.127,0.123,0.119,0.118,0.11748,0.12066, - 0.12259,0.12326,0.12269,0.12086,0.11779,0.11372,0.10963,0.1056,0.10165, - 0.09776,0.09393,0.09018,0.08649,0.08287,0.07932,0.07584,0.07242,0.06907, - 0.06579,0.06257,0.05943,0.05635,0.05341,0.05072,0.04829,0.04611,0.04419, - 0.04253,0.04111,0.03996,0.039,0.0375,0.036,0.034,0.033,0.0328,0.0325, - 0.033,0.034,0.035,0.036,0.0375,0.0385,0.04,0.042,0.043,0.044,0.0445, - 0.045,0.046,0.0475,0.049,0.0515,0.052,0.0505,0.044,0.039,0.034,0.03]) - - return interp(Kw), interp(eChl), interp(XChl) - -def betasw_ZHH2009(Tc,S,wavelength=700,theta=117,delta=0.039): - # Xiaodong Zhang, Lianbo Hu, and Ming-Xia He (2009), Scatteirng by pure - # seawater: Effect of salinity, Optics Express, Vol. 17, No. 7, 5698-5710 - # - # wavelength: backscatter wavelength (nm) - # Tc: temperauter in degree Celsius - # S: salinity - # delta: depolarization ratio, if not provided, default = 0.039 will be used. - # theta = beam scattering angle in degrees - # betasw: volume scattering at angles defined by theta. Its size is [x y], - # where x is the number of angles (x = length(theta)) and y is the number - # of wavelengths in wavelength (y = length(wavelength)) - # beta90sw: volume scattering at 90 degree. Its size is [1 y] - # bw: total scattering coefficient. Its size is [1 y] - # for backscattering coefficients, divide total scattering by 2 - # - # Xiaodong Zhang, March 10, 2009 - - def RInw(): - # refractive index of air is from Ciddor (1996,Applied Optics) - n_air = 1.0 + (5792105.0 / (238.0185 - 1 / (wavelength/1e3)**2) + 167917.0 / (57.362 - 1/(wavelength/1e3)**2)) / 1e8 - - # refractive index of seawater is from Quan and Fry (1994, Applied Optics) - n0 = 1.31405 - n1 = 1.779e-4 - n2 = -1.05e-6 - n3 = 1.6e-8 - n4 = -2.02e-6 - n5 = 15.868 - n6 = 0.01155 - n7 = -0.00423 - n8 = -4382 - n9 = 1.1455e6 - - nsw = n0 + (n1 + n2*Tc + n3*Tc**2)*S + n4*Tc**2 + (n5 + n6*S + n7*Tc)/wavelength + n8/wavelength**2 + n9/wavelength**3 # pure seawater - nsw = nsw*n_air - dnswds = (n1 + n2*Tc + n3*Tc**2 + n6/wavelength) * n_air - return nsw, dnswds - - def BetaT(): - # pure water secant bulk Millero (1980, Deep-sea Research) - kw = 19652.21 + 148.4206*Tc - 2.327105*Tc**2 + 1.360477e-2*Tc**3 - 5.155288e-5*Tc**4 - Btw_cal = 1/kw - # isothermal compressibility from Kell sound measurement in pure water - # Btw = (50.88630+0.717582*Tc+0.7819867e-3*Tc**2+31.62214e-6*Tc**3-0.1323594e-6*Tc**4+0.634575e-9*Tc**5)./(1+21.65928e-3*Tc)*1e-6; - # seawater secant bulk - a0 = 54.6746 - 0.603459*Tc + 1.09987e-2*Tc**2 - 6.167e-5*Tc**3 - b0 = 7.944e-2 + 1.6483e-2*Tc - 5.3009e-4*Tc**2 - Ks = kw + a0*S + b0*S**1.5 - - # calculate seawater isothermal compressibility from the secant bulk - IsoComp = 1/Ks*1e-5 # unit is pa - return IsoComp - - def rhou_sw(): - # density of water and seawater,unit is Kg/m^3, from UNESCO,38,1981 - a0 = 8.24493e-1 - a1 = -4.0899e-3 - a2 = 7.6438e-5 - a3 = -8.2467e-7 - a4 = 5.3875e-9 - a5 = -5.72466e-3 - a6 = 1.0227e-4 - a7 = -1.6546e-6 - a8 = 4.8314e-4 - - b0 = 999.842594 - b1 = 6.793952e-2 - b2 = -9.09529e-3 - b3 = 1.001685e-4 - b4 = -1.120083e-6 - b5 = 6.536332e-9 - - # density for pure water - density_w = b0+b1*Tc+b2*Tc**2+b3*Tc**3+b4*Tc**4+b5*Tc**5 - # density for pure seawater - density_sw = density_w +((a0+a1*Tc+a2*Tc**2+a3*Tc**3+a4*Tc**4)*S+(a5+a6*Tc+a7*Tc**2)*S**1.5+a8*S**2) - return density_sw - - - def dlnasw_ds(): - # water activity data of seawater is from Millero and Leung (1976,American - # Journal of Science,276,1035-1077). Table 19 was reproduced using - # Eqs.(14,22,23,88,107) then were fitted to polynominal equation. - # dlnawds is partial derivative of natural logarithm of water activity - # w.r.t.salinity - # lnaw = (-1.64555e-6-1.34779e-7*Tc+1.85392e-9*Tc**2-1.40702e-11*Tc**3)+...... - # (-5.58651e-4+2.40452e-7*Tc-3.12165e-9*Tc**2+2.40808e-11*Tc**3)*S+...... - # (1.79613e-5-9.9422e-8*Tc+2.08919e-9*Tc**2-1.39872e-11*Tc**3)*S**1.5+...... - # (-2.31065e-6-1.37674e-9*Tc-1.93316e-11*Tc**2)*S**2; - - dlnawds = (-5.58651e-4+2.40452e-7*Tc-3.12165e-9*Tc**2+2.40808e-11*Tc**3) + 1.5*(1.79613e-5-9.9422e-8*Tc+2.08919e-9*Tc**2-1.39872e-11*Tc**3)*S**0.5 + 2*(-2.31065e-6-1.37674e-9*Tc-1.93316e-11*Tc**2)*S - return dlnawds - - - # density derivative of refractive index from PMH model - def PMH(n_wat): - n_wat2 = n_wat**2 - n_density_derivative = (n_wat2 - 1) * (1+2/3*(n_wat2+2) * (n_wat/3-1/3/n_wat)**2) - return n_density_derivative - - # values of the constants - Na = 6.0221417930e23 # Avogadro's constant - Kbz = 1.3806503e-23 # Boltzmann constant - Tk = Tc+273.15 # Absolute tempearture - M0 = 18e-3 # Molecular weigth of water in kg/mol - - rad = theta*np.pi/180 # angle in radian as a colum variable - - # nsw: absolute refractive index of seawater - # dnds: partial derivative of seawater refractive index w.r.t. salinity - nsw, dnds = RInw() - - # isothermal compressibility is from Lepple & Millero (1971,Deep - # Sea-Research), pages 10-11 - # The error ~ +/-0.004e-6 bar^-1 - IsoComp = BetaT() - - # density of water and seawater,unit is Kg/m^3, from UNESCO,38,1981 - density_sw = rhou_sw() - - # water activity data of seawater is from Millero and Leung (1976,American - # Journal of Science,276,1035-1077). Table 19 was reproduced using - # Eq.(14,22,23,88,107) then were fitted to polynominal equation. - # dlnawds is partial derivative of natural logarithm of water activity - # w.r.t.salinity - dlnawds = dlnasw_ds() - - # density derivative of refractive index from PMH model - DFRI = PMH(nsw) ## PMH model - - # volume scattering at 90 degree due to the density fluctuation - beta_df = np.pi*np.pi/2*((wavelength*1e-9)**(-4))*Kbz*Tk*IsoComp*DFRI**2*(6+6*delta)/(6-7*delta) - - # volume scattering at 90 degree due to the concentration fluctuation - flu_con = S*M0*dnds**2/density_sw/(-dlnawds)/Na - beta_cf = 2*np.pi*np.pi*((wavelength*1e-9)**(-4))*nsw**2*(flu_con)*(6+6*delta)/(6-7*delta) - - # total volume scattering at 90 degree - beta90sw = beta_df+beta_cf - bsw=8*np.pi/3*beta90sw*(2+delta)/(1+delta) - - betasw = beta90sw * (1+((np.cos(rad))**2)*(1-delta)/(1+delta)) - - return betasw,beta90sw,bsw \ No newline at end of file diff --git a/seaexplorertools/SXBQ.py b/seaexplorertools/SXBQ.py index 4b90d61..51f544c 100644 --- a/seaexplorertools/SXBQ.py +++ b/seaexplorertools/SXBQ.py @@ -1,12 +1,6 @@ -import os, gzip - -from scipy.interpolate import interp1d, interp2d -from scipy.optimize import fsolve, fmin -from scipy import signal +import gzip import pandas as pd import numpy as np -import gsw - from glob import glob from tqdm import tqdm @@ -25,10 +19,10 @@ def __init__(self, *args, gzipped=True): def load_files(self, file_dir, gzipped=True): ''' Load payload data into one dataframe and convert main timestamp to datetime format. ''' - + file_list = glob(file_dir, recursive=True) #print(file_list) - + def extract_gzipped(filename): f = gzip.open(fileName) try: @@ -49,12 +43,12 @@ def extract_plaintext(filename): print('Failed to load : '+filename) dive = pd.DataFrame() return dive - + if gzipped: extract_fn = extract_gzipped else: extract_fn = extract_plaintext - + _tmp = [] for fileName in tqdm(file_list): dive = extract_fn(fileName) @@ -70,13 +64,13 @@ def extract_plaintext(filename): dive['Timestamp'] = dive['timeindex'].values.astype("float") dive.set_index('timeindex', inplace=True) _tmp.append(dive[dive.index > '2020-01-01'].resample('S').mean()) - + for d in range(len(_tmp)): _tmp[d]['Timestamp'] = pd.to_datetime(_tmp[d]['Timestamp'].interpolate('linear'), utc=True, origin='unix', cache='False') - + self.data = self.data.append(pd.concat(_tmp, ignore_index=True), sort=True) self.data.sort_values('Timestamp', ignore_index=True, inplace=True) - + def save(self, file_name): if file_name: print('Saving to '+file_name) @@ -88,7 +82,7 @@ def median_resample(self): self.data = self.data.resample('S').mean() #self.data = self.data.resample('S').median() self.data['Timestamp'] = pd.to_datetime(self.data['Timestamp'].interpolate('linear'), utc=True, origin='unix', cache='False') - + def process_basic_variables(self): # TODO: move basic parsing from ipynb to here if ('Lon' in self.data.columns) and ('Lat' in self.data.columns) and ('DeadReckoning' in self.data.columns): @@ -104,7 +98,7 @@ def process_basic_variables(self): # Include printed output so user know what is going on. - + ########################################################################################################################################################## ########################################################################################################################################################## ########################################################################################################################################################## @@ -123,8 +117,7 @@ def load(parquet_file): def parseGPS(sxGPS): ## CALCULATE SUBSURFACE LAT / LON (TODO: USING DEAD RECKONING) return np.sign(sxGPS) * (np.fix(np.abs(sxGPS)/100) + np.mod(np.abs(sxGPS),100)/60) -def date2float(d, epoch=pd.to_datetime(0, utc=True, origin='unix', cache='False')): - return (d - epoch).dt.total_seconds() + def grid2d(x, y, v, xi=1, yi=1, fn='median'): if np.size(xi) == 1: @@ -135,7 +128,7 @@ def grid2d(x, y, v, xi=1, yi=1, fn='median'): raw = pd.DataFrame({'x':x,'y':y,'v':v}).dropna() grid = np.full([np.size(yi),np.size(xi)], np.nan) - + raw['xbins'],xbin_iter = pd.cut(raw.x, xi,retbins=True,labels=False) raw['ybins'],ybin_iter = pd.cut(raw.y, yi,retbins=True,labels=False) @@ -145,460 +138,3 @@ def grid2d(x, y, v, xi=1, yi=1, fn='median'): XI,YI = np.meshgrid(xi, yi, indexing='ij') return grid,XI.T,YI.T -def applyLag(x,time,N): - try: - time = date2float(time) - print('Applying lag after time conversion:') - except: - print('Applying lag:') - _gg = np.isfinite(time+x) - i_time = np.arange(time[0],time[-1]) - - i_fn = interp1d(time[_gg], x[_gg], bounds_error=False, fill_value=np.NaN) - i_x = i_fn(i_time) - filt = np.concatenate([np.full(N-1,0),np.flip(np.arange(N))**2]) - filt = filt/sum(filt) - - i_x = np.convolve(i_x,filt,mode='same') - i_x[0:N+1] = x[0:N+1] - - _fn = interp1d(i_time, i_x, bounds_error=False, fill_value=np.NaN) - return _fn(time) - -def correctSalinityGarau(data,coefs=None): - if 'speed' in data.data: - flowSpeed = data.data.speed.interpolate('index').values - else: - print('Estimating flow speed...') - flowSpeed = np.abs(np.gradient(data.data.LEGATO_PRESSURE,date2float(data.data.Timestamp)) / np.sin(np.deg2rad(data.data.Pitch.interpolate('index')))).fillna(value=0.0001) - - flowSpeed[flowSpeed > 1] = 1 - flowSpeed[flowSpeed < 0.03] = 0.03 - - time = date2float(data.data['Timestamp']) - temp = data.data['LEGATO_TEMPERATURE'].values - _gd = np.isfinite(time+temp) - data.data['temperature'] = temp - - print('Performing thermal mass correction...') - Fs = np.mean(1/np.gradient(time)) - print(' Assuming a sampling frequency of '+str(Fs)+' Hz.') - - def _calcSal(temp, cond, pres, coefs): - a_offset = coefs[0] - a_slope = coefs[1] - t_offset = coefs[2] - t_slope = coefs[3] - alpha = a_offset + a_slope / flowSpeed - tau = t_offset + t_slope / np.sqrt(flowSpeed) - #tau = 11 # Parameter for Lueck and Picklo (1990) - #alpha = 0.57/tau + 0.03122 # Parameter for Lueck and Picklo (1990) - - alpha[~np.isfinite(alpha)] = a_offset - tau[~np.isfinite(tau)] = t_offset - - beta = 1/tau # Parameter for Lueck and Picklo (1990) - fn = Fs/2 # Nyquist frequency - - a = 4*fn*alpha*tau/(1+4*fn*tau) # Parameter for Lueck and Picklo (1990) - b = 1-2*a/alpha # Parameter for Lueck and Picklo (1990) - - # Compute temperature correction to obtained water temperature within the conductivity cell (Morison, 1994) - _internal_bias = np.full_like(temp,0) - - for sample in np.arange(1,len(_internal_bias)): - # Recursive filter from Morison (1994) - # if np.isfinite(temp[sample-1]): - _internal_bias[sample] = -b[sample] * _internal_bias[sample-1] + a[sample] * (temp[sample] - temp[sample-1]) - return gsw.SP_from_C(cond, temp + _internal_bias, pres) # Practical salinity - - def _regressSal(): - _dives = data.data.diveNum.values[np.isfinite(data.data.diveNum.values)] - _dives = _dives[np.unique(np.linspace(1,len(_dives),30).astype('int') - 1)] - _dives = (np.isin(data.data.diveNum.values, _dives)) & (np.isfinite(data.data.temperature.values)) & (data.data.LEGATO_PRESSURE.values > 7) & (data.data.LEGATO_PRESSURE.values < 250) - _temp = data.data.temperature.values[_dives] - _cond = data.data.LEGATO_CONDUCTIVITY.values[_dives] - _pres = data.data.LEGATO_PRESSURE.values[_dives] - _dnum = data.data.diveNum.values[_dives] ###NEW - - # Use these generic values to scale parameters to same order of magnitude to help fmin. - scaler = np.array([0.0135, 0.0264, 7.1499, 2.7858]) - - T,_,_ = grid2d(data.data.profileNum.values[_dives],data.data.LEGATO_PRESSURE.values[_dives],data.data.LEGATO_TEMPERATURE.values[_dives],xi=1,yi=0.3,fn='mean') - - def _PolyArea_triangles(sal): - S,X,Y = grid2d(data.data.profileNum.values[_dives],data.data.LEGATO_PRESSURE.values[_dives],sal,xi=1,yi=0.3,fn='mean') - def merge_profiles(idx): - x1 = S[:,idx] - x2 = S[:,idx+1] - y1 = T[:,idx] - y2 = T[:,idx+1] - _gd = np.isfinite(x1+x2+y1+y2) - x1 = x1[_gd] - x2 = x2[_gd] - y1 = y1[_gd] - y2 = y2[_gd] - return x1,x2,y1,y2 - def triangleArea(x1,y1,x2,y2,x3,y3): - return abs((0.5)*(x1*(y2-y3)+x2*(y3-y1)+x3*(y1-y2))) - area = 0 - for idd in range(np.shape(T)[1]-1): - x1,x2,y1,y2 = merge_profiles(idd) - area = area + np.nansum(triangleArea( - np.append(x1,np.NaN),np.append(y1,np.NaN), - np.append(x2,np.NaN),np.append(y2,np.NaN), - np.append(np.NaN,x1),np.append(np.NaN,y1) )) - area = area + np.nansum(triangleArea( - np.append(np.NaN,x1),np.append(np.NaN,y1), - np.append(x2,np.NaN),np.append(y2,np.NaN), - np.append(np.NaN,x2),np.append(np.NaN,y2) )) - return area - - def _PolyArea_shoelace(x,y): - _gg = np.isfinite(x+y) - return 0.5*np.abs(np.dot(x[_gg],np.roll(y[_gg],1))-np.dot(y[_gg],np.roll(x[_gg],1))) - - def _scoreFunction(x_vals): - return _PolyArea_triangles(_calcSal(_temp, _cond, _pres, x_vals*scaler)) - -# def _scoreFunction(x_vals): -# return _PolyArea_shoelace( _calcSal(_temp, _cond, _pres, x_vals*scaler) , _temp) - - - print('Beginning regression (slow)...') - print(' Initial minimisation score: '+str(_scoreFunction([1,1,1,1]))+'.') - - with tqdm(total=100) as pbar: - def callbackF(Xi): - pbar.update(1) - R = fmin(_scoreFunction, [1,1,1,1], callback=callbackF, disp=True, full_output=True, maxiter=200, ftol=0.01) - print(' Final minimisation score: '+str(_scoreFunction(R[0]))+'.') - return R[0]*scaler - - if coefs is None: - coefs = _regressSal() - print('Regressed coefficients:') - else: - print('Correcting using supplied coefficients:') - - print(coefs) - print('Applying correction to all data (slow)...') - # We have to interpolate temperature here as it's a recursive filter and NaNs propagate. - # Both temperature and the alpha and tau parameters have to be finite. - # This means interpolate and fill temperature and flowSpeed, - # and avoid divisions by zero (with flowSpeed = 0 or pitch = 0) - data.data['salinity'] = _calcSal( - data.data.temperature.interpolate('index').fillna(method='backfill'), - data.data.LEGATO_CONDUCTIVITY, - data.data.LEGATO_PRESSURE, - coefs) - - data.data.loc[~(abs(data.data['salinity']-data.data['LEGATO_SALINITY']) < 1),'salinity'] = np.NaN # Can definitely play around with this threshold - return data - -def lagCorrection(variable,referencevariable,time,lag=None): - try: - time = date2float(time) - print('Applying lag after time conversion:') - except: - print('Applying lag:') - - def _Correction(x): - return variable + (x * np.gradient(variable,time)) - - if lag is None: - from scipy import optimize - - -# def _PolyArea_triangles(sal): -# S,X,Y = grid2d(data.data.profileNum.values[_dives],data.data.LEGATO_PRESSURE.values[_dives],sal,xi=1,yi=0.3,fn='mean') -# def merge_profiles(idx): -# x1 = S[:,idx] -# x2 = S[:,idx+1] -# y1 = T[:,idx] -# y2 = T[:,idx+1] -# _gd = np.isfinite(x1+x2+y1+y2) -# x1 = x1[_gd] -# x2 = x2[_gd] -# y1 = y1[_gd] -# y2 = y2[_gd] -# return x1,x2,y1,y2 -# def triangleArea(x1,y1,x2,y2,x3,y3): -# return abs((0.5)*(x1*(y2-y3)+x2*(y3-y1)+x3*(y1-y2))) -# area = 0 -# for idd in range(np.shape(T)[1]-1): -# x1,x2,y1,y2 = merge_profiles(idd) -# area = area + np.nansum(triangleArea( -# np.append(x1,np.NaN),np.append(y1,np.NaN), -# np.append(x2,np.NaN),np.append(y2,np.NaN), -# np.append(np.NaN,x1),np.append(np.NaN,y1) )) -# area = area + np.nansum(triangleArea( -# np.append(np.NaN,x1),np.append(np.NaN,y1), -# np.append(x2,np.NaN),np.append(y2,np.NaN), -# np.append(np.NaN,x2),np.append(np.NaN,y2) )) -# return area - - def _PolyArea(x,y): - _gg = np.isfinite(x+y) - return 0.5*np.abs(np.dot(x[_gg],np.roll(y[_gg],1))-np.dot(y[_gg],np.roll(x[_gg],1))) - # Doesn't work well with intersecting polygon.... Problem?? - - - def _f(x): - return _PolyArea(_Correction(x),referencevariable) - print('Regressing lag coefficient:') - minimum = optimize.fmin(_f, 0.5, disp=True, full_output=True) - print('Calculated lag coefficient = '+np.array2string(minimum[0])) - lag = minimum[0] - - return _Correction(lag) - - -########################################################################################################################################################## -########################################################################################################################################################## -########################################################################################################################################################## -# # -# SEMI-DYNAMIC FLIGHT MODEL CLASS WITH DUAL CONSTRAINTS # -# # -########################################################################################################################################################## -########################################################################################################################################################## -########################################################################################################################################################## - -#New Steady State model (With adcp constraint) -class SemiDynamicModel(object): - def __init__(self, time, sal, temp_ext, pres_ext, lon, lat, ballast, pitch, profile, navresource, horz_vert_weighting, speed_through_water, **param): - # Questions: - # is dzdt spiky? - # do values need to be interpolated? - - # Parse input variables: - self.timestamp = time - self.time = date2float(self.timestamp) - - self.external_pressure = pres_ext - self.external_temperature = temp_ext - - self.longitude = lon - self.latitude = lat - self.profile = profile - - self.salinity = sal - - self.ballast = ballast/1000000 # m^3 - self.pitch = np.deg2rad(pitch) # rad - - self.navresource=navresource - - self.hv_weighting=horz_vert_weighting - self.speed_through_water=speed_through_water - - - # Rerefence parameters for scaling and initial guesses: - self.param_reference = dict({ - 'mass': 60, # Vehicle mass in kg - 'vol0': 0.06, # Reference volume in m**3, with ballast at 0 (with -500 to 500 range), at surface pressure and 20 degrees C - 'area_w': 0.24, # Wing surface area, m**2 - - 'Cd_0': 0.046, # - 'Cd_1': 2.3, # - 'Cl': 2.0, # Negative because wrong convention on theta - - 'comp_p': 4.7e-06, #1.023279317627415e-06, #Pressure dependent hull compression factor - 'comp_t': 8.9e-05, #1.5665248101730484e-04, # Temperature dependent hull compression factor - - 'SS_tau': 13, #characteristic response time of the glider in sec - - }) - - self.param = self.param_reference.copy() - for k,v in param.items(): - self.param[k] = v - self.param_initial = self.param - - self.depth = gsw.z_from_p(self.external_pressure,self.latitude) # m . Note depth (Z) is negative, so diving is negative dZdt - self.dZdt = np.gradient(self.depth,self.time) # m.s-1 - - self.g = gsw.grav(self.latitude,self.external_pressure) - - self.SA = gsw.SA_from_SP(self.salinity, self.external_pressure, self.longitude, self.latitude) - self.CT = gsw.CT_from_t(self.SA, self.external_temperature, self.external_pressure) - self.rho = gsw.rho(self.SA, self.CT, self.external_pressure) - - ### Basic model - # Relies on steady state assumption that buoyancy, weight, drag and lift cancel out when not accelerating. - # F_B - cos(glide_angle)*F_L - sin(glide_angle)*F_D - F_g = 0 - # cos(glide_angle)*F_L + sin(glide_angle)*F_D = 0 - - # Begin with an initial computation of angle of attack and speed through water: - self.model_function() - - ### Get good datapoints to regress over - self._valid = np.full(np.shape(self.pitch),True) - self._valid[self.external_pressure < 3] = False - self._valid[np.abs(np.rad2deg(self.pitch)) < 10] = False - self._valid[np.abs(np.rad2deg(self.pitch)) > 60] = False - self._valid[np.abs(np.gradient(self.dZdt,self.time)) > 0.005] = False # Accelerations - self._valid[np.gradient(self.pitch,self.time)==0] = False # Rotation - # self._valid = self._valid & ((navresource == 100) | (navresource == 117) ) #100=glider going down & 117=glider going up (=> not at surface or inflecting) - - # Do first pass regression on vol parameters, or volume and hydro? - self.regression_parameters = ('vol0','Cd_0','Cd_1','Cl','comp_p','comp_t') - # ('vol0','comp_p','comp_t','Cd_0','Cd_1','Cl', 'SS_tau') # Has to be a tuple, requires a trailing comma if single element - - print('Number of valid points: '+str(np.count_nonzero(self._valid))+' (out of '+str(len(self._valid))+')') - - - ### Principal forces - @property - def F_B(self): - return self.g * self.rho * (self.ballast + self.vol0 * - (1 - self.comp_p*(self.external_pressure) + self.comp_t*(self.external_temperature))) - @property - def F_g(self): - return self.mass * self.g - - @property - def vert_dir(self): - return np.sign(self.F_B-self.F_g) #Positive is buoyancy force up & negative is buoyancy force down - -# @property -# def F_L(self): -# return self.dynamic_pressure * self.area_w * (self.Cl) * self.alpha - -# @property -# def F_D(self): -# return self.dynamic_pressure * self.area_w * (self.Cd_0 + self.Cd_1*(self.alpha)**2) - - ### Important variables - @property - def glide_angle(self): - return self.pitch + self.alpha - - @property - def dynamic_pressure(self): - return self.rho * self.speed**2 / 2 - - @property - def w_H2O(self): - # Water column upwelling - return self.dZdt - self.speed_vert - - @property - def R1(self): - return np.sqrt(np.nanmean( (1-self.hv_weighting)*self.w_H2O[self._valid]**2 + self.hv_weighting*((self.speed_through_water - self.speed)[self._valid]**2) )) - - ### Basic equations - def _solve_alpha(self): - _pitch_range = np.linspace( np.deg2rad(0), np.deg2rad(90) , 90) - _alpha_range1 = np.zeros_like(_pitch_range) - _alpha_range2 = np.zeros_like(_pitch_range) - - # Resolve for normal flight - for _istep, _pitch in enumerate(_pitch_range): - _tmp = fsolve(self._equation_alpha, 0.01, args=(_pitch), full_output=True) - _alpha_range1[_istep] = _tmp[0] - - # Resolve for stall ## VERIFY WHAT THEY DID HERE - for _istep, _pitch in enumerate(_pitch_range): - if (np.sign(_pitch)>0) : - _tmp = fsolve(self._equation_alpha, (-np.pi/2 -_pitch + 0.01), args=(_pitch), full_output=True) - _alpha_range2[_istep] = _tmp[0] - else : - _tmp = fsolve(self._equation_alpha, (np.pi/2 -_pitch - 0.01), args=(_pitch), full_output=True) - _alpha_range2[_istep] = _tmp[0] - - _interp_fn1 = interp1d(_pitch_range,_alpha_range1) - _interp_fn2 = interp1d(_pitch_range,_alpha_range2) - - Res=_interp_fn1(np.abs(self.pitch)) * np.sign(self.pitch) #Résolution noramle - Res[self.vert_dir*np.sign(self.pitch)<0]=(_interp_fn2(np.abs(self.pitch)) * np.sign(self.pitch))[self.vert_dir*np.sign(self.pitch)<0] #Résolution décrochage - return Res - - def _equation_alpha(self, _alpha, _pitch): - return ( self.Cd_0 + self.Cd_1 * (_alpha**2) ) / ( (self.Cl) * np.tan(_alpha + _pitch) ) - _alpha - - def _solve_speed(self): - _dynamic_pressure = (self.F_B - self.F_g) * np.sin(self.glide_angle) / ( self.Cd_0 + self.Cd_1 * (self.alpha)**2 ) - return np.sqrt(2 * _dynamic_pressure / self.rho / self.area_w) - - def model_function(self): - - def fillGaps(x,y): - fill = lambda arr: pd.DataFrame(arr).fillna(method='bfill').fillna(method='ffill').values.flatten() - f = interp1d(x[np.isfinite(x+y)],y[np.isfinite(x+y)], bounds_error=False, fill_value=np.NaN) - return(fill(f(x))) - - sos = signal.butter(1, 1/(2*np.pi*self.SS_tau), 'lowpass', fs=1, output='sos') #order, cutoff frequency, "lowpass", sampling frequency, - _lowpassDynamics = lambda arr: signal.sosfilt(sos, arr) - - - self.alpha = fillGaps(self.time, self._solve_alpha()) # - - self.raw_speed = self._solve_speed() # fillGaps(self.time, self._solve_speed()) - self.raw_speed [ (self.external_pressure < 1) & (self.navresource == 116) ] = 0 #Set speed to be zero at the surface ## TODO: Not a good way to do - - self.speed = _lowpassDynamics(fillGaps(self.time, self.raw_speed)) - self.speed_vert = np.sin(self.glide_angle)*self.speed - self.speed_horz = np.cos(self.glide_angle)*self.speed - - def cost_function(self,x_initial): - for _istep, _key in enumerate(self.regression_parameters): - self.param[_key] = x_initial[_istep] * self.param_reference[_key] - self.model_function() - return self.R1 - - def regress(self, maxiter = 300): - x_initial = [self.param[_key] / self.param_reference[_key] for _istep,_key in enumerate(self.regression_parameters)] - print('Initial parameters: ', self.param) - print('Non-optimised score: '+str(self.cost_function(x_initial)) ) - print('Regressing...') - - with tqdm(total=maxiter) as pbar: - def callbackF(Xi): - pbar.update(1) - R = fmin(self.cost_function, x_initial, callback=callbackF, disp=True, full_output=True, maxiter=maxiter, ftol=0.00001) - - for _istep,_key in enumerate(self.regression_parameters): - self.param[_key] = R[0][_istep] * self.param_reference[_key] - - x_initial = [self.param[_key] / self.param_reference[_key] for _istep,_key in enumerate(self.regression_parameters)] - print('Optimised parameters: ', self.param) - print('Final Optimised score: '+str(self.cost_function(x_initial)) ) - self.model_function() - - ### Coefficients - @property - def mass(self): - return self.param['mass'] - - @property - def vol0(self): - return self.param['vol0'] - - @property - def comp_p(self): - return self.param['comp_p'] - - @property - def comp_t(self): - return self.param['comp_t'] - - @property - def area_w(self): - return self.param['area_w'] - - @property - def Cd_0(self): - return self.param['Cd_0'] - - @property - def Cd_1(self): - return self.param['Cd_1'] - - @property - def Cl(self): - return self.param['Cl'] - - @property - def SS_tau(self): - return self.param['SS_tau'] \ No newline at end of file