#! /usr/bin/env python # # =========================================================== # # Alexander R Davies # # Ocean Exploration, Remote Sensing, and Biogeography Lab # School of Marine Science and Policy # College of Earth, Ocean, and Environment # University of Delaware # ardavies@udel.edu # # Latest Update: 01/13/2014 # # NOTES: 1) Does all the OSCAR Data Processing # # # # =========================================================== # Import Modules # =========================================================== import numpy as np import matplotlib matplotlib.use('Agg') import matplotlib.pyplot as plt from datetime import datetime from mpl_toolkits.basemap import Basemap, addcyclic from mpl_toolkits.basemap import shiftgrid, cm from scipy.ndimage.filters import minimum_filter, maximum_filter from netCDF4 import Dataset import glob import os from scipy.io import netcdf import matplotlib.transforms as mtransforms from matplotlib.patches import FancyBboxPatch # from matplotlib.transforms import Bboxcx from matplotlib.path import Path import matplotlib.patches as patches import math as ma import csv from matplotlib.ticker import MaxNLocator from matplotlib import rc, rcParams from numpy.random import uniform, seed from matplotlib.mlab import griddata from numpy import arange,array,ones#,random,linalg from pylab import plot,show from scipy import interpolate from matplotlib.colors import LogNorm from matplotlib.backends.backend_pdf import PdfPages # # =========================================================== # # FUNCTIONS # # =========================================================== # # =========================================================== # Function to input delX, delY, return bearing and magnitude # =========================================================== # def find_bearing_magnitude(X, Y): import math as m if ((Y >= 0) and (X >=0 )): bearing_rads = m.atan(abs(X)/abs(Y)) bearing = m.degrees(bearing_rads) magnitude = abs(X)/m.sin(bearing_rads) elif ((Y < 0) and (X >=0)): bearing_rads = m.atan(abs(Y)/abs(X)) bearing = m.degrees(bearing_rads) + 90 magnitude = abs(Y)/m.sin(bearing_rads) elif ((Y < 0) and (X < 0)): bearing_rads = m.atan(abs(X)/abs(Y)) bearing = m.degrees(bearing_rads) + 180 magnitude = abs(X)/m.sin(bearing_rads) elif ((Y >= 0) and (X < 0)): bearing_rads = m.atan(abs(Y)/abs(X)) bearing = m.degrees(bearing_rads) + 270 magnitude = abs(Y)/m.sin(bearing_rads) else: print "#### ERROR IN: FIND BEARING AND MAG FUNCTION ####" return bearing_rads, bearing, magnitude # # =========================================================== # Function Determine new Lat/Lon from Old + Bearing + dist # =========================================================== # def destpoint(lat1, lon1, brng, d): import math R = 6374 #Radius of the Earth # brng = Bearing converted to radians. # d = Distance in km # lat1 = Current lat point converted to radians # lon1 = Current lon point converted to radians lat2 = math.asin(math.sin(lat1)*math.cos(d/R) + math.cos(lat1)*math.sin(d/R)*math.cos(brng)) lon2 = lon1 + math.atan2(math.sin(brng)*math.sin(d/R)*math.cos(lat1), math.cos(d/R)-math.sin(lat1)*math.sin(lat2)) lat22 = math.degrees(lat2) lon22 = math.degrees(lon2) return lat22, lon22 # # =========================================================== # Function to read csv files # =========================================================== # def csvread(filename): with open(filename, 'rb') as f: reader = csv.reader(f, delimiter=',') nrow = 0; for row in reader: nrow = nrow +1 ncol = len(row) array = np.zeros([ncol,nrow]) with open(filename, 'rb') as f: reader = csv.reader(f, delimiter=',') counter = 0 for row in reader: for jj in range(0,ncol): array[jj,counter] = row[jj] counter = counter + 1 return array, nrow, ncol # # =========================================================== # Function to Find the nearest value # =========================================================== # def find_nearest(array,value): idx = (np.abs(array-value)).argmin() return array[idx], idx # # =========================================================== # Function to Calculate Distances on a Sphere in kilometers # =========================================================== # # def pos2dist(pos): # # Pos = [lon1,lat1,lon2,lat2] # import math as m # R_aver = 6374 # pi = 3.14 # deg2rad = pi/180 # lag3 = pos[1] * deg2rad # lon3 = pos[0] * deg2rad # lag4 = pos[3] * deg2rad # lon4 = pos[2] * deg2rad # dist = R_aver * m.acos(m.cos(lag3)*m.cos(lag4)*m.cos(lon3-lon4) + m.sin(lag3)*m.sin(lag4)) # # Returns a distance in KM # return dist def pos2dist(pos): # Pos = [lon1,lat1,lon2,lat2] import math as m R_aver = 6374 pi = 3.14 deg2rad = pi/180 lat1 = pos[1] * deg2rad #lat1 lon1 = pos[0] * deg2rad # lon1 lat2 = pos[3] * deg2rad #lat2 lon2 = pos[2] * deg2rad #lon2 dlon = lon2 - lon1 dlat = lat2 - lat1 e = m.cos(lat1)*m.cos(lat2)*m.sin(dlon/2)**2 d = m.sin(dlat/2)**2 c = d + e b = c**0.5 dist = 2*R_aver*m.asin(b) # Returns a distance in KM return dist # # =========================================================== # Function to new lon from zonal displacemet in km # =========================================================== # def lonchange(input): # input = [lat,lon1,dzon] import math as m R_aver = 6374 pi = 3.14 deg2rad = pi/180 rad2deg = 180/pi lat = input[0] * deg2rad #lat1 lon1 = input[1] * deg2rad # lon1 dzon = input[2] # zonal displacement in km2 c = m.cos(lat) b = m.sin(dzon/2/R_aver) a = b/c lon2rad = 2*m.asin(a) + lon1 lon2 = lon2rad*rad2deg return lon2 # # =========================================================== # Function to new lat from meridional displacemet in km # =========================================================== # def latchange(input): # input = [lat1,dmer] import math as m R_aver = 6374 pi = 3.14 deg2rad = pi/180 rad2deg = 180/pi lat1 = input[0]*deg2rad dmer = input[1] # zonal displacement in km2 lat2rad = dmer/R_aver + lat1 lat2 = lat2rad*rad2deg return lat2 # # =========================================================== # Function to Do Regressions # =========================================================== # def linearregress(x, y): from scipy import stats # # Linearly Regression slope, intercept, rvalue, pvalue, stderr = stats.linregress(x,y) # # Regression Line regressline = slope*x+intercept # return regressline, rvalue, pvalue, stderr # # # =================================================================================== # Function for Bilinear Interpolation # =================================================================================== # def bilinear_interpolation(x, y, points): # # Order Points vy x, then by y points = sorted(points) # order points by x, then by y (x1, y1, q11), (_x1, y2, q12), (x2, _y1, q21), (_x2, _y2, q22) = points # # Are the points a rectangle? if x1 != _x1 or x2 != _x2 or y1 != _y1 or y2 != _y2: raise ValueError('points do not form a rectangle') if not x1 <= x <= x2 or not y1 <= y <= y2: raise ValueError('(x, y) not within the rectangle') # # Bilinear Interp and return value. return (q11 * (x2 - x) * (y2 - y) + q21 * (x - x1) * (y2 - y) + q12 * (x2 - x) * (y - y1) + q22 * (x - x1) * (y - y1) ) / ((x2 - x1) * (y2 - y1) + 0.0) # # # =================================================================================== # Function to Calculate the Backscatter of Seawater # =================================================================================== # def betasw_ZHH2009(Lambda,Tc,theta,S): # # Xiaodong Zhang, Lianbo Hu, and Ming-Xia He (2009), Scatteirng by pure seawater: Effect of salinity, Optics Express, Vol. 17, No. 7, 5698-5710 # # Lambda (nm): wavelength # Tc: temperauter in degree Celsius, must be a scalar # S: salinity, must be scalar # 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 Lambda (y = length(Lambda)) # 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 # # Originally By: Xiaodong Zhang, March 10, 2009 (Modified for Python by: Alexander R. Davies, April 13, 2014) # import numpy as np # # =================================================================================== # Values and Constants # =================================================================================== # #Avogadro's constant Na = 6.0221417930e23 # # Boltzmann constant Kbz = 1.3806503e-23 # # Absolute tempearture Tk = Tc+273.15 # # Molecular weigth of water in kg/mol M0 = 18e-3 # # Depolarization ratio, default = 0.039 will be used from Farinato and Roswell (1976). You can change this!! delta = 0.039 # # Degrees to Radians rad = np.radians(theta) # # =================================================================================== # absolute refractive index of seawater (nsw); d(nsw) w.r.t. salinity # =================================================================================== # # refractive index of air is from Ciddor (1996,Applied Optics) #n_air = 1.0 + (5792105.0./(238.0185 -1./(Lambda/1e3).^2)+167917.0./(57.362-1./(Lambda/1e3 ).^2))/1e8; << Original Matlab Code n_air = 1.0 + (5792105.0/(238.0185 -1 /(Lambda/1e3)**2)+167917.0 /(57.362-1 /(Lambda/1e3)**2))/10**8 # # 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 # # For Pure Water First # nsw = n0+(n1+n2*Tc+n3*Tc^2 )*S+n4*Tc^2 +(n5+n6*S+n7*Tc)./Lambda + n8./Lambda.^2 + n9./Lambda.^3; << Original Matlab Code nsw = n0+(n1+n2*Tc+n3*Tc**2)*S+n4*Tc**2+(n5+n6*S+n7*Tc) /Lambda + n8 /Lambda**2 + n9 /Lambda**3 # # Saltwater nsw = nsw*n_air # # Derivative #dnswds = (n1+n2*Tc+n3*Tc^2 +n6./Lambda).*n_air; << Original Matlab Code dnswds = (n1+n2*Tc+n3*Tc**2+n6/Lambda)*n_air dnds = dnswds # # # =================================================================================== # isothermal compressibility # =================================================================================== # # From Lepple & Millero (1971,Deep Sea-Research), pages 10-11; error ~ +/-0.004e-6 bar^-1 # # 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; << Original Matlab Code kw = 19652.21 + 148.4206*Tc - 2.327105*Tc**2 + 1.360477e-2*Tc**3 - 5.155288e-5*Tc**4; Btw_cal = 1/kw # # # seawater secant bulk #a0 = 54.6746 - 0.603459*Tc +1.09987e-2*Tc.^2 - 6.167e-5*Tc.^3; << Original Matlab Code 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; << Original Matlab Code b0 = 7.944e-2 + 1.6483e-2*Tc -5.3009e-4*Tc**2 # #Ks = kw + a0*S + b0*S.^1.5; << Original Matlab Code Ks = kw + a0*S + b0*S**1.5 # # calculate seawater isothermal compressibility from the secant bulk IsoComp = 1./Ks*1e-5; # unit is pa # # =================================================================================== # Density of Water # =================================================================================== # # 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; << Original Matlab Code 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); << Original Matlab Code 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); # # =================================================================================== # Partial derivative of natural logarithm of water activity w.r.t.salinity # =================================================================================== # # 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; << Original Matlab Code 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 # # =================================================================================== # Density derivative of refractive index from PMH model # =================================================================================== # n_wat = nsw 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); << Original Matlab Code c1 = (n_wat2-1) c21 = 2./3.*(n_wat2+2.) c22 = n_wat/3. c23 = 1./3./n_wat c2 = (1. + c21*(c22-c23)**2.) n_density_derivative = c1*c2 DFRI = n_density_derivative # # =================================================================================== # Volume Scattering # =================================================================================== # # volume scattering at 90 degree due to the density fluctuation #beta_df = pi*pi/2*((Lambda*1e-9).^(-4))*Kbz*Tk*IsoComp.*DFRI.^2*(6+6*delta)/(6-7*delta); << Original Matlab Code beta_df = np.pi*np.pi/2*((Lambda*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; << Matlab Code flu_con = S*M0*dnds**2/density_sw/(-dlnawds)/Na # #beta_cf = 2*pi*pi*((Lambda*1e-9).^(-4)).*nsw.^2.*(flu_con)*(6+6*delta)/(6-7*delta); << Original Matlab Code beta_cf = 2*np.pi*np.pi*((Lambda*1e-9)**(-4))*nsw**2*(flu_con)*(6 + 6*delta)/(6 -7*delta) # # total volume scattering at 90 degree beta90sw = beta_df+beta_cf # # volume scattering coefficient bsw=8*np.pi/3*beta90sw*(2 + delta)/(1 + delta) # # Volume scattering at the angle defined by theta betasw=beta90sw*(1+((np.cos(rad))**2)*(1-delta)/(1+delta)) return betasw,beta90sw,bsw # # =========================================================== # # READING FLOAT DATA # # =========================================================== # # # =================================================================================== # Initial Set-up to read float data # =================================================================================== # # Get to the float data... os.chdir('/data/orbprocess_mail/alex/Jan01_Jun04_2013/data/type/noday71') # # Compile a list of the filenames and sort them based on the time they were generated fullnames = glob.glob('Full*') gradnames = glob.glob('Grad*') infonames = glob.glob('Text*') fullnames.sort(key = lambda x: os.path.getmtime(x)) gradnames.sort(key = lambda x: os.path.getmtime(x)) infonames.sort(key = lambda x: os.path.getmtime(x)) # # Array Initialization arraylen = len(infonames) floatdate = np.zeros(arraylen) floatlat = [None]*arraylen floatlon = [None]*arraylen day = np.zeros(arraylen) month=[None]*arraylen year=[None]*arraylen daystr=[None]*arraylen listdate =[None]*arraylen # # =========================================================== # Processing Float Time and Location Data # =========================================================== # # Retrive the Information for Each Data Profile for d in range(0,int(arraylen)): # # Get the Date Info f = open('/data/orbprocess_mail/alex/Jan01_Jun04_2013/data/type/noday71/' + infonames[d], 'r') while 1: line = f.readline() if line[:4] == "Date": month[d] = str(line[5:7]) day[d] = int(line[8:10]) daystr[d] = str(line[8:10]) year[d] = str(line[11:15]) break # # Get the LaT and Lon f = open('/data/orbprocess_mail/alex/Jan01_Jun04_2013/data/type/noday71/' + infonames[d], 'r') while 1: line = f.readline() if line[:3] == "Lat": floatlat[d] = str(line[5:12]) break f = open('/data/orbprocess_mail/alex/Jan01_Jun04_2013/data/type/noday71/' + infonames[d], 'r') while 1: line = f.readline() if line[:3] == "Lon": floatlon[d] = str(line[5:12]) break # # Get the DOY for each profile if (year[d] == '2013'): yr = 0 elif (year[d] == '2014'): yr = 365 elif (year[d] == '2015'): yr = 365*2 if(month[d] == '01'): mon = 0 elif(month[d] == '02'): mon = 31 elif(month[d] == '03'): mon = 31 + 28 elif(month[d] == '04'): mon = (31 + 28 + 31) elif(month[d] == '05'): mon = (31 + 28 + 31 + 30) elif(month[d] == '06'): mon = (31 + 28 + 31 + 30 + 31) elif(month[d] == '07'): mon = (31 + 28 + 31 + 30 + 31 + 30) elif(month[d] == '08'): mon = (31 + 28 + 31 + 30 + 31 + 30 + 31) elif(month[d] == '09'): mon = (31 + 28 + 31 + 30 + 31 + 30 + 31 + 31) elif(month[d] == '10'): mon = (31 + 28 + 31 + 30 + 31 + 30 + 31 + 31 + 30) elif(month[d] == '11'): mon = (31 + 28 + 31 + 30 + 31 + 30 + 31 + 31 + 30 + 31) elif(month[d] == '12'): mon = (31 + 28 + 31 + 30 + 31 + 30 + 31 + 31 + 30 + 31 + 30) else: print 'invalid month' floatdate[d] = yr + mon + day[d] # # Change the longitude format and make lat/lon floats for i in range(0,arraylen): floatlon[i] = 360 + float(floatlon[i]) floatlat[i] = float(floatlat[i]) # # Get Back to start directory os.chdir('/home/ardavies/satdata/OSCAR') # # =========================================================== # # READING OSCAR DATA # # =========================================================== # # ncdump file header to .txt file os.system('ncdump -h oscar-third-Jan5-jun20.nc > oscar-third-Jan5-jun20.txt') # # open the netcdf file to read ncf = netcdf.netcdf_file('oscar-third-Jan5-jun20.nc','r') # # Getting the netcdf data into the corresponding arrays datadates = ncf.variables['time'][:].squeeze() v = ncf.variables['v'][:].squeeze() u = ncf.variables['u'][:].squeeze() rawlon = ncf.variables['lon'][:].squeeze() # degrees E rawlat = ncf.variables['lat'][:].squeeze() # degrees N # # Array Initailization numdates = len(u[:,0,0]) lonlen = len(rawlon) latlen = len(rawlat) lat = np.zeros([latlen,lonlen]) lon = np.zeros([latlen,lonlen]) uavg = np.zeros([numdates-1,latlen,lonlen]) vavg = np.zeros([numdates-1,latlen,lonlen]) # # Getting lat/lon's in i,j arrays for i in range(0,lonlen): lon[:,i] = rawlon[i] for j in range(0,latlen): lat[j,:] = rawlat[j] # # Average Velocities for k in range(0,numdates-1): for i in range(0,lonlen): for j in range(0,latlen): uavg[k,j,i] = (u[k,j,i] + u[k+1,j,i])/2 vavg[k,j,i] = (v[k,j,i] + v[k+1,j,i])/2 # # DOY of 5 day avg Oscar data centerdates = np.zeros(numdates) for k in range(0,numdates): centerdates[k] = 6 + datadates[k] # # =========================================================== # # DAILY AVERAGED OSCAR DATA # # =========================================================== # # Daily Averaged OSCAR Data Array Initialization oscardates = np.linspace(centerdates[0],centerdates[numdates-1],(centerdates[numdates-1]-centerdates[0])+1) # udaily = np.zeros([len(oscardates)-8,latlen,lonlen]) vdaily = np.zeros([len(oscardates)-8,latlen,lonlen]) # for ii in range(0,len(oscardates)-8): # # Get the closest OSCAR data date to the current day in the loop. k is index of the closest centerdate closestdate, k = find_nearest(centerdates, oscardates[ii]) # # If the closest data day and the current day in the loop are the same, the current day take 100% of the u and v velocities from the oscar data day if (closestdate == oscardates[ii]): udaily[ii,:,:] = u[k,:,:] vdaily[ii,:,:] = v[k,:,:] # # What if k = 0? We must "weigh" the u and v currents on the current day of the loop between centerdates[k=0] and centerdates[(k=0) +1] as k-1 doesn't exist elif (k == 0): # # next closest date must be at k+1 nextclosestdate = centerdates[k+1] kk = k + 1 # # The u and v components at the current day in the loop are weighed by how close they are to the oscar data dates. udaily[ii,:,:] = u[k,:,:]*(abs(oscardates[ii] - nextclosestdate)/abs(nextclosestdate-closestdate)) + u[kk,:,:]*(abs(oscardates[ii] - closestdate)/abs(nextclosestdate-closestdate)) vdaily[ii,:,:] = v[k,:,:]*(abs(oscardates[ii] - nextclosestdate)/abs(nextclosestdate-closestdate)) + v[kk,:,:]*(abs(oscardates[ii] - closestdate)/abs(nextclosestdate-closestdate)) # # If the dates do not match exactly, and k doesn't equal 0 then we need to figure out if the next closest oscar data dates are before or after centerdates[k] else: # # Next closest day to the current day of the loop is AFTER centerdates[k] if ((abs(oscardates[ii]-centerdates[k+1]))<=(abs(oscardates[ii]-centerdates[k-1]))): nextclosestdate = centerdates[k+1] kk = k+1 # # Next closest day to the current day of the loop is BEFORE centerdates[k] else: nextclosestdate = centerdates[k-1] kk = k-1 # # The u and v components at the current day in the loop are weighed by how close they are to the oscar data dates. udaily[ii,:,:] = u[k,:,:]*(abs(oscardates[ii] - nextclosestdate)/abs(nextclosestdate-closestdate)) + u[kk,:,:]*(abs(oscardates[ii] - closestdate)/abs(nextclosestdate-closestdate)) vdaily[ii,:,:] = v[k,:,:]*(abs(oscardates[ii] - nextclosestdate)/abs(nextclosestdate-closestdate)) + v[kk,:,:]*(abs(oscardates[ii] - closestdate)/abs(nextclosestdate-closestdate)) # # =========================================================== # # MKE # # =========================================================== # import math as ma mke = np.zeros([len(oscardates)-8,latlen,lonlen]) # # MKE for Daily Weighed OSCAR Currents for k in range(0,len(oscardates)-8): for i in range(0,lonlen): for j in range(0,latlen): mke[k,j,i] = 0.5*((udaily[k,j,i]*100)**2 + (vdaily[k,j,i]*100)**2) # # =========================================================== # # OSCAR CURRENT AND MAKE INTERPOLATION TO FLOAT LOCATION # # =========================================================== # # # =================================================================================== # Find Nearest Lat's and Lon's to float location for U and V components # =================================================================================== # latvalue1 = np.zeros([arraylen]) latvalue2 = np.zeros([arraylen]) lonvalue1 = np.zeros([arraylen]) lonvalue2 = np.zeros([arraylen]) # latindex1 = np.zeros([arraylen]) latindex2 = np.zeros([arraylen]) lonindex1 = np.zeros([arraylen]) lonindex2 = np.zeros([arraylen]) # for jj in range(0,arraylen): # # Get the closest OSCAR data latitude with respect to the float latitude latvalue1[jj], latindex1[jj] = find_nearest(lat[:,0], floatlat[jj]) # # Find the next closest OSCAR data latitude with respect to the float latitude upperlat = lat[latindex1[jj]+1,0] lowerlat = lat[latindex1[jj]-1,0] # # Closest lat + 0.33 degrees? if (abs(floatlat[jj]-upperlat)<=abs(floatlat[jj]-lowerlat)): latindex2[jj] = latindex1[jj] + 1 # # Closest lat - 0.33 degrees? else: latindex2[jj] = latindex1[jj] - 1 latvalue2[jj] = lat[latindex2[jj],0] # # Get the closest OSCAR data longitude with respect to the float longitude lonvalue1[jj], lonindex1[jj] = find_nearest(lon[0,:], floatlon[jj]) # # Find the next closest OSCAR data longitude with respect to the float longitude upperlon = lon[0,lonindex1[jj]+1] lowerlon = lon[0,lonindex1[jj]-1] # # Closest lon + 0.33 degrees? if (abs(floatlon[jj]-upperlon)<=abs(floatlon[jj]-lowerlon)): lonindex2[jj] = lonindex1[jj] + 1 # # Closest lon - 0.33 degrees? else: lonindex2[jj] = lonindex1[jj] - 1 lonvalue2[jj] = lon[0,lonindex2[jj]] # # =================================================================================== # Daily Bi-Linear Interpolation # =================================================================================== # # Initializing Arrays for Interpolation floatU = np.zeros([arraylen]) floatV = np.zeros([arraylen]) floatMKE = np.zeros([arraylen]) # # Loop through the float data for jj in range(0,arraylen): # # Loop though the daily oscar data for k in range(0,len(oscardates)-8): # # Loop through the daily oscar data ates to find where it equals the float data dates if (oscardates[k] == floatdate[jj]): # # Interpolate Linearly Along Closest Latitude (latvalue1) x1_closest = abs(pos2dist([lonvalue1[jj],latvalue1[jj],floatlon[jj],latvalue1[jj]])) x2_closest = abs(pos2dist([lonvalue2[jj],latvalue1[jj],floatlon[jj],latvalue1[jj]])) xtot_closest = abs(pos2dist([lonvalue1[jj],latvalue1[jj],lonvalue2[jj],latvalue1[jj]])) Uinterp_closest = udaily[k,latindex1[jj], lonindex1[jj]]*(x2_closest/xtot_closest) + udaily[k, latindex1[jj], lonindex2[jj]]*(x1_closest/xtot_closest) Vinterp_closest = vdaily[k,latindex1[jj], lonindex1[jj]]*(x2_closest/xtot_closest) + vdaily[k, latindex1[jj], lonindex2[jj]]*(x1_closest/xtot_closest) MKEinterp_closest = mke[k,latindex1[jj], lonindex1[jj]]*(x2_closest/xtot_closest) + mke[k, latindex1[jj], lonindex2[jj]]*(x1_closest/xtot_closest) # # Interpolate Linearly Along Furthest Latitude (latvalue2) x1_furthest = abs(pos2dist([lonvalue1[jj],latvalue2[jj],floatlon[jj],latvalue2[jj]])) x2_furthest = abs(pos2dist([lonvalue2[jj],latvalue2[jj],floatlon[jj],latvalue2[jj]])) xtot_furthest = abs(pos2dist([lonvalue1[jj],latvalue2[jj],lonvalue2[jj],latvalue2[jj]])) Uinterp_furthest = udaily[k,latindex2[jj], lonindex1[jj]]*(x2_furthest/xtot_furthest) + udaily[k, latindex2[jj], lonindex2[jj]]*(x1_furthest/xtot_furthest) Vinterp_furthest = vdaily[k,latindex2[jj], lonindex1[jj]]*(x2_furthest/xtot_furthest) + vdaily[k, latindex2[jj], lonindex2[jj]]*(x1_furthest/xtot_furthest) MKEinterp_furthest = mke[k,latindex2[jj], lonindex1[jj]]*(x2_furthest/xtot_furthest) + mke[k, latindex2[jj], lonindex2[jj]]*(x1_furthest/xtot_furthest) # # Tnterpolate Along Londitude y_closest = abs(pos2dist([floatlon[jj],latvalue1[jj],floatlon[jj],floatlat[jj]])) y_furthest = abs(pos2dist([floatlon[jj],latvalue2[jj],floatlon[jj],floatlat[jj]])) ytot = abs(pos2dist([floatlon[jj],latvalue1[jj],floatlon[jj],latvalue2[jj]])) floatU[jj]= Uinterp_furthest*(y_closest/ytot) + Uinterp_closest*(y_furthest/ytot) floatV[jj]= Vinterp_furthest*(y_closest/ytot) + Vinterp_closest*(y_furthest/ytot) floatMKE[jj]= MKEinterp_furthest*(y_closest/ytot) + MKEinterp_closest*(y_furthest/ytot) # # =========================================================== # # FLOAT & OCSAR CURRENT TRACERS # # =========================================================== # # # =================================================================================== # The Tracer U and V components for each float surfacing # =================================================================================== # # Initialize Arrays tracerU_kmday = np.zeros([arraylen]) tracerV_kmday = np.zeros([arraylen]) tracerX_km = np.zeros([arraylen]) tracerY_km = np.zeros([arraylen]) for jj in range(0,arraylen): # # U and V Velocitues in km/day tracerU_kmday[jj] = floatU[jj]*86400/1000 tracerV_kmday[jj] = floatV[jj]*86400/1000 # # Distance Tracer will go in 2 days tracerX_km[jj] = tracerU_kmday[jj]*2 tracerY_km[jj] = tracerV_kmday[jj]*2 # # =========================================================== # Use haversine's to calculate new lat and lon of tracer # =========================================================== # tracer_newlat_U_V = np.zeros([arraylen]) tracer_newlon_U_V = np.zeros([arraylen]) magnitude_U_V = np.zeros([arraylen]) for jj in range(0,arraylen): tracer_newlat_U_V[jj] = latchange([floatlat[jj],tracerY_km[jj]]) tracer_newlon_U_V[jj] = lonchange([floatlat[jj],floatlon[jj],tracerX_km[jj]]) for jj in range(0,arraylen): magnitude_U_V[jj] = pos2dist([floatlon[jj],floatlat[jj],tracer_newlon_U_V[jj],tracer_newlat_U_V[jj]]) difnorm = np.zeros([arraylen]) normtrac = np.zeros([arraylen]) for jj in range(0,arraylen): normtrac[jj] = np.sqrt(tracerX_km[jj]**2 + tracerY_km[jj]**2) difnorm[jj] = magnitude_U_V[jj]-normtrac[jj] # # =========================================================== # Tracer bearing and magnitude # =========================================================== # origtotracer = np.zeros([arraylen-1]) origtofloat = np.zeros([arraylen-1]) tracertofloat = np.zeros([arraylen-1]) tracerfloatratio = np.zeros([arraylen-1]) # for jj in range(0,arraylen-1): origtotracer[jj] = pos2dist([floatlon[jj],floatlat[jj], tracer_newlon_U_V[jj], tracer_newlat_U_V[jj]]) origtofloat[jj] = pos2dist([floatlon[jj],floatlat[jj], floatlon[jj+1],floatlat[jj+1]]) tracertofloat[jj] = pos2dist([floatlon[jj+1],floatlat[jj+1], tracer_newlon_U_V[jj], tracer_newlat_U_V[jj]]) tracerfloatratio[jj] = abs(tracertofloat[jj])/origtotracer[jj] # # =========================================================== # # LINEARLY INTERPOLATING PROFILES TO KNOWN VERTICAL GRID # # =========================================================== # # os.chdir('/data/orbprocess_mail/alex/Jan01_Jun04_2013/data/type/noday71/') rowslen = np.zeros(arraylen) for i in range(0,arraylen): dataout, rowslen[i], cols = csvread(fullnames[i]) print dataout[0,0] profilelen = np.amin(rowslen) profilelen = int(profilelen) print rowslen #rowslen = np.zeros(arraylen) RawData = np.zeros([profilelen,10,arraylen]) for i in range(0,arraylen): datas, rows, cols = csvread(fullnames[i]) if rows == 120: for j in range(1,120): RawData[j-1,0,i] = datas[0,j] # Depth RawData[j-1,1,i] = datas[1,j] # Temperature RawData[j-1,2,i] = datas[2,j] # Density RawData[j-1,3,i] = datas[3,j] # Salinity RawData[j-1,4,i] = datas[4,j] # Pressure RawData[j-1,5,i] = datas[7,j] # Oxygen RawData[j-1,6,i] = datas[11,j] # Chlorophyll RawData[j-1,7,i] = datas[12,j] # Backscatter RawData[j-1,8,i] = datas[13,j] # CDOM RawData[j-1,9,i] = datas[2,j] - datas[2,profilelen-1] # Density minus sfc density if rows == 119: for j in range(0,119): RawData[j,0,i] = datas[0,j] # Depth RawData[j,1,i] = datas[1,j] # Temperature RawData[j,2,i] = datas[2,j] # Density RawData[j,3,i] = datas[3,j] # Salinity RawData[j,4,i] = datas[4,j] # Pressure RawData[j,5,i] = datas[7,j] # Oxygen RawData[j,6,i] = datas[11,j] # Chlorophyll RawData[j,7,i] = datas[12,j] # Backscatter RawData[j,8,i] = datas[13,j] # CDOM RawData[j,9,i] = datas[2,j] - datas[2,profilelen-1] # Density minus sfc density # #rowslen = np.zeros(arraylen) RawDataAvg3 = np.zeros([profilelen-2,10,arraylen]) RawDataAvg5 = np.zeros([profilelen-4,10,arraylen]) for i in range(0,arraylen): # # Filtering over 3 points for j in range(profilelen-2): RawDataAvg3[j,0,i] = (RawData[j,0,i] + RawData[j+1,0,i] + RawData[j+2,0,i])/3 RawDataAvg3[j,1,i] = (RawData[j,1,i] + RawData[j+1,1,i] + RawData[j+2,1,i])/3 RawDataAvg3[j,2,i] = (RawData[j,2,i] + RawData[j+1,2,i] + RawData[j+2,2,i])/3 RawDataAvg3[j,3,i] = (RawData[j,3,i] + RawData[j+1,3,i] + RawData[j+2,3,i])/3 RawDataAvg3[j,4,i] = (RawData[j,4,i] + RawData[j+1,4,i] + RawData[j+2,4,i])/3 RawDataAvg3[j,5,i] = (RawData[j,5,i] + RawData[j+1,5,i] + RawData[j+2,5,i])/3 RawDataAvg3[j,6,i] = (RawData[j,6,i] + RawData[j+1,6,i] + RawData[j+2,6,i])/3 RawDataAvg3[j,7,i] = (RawData[j,7,i] + RawData[j+1,7,i] + RawData[j+2,7,i])/3 RawDataAvg3[j,8,i] = (RawData[j,8,i] + RawData[j+1,8,i] + RawData[j+2,8,i])/3 RawDataAvg3[j,9,i] = (RawData[j,9,i] + RawData[j+1,9,i] + RawData[j+2,9,i])/3 # # Filtering over 5 points for j in range(profilelen-4): RawDataAvg5[j,0,i] = (RawData[j,0,i] + RawData[j+1,0,i] + RawData[j+2,0,i] + RawData[j+3,0,i] + RawData[j+4,0,i])/5 RawDataAvg5[j,1,i] = (RawData[j,1,i] + RawData[j+1,1,i] + RawData[j+2,1,i] + RawData[j+3,1,i] + RawData[j+4,1,i])/5 RawDataAvg5[j,2,i] = (RawData[j,2,i] + RawData[j+1,2,i] + RawData[j+2,2,i] + RawData[j+3,2,i] + RawData[j+4,2,i])/5 RawDataAvg5[j,3,i] = (RawData[j,3,i] + RawData[j+1,3,i] + RawData[j+2,3,i] + RawData[j+3,3,i] + RawData[j+4,3,i])/5 RawDataAvg5[j,4,i] = (RawData[j,4,i] + RawData[j+1,4,i] + RawData[j+2,4,i] + RawData[j+3,4,i] + RawData[j+4,4,i])/5 RawDataAvg5[j,5,i] = (RawData[j,5,i] + RawData[j+1,5,i] + RawData[j+2,5,i] + RawData[j+3,5,i] + RawData[j+4,5,i])/5 RawDataAvg5[j,6,i] = (RawData[j,6,i] + RawData[j+1,6,i] + RawData[j+2,6,i] + RawData[j+3,6,i] + RawData[j+4,6,i])/5 RawDataAvg5[j,7,i] = (RawData[j,7,i] + RawData[j+1,7,i] + RawData[j+2,7,i] + RawData[j+3,7,i] + RawData[j+4,7,i])/5 RawDataAvg5[j,8,i] = (RawData[j,8,i] + RawData[j+1,8,i] + RawData[j+2,8,i] + RawData[j+3,8,i] + RawData[j+4,8,i])/5 RawDataAvg5[j,9,i] = (RawData[j,9,i] + RawData[j+1,9,i] + RawData[j+2,9,i] + RawData[j+3,9,i] + RawData[j+4,9,i])/5 # # =========================================================== # # SUB-SETTING PROFILES by 400 CM2/S2 # # =========================================================== # floatMKEle400 = np.array([value for value in floatMKE if value <= 400.0]) floatMKEg400 = np.array([value for value in floatMKE if value > 400.0]) lenfloatMKEle400 = float(len(floatMKEle400)) lenfloatMKEg400 = float(len(floatMKEg400)) RawData_g400 = np.zeros([profilelen,10,lenfloatMKEg400]) RawData_le400 = np.zeros([profilelen,10,lenfloatMKEle400]) RawDataAvg3_g400 = np.zeros([profilelen-2,10,lenfloatMKEg400]) RawDataAvg3_le400 = np.zeros([profilelen-2,10,lenfloatMKEle400]) RawDataAvg5_g400 = np.zeros([profilelen-4,10,lenfloatMKEg400]) RawDataAvg5_le400 = np.zeros([profilelen-4,10,lenfloatMKEle400]) floatdate_le400 = np.zeros(lenfloatMKEle400) floatdate_g400 = np.zeros(lenfloatMKEg400) ii = 0 iii = 0 for i in range(0,arraylen): if floatMKE[i] <= 400.0: RawData_le400[:,0,ii] = RawData[:,0,i] RawData_le400[:,1,ii] = RawData[:,1,i] RawData_le400[:,2,ii] = RawData[:,2,i] RawData_le400[:,3,ii] = RawData[:,3,i] RawData_le400[:,4,ii] = RawData[:,4,i] RawData_le400[:,5,ii] = RawData[:,5,i] RawData_le400[:,6,ii] = RawData[:,6,i] RawData_le400[:,7,ii] = RawData[:,7,i] RawData_le400[:,8,ii] = RawData[:,8,i] RawData_le400[:,9,ii] = RawData[:,9,i] RawDataAvg3_le400[:,0,ii] = RawDataAvg3[:,0,i] RawDataAvg3_le400[:,1,ii] = RawDataAvg3[:,1,i] RawDataAvg3_le400[:,2,ii] = RawDataAvg3[:,2,i] RawDataAvg3_le400[:,3,ii] = RawDataAvg3[:,3,i] RawDataAvg3_le400[:,4,ii] = RawDataAvg3[:,4,i] RawDataAvg3_le400[:,5,ii] = RawDataAvg3[:,5,i] RawDataAvg3_le400[:,6,ii] = RawDataAvg3[:,6,i] RawDataAvg3_le400[:,7,ii] = RawDataAvg3[:,7,i] RawDataAvg3_le400[:,8,ii] = RawDataAvg3[:,8,i] RawDataAvg3_le400[:,9,ii] = RawDataAvg3[:,9,i] RawDataAvg5_le400[:,0,ii] = RawDataAvg5[:,0,i] RawDataAvg5_le400[:,1,ii] = RawDataAvg5[:,1,i] RawDataAvg5_le400[:,2,ii] = RawDataAvg5[:,2,i] RawDataAvg5_le400[:,3,ii] = RawDataAvg5[:,3,i] RawDataAvg5_le400[:,4,ii] = RawDataAvg5[:,4,i] RawDataAvg5_le400[:,5,ii] = RawDataAvg5[:,5,i] RawDataAvg5_le400[:,6,ii] = RawDataAvg5[:,6,i] RawDataAvg5_le400[:,7,ii] = RawDataAvg5[:,7,i] RawDataAvg5_le400[:,8,ii] = RawDataAvg5[:,8,i] RawDataAvg5_le400[:,9,ii] = RawDataAvg5[:,9,i] floatdate_le400[ii] = floatdate[i] ii = ii + 1 else: RawData_g400[:,0,iii] = RawData[:,0,i] RawData_g400[:,1,iii] = RawData[:,1,i] RawData_g400[:,2,iii] = RawData[:,2,i] RawData_g400[:,3,iii] = RawData[:,3,i] RawData_g400[:,4,iii] = RawData[:,4,i] RawData_g400[:,5,iii] = RawData[:,5,i] RawData_g400[:,6,iii] = RawData[:,6,i] RawData_g400[:,7,iii] = RawData[:,7,i] RawData_g400[:,8,iii] = RawData[:,8,i] RawData_g400[:,9,iii] = RawData[:,9,i] RawDataAvg3_g400[:,0,iii] = RawDataAvg3[:,0,i] RawDataAvg3_g400[:,1,iii] = RawDataAvg3[:,1,i] RawDataAvg3_g400[:,2,iii] = RawDataAvg3[:,2,i] RawDataAvg3_g400[:,3,iii] = RawDataAvg3[:,3,i] RawDataAvg3_g400[:,4,iii] = RawDataAvg3[:,4,i] RawDataAvg3_g400[:,5,iii] = RawDataAvg3[:,5,i] RawDataAvg3_g400[:,6,iii] = RawDataAvg3[:,6,i] RawDataAvg3_g400[:,7,iii] = RawDataAvg3[:,7,i] RawDataAvg3_g400[:,8,iii] = RawDataAvg3[:,8,i] RawDataAvg3_g400[:,9,iii] = RawDataAvg3[:,9,i] RawDataAvg5_g400[:,0,iii] = RawDataAvg5[:,0,i] RawDataAvg5_g400[:,1,iii] = RawDataAvg5[:,1,i] RawDataAvg5_g400[:,2,iii] = RawDataAvg5[:,2,i] RawDataAvg5_g400[:,3,iii] = RawDataAvg5[:,3,i] RawDataAvg5_g400[:,4,iii] = RawDataAvg5[:,4,i] RawDataAvg5_g400[:,5,iii] = RawDataAvg5[:,5,i] RawDataAvg5_g400[:,6,iii] = RawDataAvg5[:,6,i] RawDataAvg5_g400[:,7,iii] = RawDataAvg5[:,7,i] RawDataAvg5_g400[:,8,iii] = RawDataAvg5[:,8,i] RawDataAvg5_g400[:,9,iii] = RawDataAvg5[:,9,i] floatdate_g400[iii] = floatdate[i] iii = iii + 1 # # =========================================================== # # Derivatives # # =========================================================== # lenfloatMKEle400 = int(lenfloatMKEle400) lenfloatMKEg400 = int(lenfloatMKEg400) # drhodz_RawData_le400 = np.zeros([profilelen-1,2,lenfloatMKEle400]) drhodz_depth_RawData_le400 = np.zeros([profilelen-1,lenfloatMKEle400]) drhodz_RawData_g400 = np.zeros([profilelen-1,2,lenfloatMKEg400]) drhodz_depth_RawData_g400 = np.zeros([profilelen-1,lenfloatMKEg400]) # drhodz_RawDataAvg3_le400 = np.zeros([profilelen-3,2,lenfloatMKEle400]) drhodz_depth_RawDataAvg3_le400 = np.zeros([profilelen-3,lenfloatMKEle400]) drhodz_RawDataAvg3_g400 = np.zeros([profilelen-3,2,lenfloatMKEg400]) drhodz_depth_RawDataAvg3_g400 = np.zeros([profilelen-3,lenfloatMKEg400]) # drhodz_RawDataAvg5_le400 = np.zeros([profilelen-5,2,lenfloatMKEle400]) drhodz_depth_RawDataAvg5_le400 = np.zeros([profilelen-5,lenfloatMKEle400]) drhodz_RawDataAvg5_g400 = np.zeros([profilelen-5,2,lenfloatMKEg400]) drhodz_depth_RawDataAvg5_g400 = np.zeros([profilelen-5,lenfloatMKEg400]) # # <= 400 for i in range(0,lenfloatMKEle400): # # drhodz for j in range(profilelen-1): drhodz_RawData_le400[j,0,i] = (RawData_le400[j+1,2,i] - RawData_le400[j,2,i])/(RawData_le400[j+1,0,i] - RawData_le400[j,0,i]) drhodz_RawData_le400[j,1,i] = (RawData_le400[j+1,9,i] - RawData_le400[j,9,i])/(RawData_le400[j+1,0,i] - RawData_le400[j,0,i]) drhodz_depth_RawData_le400[j,i] = (RawData_le400[j+1,0,i] + RawData_le400[j,0,i])/2 # # drhodz 3 point filter for j in range(profilelen-3): drhodz_RawDataAvg3_le400[j,0,i] = (RawDataAvg3_le400[j+1,2,i] - RawDataAvg3_le400[j,2,i])/(RawDataAvg3_le400[j+1,0,i] - RawDataAvg3_le400[j,0,i]) drhodz_RawDataAvg3_le400[j,1,i] = (RawDataAvg3_le400[j+1,9,i] - RawDataAvg3_le400[j,9,i])/(RawDataAvg3_le400[j+1,0,i] - RawDataAvg3_le400[j,0,i]) drhodz_depth_RawDataAvg3_le400[j,i] = (RawDataAvg3_le400[j+1,0,i] + RawDataAvg3_le400[j,0,i])/2 # # drhodz 5 point filter for j in range(profilelen-5): drhodz_RawDataAvg5_le400[j,0,i] = (RawDataAvg5_le400[j+1,2,i] - RawDataAvg5_le400[j,2,i])/(RawDataAvg5_le400[j+1,0,i] - RawDataAvg5_le400[j,0,i]) drhodz_RawDataAvg5_le400[j,1,i] = (RawDataAvg5_le400[j+1,9,i] - RawDataAvg5_le400[j,9,i])/(RawDataAvg5_le400[j+1,0,i] - RawDataAvg5_le400[j,0,i]) drhodz_depth_RawDataAvg5_le400[j,i] = (RawDataAvg5_le400[j+1,0,i] + RawDataAvg5_le400[j,0,i])/2 # # > 400 for i in range(0,lenfloatMKEg400): # # drhodz for j in range(profilelen-1): drhodz_RawData_g400[j,0,i] = (RawData_g400[j+1,2,i] - RawData_g400[j,2,i])/(RawData_g400[j+1,0,i] - RawData_g400[j,0,i]) drhodz_RawData_g400[j,1,i] = (RawData_g400[j+1,9,i] - RawData_g400[j,9,i])/(RawData_g400[j+1,0,i] - RawData_g400[j,0,i]) drhodz_depth_RawData_g400[j,i] = (RawData_g400[j+1,0,i] + RawData_g400[j,0,i])/2 # # drhodz 3 point filter for j in range(profilelen-3): drhodz_RawDataAvg3_g400[j,0,i] = (RawDataAvg3_g400[j+1,2,i] - RawDataAvg3_g400[j,2,i])/(RawDataAvg3_g400[j+1,0,i] - RawDataAvg3_g400[j,0,i]) drhodz_RawDataAvg3_g400[j,1,i] = (RawDataAvg3_g400[j+1,9,i] - RawDataAvg3_g400[j,9,i])/(RawDataAvg3_g400[j+1,0,i] - RawDataAvg3_g400[j,0,i]) drhodz_depth_RawDataAvg3_g400[j,i] = (RawDataAvg3_g400[j+1,0,i] + RawDataAvg3_g400[j,0,i])/2 # # drhodz 5 point filter for j in range(profilelen-5): drhodz_RawDataAvg5_g400[j,0,i] = (RawDataAvg5_g400[j+1,2,i] - RawDataAvg5_g400[j,2,i])/(RawDataAvg5_g400[j+1,0,i] - RawDataAvg5_g400[j,0,i]) drhodz_RawDataAvg5_g400[j,1,i] = (RawDataAvg5_g400[j+1,9,i] - RawDataAvg5_g400[j,9,i])/(RawDataAvg5_g400[j+1,0,i] - RawDataAvg5_g400[j,0,i]) drhodz_depth_RawDataAvg5_g400[j,i] = (RawDataAvg5_g400[j+1,0,i] + RawDataAvg5_g400[j,0,i])/2 # # =========================================================== # # filtered and base profile AVG AND STDEV'S # # =========================================================== # RawData_le400Avg = np.zeros([profilelen,10]) RawData_le400Std = np.zeros([profilelen,10]) RawData_g400Avg = np.zeros([profilelen,10]) RawData_g400Std = np.zeros([profilelen,10]) for j in range(0,profilelen): for k in range(0,10): RawData_le400Avg[j,k] = np.mean(RawData_le400[j,k,:]) RawData_le400Std[j,k] = np.std(RawData_le400[j,k,:]) RawData_g400Avg[j,k] = np.mean(RawData_g400[j,k,:]) RawData_g400Std[j,k] = np.std(RawData_g400[j,k,:]) # RawDataAvg3_le400Avg = np.zeros([profilelen-2,10]) RawDataAvg3_le400Std = np.zeros([profilelen-2,10]) RawDataAvg3_g400Avg = np.zeros([profilelen-2,10]) RawDataAvg3_g400Std = np.zeros([profilelen-2,10]) for j in range(0,profilelen-2): for k in range(0,10): RawDataAvg3_le400Avg[j,k] = np.mean(RawDataAvg3_le400[j,k,:]) RawDataAvg3_le400Std[j,k] = np.std(RawDataAvg3_le400[j,k,:]) RawDataAvg3_g400Avg[j,k] = np.mean(RawDataAvg3_g400[j,k,:]) RawDataAvg3_g400Std[j,k] = np.std(RawDataAvg3_g400[j,k,:]) # RawDataAvg5_le400Avg = np.zeros([profilelen-4,10]) RawDataAvg5_le400Std = np.zeros([profilelen-4,10]) RawDataAvg5_g400Avg = np.zeros([profilelen-4,10]) RawDataAvg5_g400Std = np.zeros([profilelen-4,10]) for j in range(0,profilelen-4): for k in range(0,10): RawDataAvg5_le400Avg[j,k] = np.mean(RawDataAvg5_le400[j,k,:]) RawDataAvg5_le400Std[j,k] = np.std(RawDataAvg5_le400[j,k,:]) RawDataAvg5_g400Avg[j,k] = np.mean(RawDataAvg5_g400[j,k,:]) RawDataAvg5_g400Std[j,k] = np.std(RawDataAvg5_g400[j,k,:]) # # =========================================================== # # drhodz AVG AND STDEV'S # # =========================================================== # drhodz_RawData_le400Avg = np.zeros([profilelen-1,2]) drhodz_depth_RawData_le400Avg = np.zeros([profilelen-1]) drhodz_RawData_g400Avg = np.zeros([profilelen-1,2]) drhodz_depth_RawData_g400Avg = np.zeros([profilelen-1]) # drhodz_RawData_le400Std = np.zeros([profilelen-1,2]) drhodz_depth_RawData_le400Std = np.zeros([profilelen-1]) drhodz_RawData_g400Std = np.zeros([profilelen-1,2]) drhodz_depth_RawData_g400Std = np.zeros([profilelen-1]) # for j in range(0,profilelen-1): # # Mean's drhodz_RawData_le400Avg[j,0] = np.mean(drhodz_RawData_le400[j,0,:]) drhodz_RawData_le400Avg[j,1] = np.mean(drhodz_RawData_le400[j,1,:]) drhodz_depth_RawData_le400Avg[j] = np.mean(drhodz_depth_RawData_le400[j,:]) # drhodz_RawData_g400Avg[j,0] = np.mean(drhodz_RawData_g400[j,0,:]) drhodz_RawData_g400Avg[j,1] = np.mean(drhodz_RawData_g400[j,1,:]) drhodz_depth_RawData_g400Avg[j] = np.mean(drhodz_depth_RawData_g400[j,:]) # # Std's drhodz_RawData_le400Std[j,0] = np.std(drhodz_RawData_le400[j,0,:]) drhodz_RawData_le400Std[j,1] = np.std(drhodz_RawData_le400[j,1,:]) drhodz_depth_RawData_le400Std[j] = np.std(drhodz_depth_RawData_le400[j,:]) # drhodz_RawData_g400Std[j,0] = np.std(drhodz_RawData_g400[j,0,:]) drhodz_RawData_g400Std[j,1] = np.std(drhodz_RawData_g400[j,1,:]) drhodz_depth_RawData_g400Std[j] = np.std(drhodz_depth_RawData_g400[j,:]) # drhodz_RawDataAvg3_le400Avg = np.zeros([profilelen-3,2]) drhodz_depth_RawDataAvg3_le400Avg = np.zeros([profilelen-3]) drhodz_RawDataAvg3_g400Avg = np.zeros([profilelen-3,2]) drhodz_depth_RawDataAvg3_g400Avg = np.zeros([profilelen-3]) # drhodz_RawDataAvg3_le400Std = np.zeros([profilelen-3,2]) drhodz_depth_RawDataAvg3_le400Std = np.zeros([profilelen-3]) drhodz_RawDataAvg3_g400Std = np.zeros([profilelen-3,2]) drhodz_depth_RawDataAvg3_g400Std = np.zeros([profilelen-3]) # for j in range(0,profilelen-3): # # Mean's drhodz_RawDataAvg3_le400Avg[j,0] = np.mean(drhodz_RawDataAvg3_le400[j,0,:]) drhodz_RawDataAvg3_le400Avg[j,1] = np.mean(drhodz_RawDataAvg3_le400[j,1,:]) drhodz_depth_RawDataAvg3_le400Avg[j] = np.mean(drhodz_depth_RawDataAvg3_le400[j,:]) # drhodz_RawDataAvg3_g400Avg[j,0] = np.mean(drhodz_RawDataAvg3_g400[j,0,:]) drhodz_RawDataAvg3_g400Avg[j,1] = np.mean(drhodz_RawDataAvg3_g400[j,1,:]) drhodz_depth_RawDataAvg3_g400Avg[j] = np.mean(drhodz_depth_RawDataAvg3_g400[j,:]) # # Std's drhodz_RawDataAvg3_le400Std[j,0] = np.std(drhodz_RawDataAvg3_le400[j,0,:]) drhodz_RawDataAvg3_le400Std[j,1] = np.std(drhodz_RawDataAvg3_le400[j,1,:]) drhodz_depth_RawDataAvg3_le400Std[j] = np.std(drhodz_depth_RawDataAvg3_le400[j,:]) # drhodz_RawDataAvg3_g400Std[j,0] = np.std(drhodz_RawDataAvg3_g400[j,0,:]) drhodz_RawDataAvg3_g400Std[j,1] = np.std(drhodz_RawDataAvg3_g400[j,1,:]) drhodz_depth_RawDataAvg3_g400Std[j] = np.std(drhodz_depth_RawDataAvg3_g400[j,:]) # drhodz_RawDataAvg5_le400Avg = np.zeros([profilelen-3,2]) drhodz_depth_RawDataAvg5_le400Avg = np.zeros([profilelen-3]) drhodz_RawDataAvg5_g400Avg = np.zeros([profilelen-3,2]) drhodz_depth_RawDataAvg5_g400Avg = np.zeros([profilelen-3]) # drhodz_RawDataAvg5_le400Std = np.zeros([profilelen-5,2]) drhodz_depth_RawDataAvg5_le400Std = np.zeros([profilelen-5]) drhodz_RawDataAvg5_g400Std = np.zeros([profilelen-5,2]) drhodz_depth_RawDataAvg5_g400Std = np.zeros([profilelen-5]) # for j in range(0,profilelen-5): # # Mean's drhodz_RawDataAvg5_le400Avg[j,0] = np.mean(drhodz_RawDataAvg5_le400[j,0,:]) drhodz_RawDataAvg5_le400Avg[j,1] = np.mean(drhodz_RawDataAvg5_le400[j,1,:]) drhodz_depth_RawDataAvg5_le400Avg[j] = np.mean(drhodz_depth_RawDataAvg5_le400[j,:]) # drhodz_RawDataAvg5_g400Avg[j,0] = np.mean(drhodz_RawDataAvg5_g400[j,0,:]) drhodz_RawDataAvg5_g400Avg[j,1] = np.mean(drhodz_RawDataAvg5_g400[j,1,:]) drhodz_depth_RawDataAvg5_g400Avg[j] = np.mean(drhodz_depth_RawDataAvg5_g400[j,:]) # # Std's drhodz_RawDataAvg5_le400Std[j,0] = np.std(drhodz_RawDataAvg5_le400[j,0,:]) drhodz_RawDataAvg5_le400Std[j,1] = np.std(drhodz_RawDataAvg5_le400[j,1,:]) drhodz_depth_RawDataAvg5_le400Std[j] = np.std(drhodz_depth_RawDataAvg5_le400[j,:]) # drhodz_RawDataAvg5_g400Std[j,0] = np.std(drhodz_RawDataAvg5_g400[j,0,:]) drhodz_RawDataAvg5_g400Std[j,1] = np.std(drhodz_RawDataAvg5_g400[j,1,:]) drhodz_depth_RawDataAvg5_g400Std[j] = np.std(drhodz_depth_RawDataAvg5_g400[j,:]) # # =========================================================== # # Plotting # # =========================================================== # # # =========================================================== # avg drho/dz raw mean profiles # =========================================================== # from matplotlib.font_manager import FontProperties legendfont = FontProperties() legendfont.set_name('Computer Modern Roman') legendfont.set_size('x-small') rcParams['axes.labelsize'] = 18 rcParams['xtick.labelsize'] = 18 rcParams['ytick.labelsize'] = 18 rcParams['legend.fontsize'] = 14 # from matplotlib import rcParams rcParams['font.family'] = 'serif' rcParams['font.serif'] = ['Computer Modern Roman'] rcParams['text.usetex'] = True from pylab import * # # Plot Set-up import math as ma from mpl_toolkits.axes_grid1 import make_axes_locatable fig = plt.figure() ax = fig.add_axes([0.15,0.1,0.68,0.85]) # ax.plot(abs(drhodz_RawData_g400Avg[:,1]),drhodz_depth_RawData_g400Avg,'k', linewidth = 3, zorder = 1) ax.plot(abs(drhodz_RawData_le400Avg[:,1]),drhodz_depth_RawData_le400Avg,'g', linewidth = 3, zorder = 2) # ax.set_xlabel(r'drho/dz (kgm$^{-4}$)') ax.set_ylabel(r'Depth (m)') ax.set_ylim([-300,0]) ax.set_xlim([-.005,.02]) #axvline((1027.45-1026.), linewidth=3, color='b',linestyle=':', zorder = 2) # # Save and spc to storm plt.savefig('/home/ardavies/satdata/OSCAR/drhodz_means_profile.png') # # =========================================================== # avg drho/dz raw mean profiles 5Avg # =========================================================== # from matplotlib.font_manager import FontProperties legendfont = FontProperties() legendfont.set_name('Computer Modern Roman') legendfont.set_size('x-small') rcParams['axes.labelsize'] = 18 rcParams['xtick.labelsize'] = 18 rcParams['ytick.labelsize'] = 18 rcParams['legend.fontsize'] = 14 # from matplotlib import rcParams rcParams['font.family'] = 'serif' rcParams['font.serif'] = ['Computer Modern Roman'] rcParams['text.usetex'] = True from pylab import * # # Plot Set-up import math as ma from mpl_toolkits.axes_grid1 import make_axes_locatable fig = plt.figure() ax = fig.add_axes([0.15,0.1,0.68,0.85]) # ax.plot(abs(drhodz_RawDataAvg3_g400Avg[:,1]),drhodz_depth_RawDataAvg3_g400Avg,'k', linewidth = 3, zorder = 1) ax.plot(abs(drhodz_RawDataAvg3_le400Avg[:,1]),drhodz_depth_RawDataAvg3_le400Avg,'g', linewidth = 3, zorder = 2) ax.set_xlabel(r'drho/dz (kgm$^{-4}$)') ax.set_ylabel(r'Depth (m)') ax.set_ylim([-300,0]) ax.set_xlim([-.005,.02]) #axvline((1027.45-1026.), linewidth=3, color='b',linestyle=':', zorder = 2) # # Save and spc to storm plt.savefig('/home/ardavies/satdata/OSCAR/drhodz_meansAvg3_profile.png') # # =========================================================== # avg drho/dz raw mean profiles 5Avg # =========================================================== # from matplotlib.font_manager import FontProperties legendfont = FontProperties() legendfont.set_name('Computer Modern Roman') legendfont.set_size('x-small') rcParams['axes.labelsize'] = 18 rcParams['xtick.labelsize'] = 18 rcParams['ytick.labelsize'] = 18 rcParams['legend.fontsize'] = 14 # from matplotlib import rcParams rcParams['font.family'] = 'serif' rcParams['font.serif'] = ['Computer Modern Roman'] rcParams['text.usetex'] = True from pylab import * # # Plot Set-up import math as ma from mpl_toolkits.axes_grid1 import make_axes_locatable fig = plt.figure() ax = fig.add_axes([0.15,0.1,0.68,0.85]) # ax.plot(abs(drhodz_RawDataAvg5_g400Avg[:,1]),drhodz_depth_RawDataAvg5_g400Avg,'k', linewidth = 3, zorder = 1) ax.plot(abs(drhodz_RawDataAvg5_le400Avg[:,1]),drhodz_depth_RawDataAvg5_le400Avg,'g', linewidth = 3, zorder = 2) ax.set_xlabel(r'drho/dz (kgm$^{-4}$)') ax.set_ylabel(r'Depth (m)') ax.set_ylim([-300,0]) ax.set_xlim([-.005,.02]) #axvline((1027.45-1026.), linewidth=3, color='b',linestyle=':', zorder = 2) # # Save and spc to storm plt.savefig('/home/ardavies/satdata/OSCAR/drhodz_meansAvg5_profile.png') # # =========================================================== # avg drho/dz raw ind. profiles # =========================================================== # pp1 = PdfPages('/home/ardavies/satdata/OSCAR/drhodz_indiv_le400_profile.pdf') # for i in range(0,lenfloatMKEle400): # # =========================================================== # Plotting Density Profiles and Filters (full profiiles) # =========================================================== # # Font Set-up import matplotlib.pyplot as plt import numpy as np import math as ma from mpl_toolkits.axes_grid1 import make_axes_locatable from matplotlib.font_manager import FontProperties legendfont = FontProperties() legendfont.set_name('Computer Modern Roman') legendfont.set_size('x-small') rcParams['axes.labelsize'] = 18 rcParams['xtick.labelsize'] = 18 rcParams['ytick.labelsize'] = 18 rcParams['legend.fontsize'] = 14 # from matplotlib import rcParams rcParams['font.family'] = 'serif' rcParams['font.serif'] = ['Computer Modern Roman'] rcParams['text.usetex'] = True # # fig = plt.figure() # ax = fig.add_subplot(1, 1, 1) # fig = plt.figure() host = fig.add_subplot(111) from pylab import * # # Sharing y axis par1 = host.twiny() # # Plotting p1 = host.plot(abs(drhodz_RawData_le400[:,1,i]),drhodz_depth_RawData_le400[:,i],'k',linewidth = 3) p2 = par1.plot(RawData_le400[:,6,i],RawData_le400[:,0,i],'g',linewidth = 1) # # Setting-up axis host.set_xlim([-.005,.02]) host.set_ylim([-300,0]) par1.set_xlim([0,6]) # # # Ploting horizontal lines for derivative depths # host.axhline(xmin=0.0, xmax=0.16,y=thedepthofmax1[i],linewidth=1, color='0.75') # host.axhline(xmin=0.0, xmax=0.14,y=thedepthofmax3[i],linewidth=1, color='m') # host.axhline(xmin=0.0, xmax=0.11,y=thedepthofmax5[i],linewidth=1, color='b') # host.axhline(xmin=0.0, xmax=0.10,y=thedepthofmax7[i],linewidth=1, color='k') # # # host.axhline(xmin=0.84, xmax=1,y=thedepthofmax21[i],linewidth=1, color='0.75') # host.axhline(xmin=0.86, xmax=1,y=thedepthofmax23[i],linewidth=1, color='m') # host.axhline(xmin=0.88, xmax=1,y=thedepthofmax25[i],linewidth=1, color='b') # host.axhline(xmin=0.9, xmax=1,y=thedepthofmax27[i],linewidth=1, color='k') # # axis labels host.set_xlabel(r'drho/dz (kgm$^{-4}$)') host.set_ylabel("Depth") par1.set_ylabel("Chlorophyll Fluorescence (micro g/l)") # # fig.suptitle("Date: " + str(floatdate[i])) fig.savefig(pp1, format = "pdf") pp1.close() # # =========================================================== # avg drho/dz raw ind profiles 3Avg # =========================================================== # pp1 = PdfPages('/home/ardavies/satdata/OSCAR/drhodz_indivAvg3_le400_profile.pdf') # for i in range(0,lenfloatMKEle400): # # =========================================================== # Plotting Density Profiles and Filters (full profiiles) # =========================================================== # # Font Set-up import matplotlib.pyplot as plt import numpy as np import math as ma from mpl_toolkits.axes_grid1 import make_axes_locatable from matplotlib.font_manager import FontProperties legendfont = FontProperties() legendfont.set_name('Computer Modern Roman') legendfont.set_size('x-small') rcParams['axes.labelsize'] = 18 rcParams['xtick.labelsize'] = 18 rcParams['ytick.labelsize'] = 18 rcParams['legend.fontsize'] = 14 # from matplotlib import rcParams rcParams['font.family'] = 'serif' rcParams['font.serif'] = ['Computer Modern Roman'] rcParams['text.usetex'] = True # # fig = plt.figure() # ax = fig.add_subplot(1, 1, 1) # fig = plt.figure() host = fig.add_subplot(111) from pylab import * # # Sharing y axis par1 = host.twiny() # # Plotting p1 = host.plot(abs(drhodz_RawDataAvg3_le400[:,1,i]),drhodz_depth_RawDataAvg3_le400[:,i],'k',linewidth = 3) p2 = par1.plot(RawDataAvg3_le400[:,6,i],RawDataAvg3_le400[:,0,i],'g',linewidth = 1) # # Setting-up axis host.set_xlim([-.005,.02]) host.set_ylim([-300,0]) par1.set_xlim([0,6]) # # # Ploting horizontal lines for derivative depths # host.axhline(xmin=0.0, xmax=0.16,y=thedepthofmax1[i],linewidth=1, color='0.75') # host.axhline(xmin=0.0, xmax=0.14,y=thedepthofmax3[i],linewidth=1, color='m') # host.axhline(xmin=0.0, xmax=0.11,y=thedepthofmax5[i],linewidth=1, color='b') # host.axhline(xmin=0.0, xmax=0.10,y=thedepthofmax7[i],linewidth=1, color='k') # # # host.axhline(xmin=0.84, xmax=1,y=thedepthofmax21[i],linewidth=1, color='0.75') # host.axhline(xmin=0.86, xmax=1,y=thedepthofmax23[i],linewidth=1, color='m') # host.axhline(xmin=0.88, xmax=1,y=thedepthofmax25[i],linewidth=1, color='b') # host.axhline(xmin=0.9, xmax=1,y=thedepthofmax27[i],linewidth=1, color='k') # # axis labels host.set_xlabel(r'drho/dz (kgm$^{-4}$)') host.set_ylabel("Depth") par1.set_ylabel("Chlorophyll Fluorescence (micro g/l)") # # fig.suptitle("Date: " + str(floatdate[i])) fig.savefig(pp1, format = "pdf") pp1.close() # # =========================================================== # avg drho/dz raw ind. profiles 3Avg # =========================================================== # pp1 = PdfPages('/home/ardavies/satdata/OSCAR/drhodz_indivAvg5_le400_profile.pdf') # for i in range(0,lenfloatMKEle400): # # =========================================================== # Plotting Density Profiles and Filters (full profiiles) # =========================================================== # # Font Set-up import matplotlib.pyplot as plt import numpy as np import math as ma from mpl_toolkits.axes_grid1 import make_axes_locatable from matplotlib.font_manager import FontProperties legendfont = FontProperties() legendfont.set_name('Computer Modern Roman') legendfont.set_size('x-small') rcParams['axes.labelsize'] = 18 rcParams['xtick.labelsize'] = 18 rcParams['ytick.labelsize'] = 18 rcParams['legend.fontsize'] = 14 # from matplotlib import rcParams rcParams['font.family'] = 'serif' rcParams['font.serif'] = ['Computer Modern Roman'] rcParams['text.usetex'] = True # # fig = plt.figure() # ax = fig.add_subplot(1, 1, 1) # fig = plt.figure() host = fig.add_subplot(111) from pylab import * # # Sharing y axis par1 = host.twiny() # # Plotting p1 = host.plot(abs(drhodz_RawDataAvg5_le400[:,1,i]),drhodz_depth_RawDataAvg5_le400[:,i],'k',linewidth = 3) p2 = par1.plot(RawDataAvg5_le400[:,6,i],RawDataAvg5_le400[:,0,i],'g',linewidth = 1) # # Setting-up axis host.set_xlim([-.005,.02]) host.set_ylim([-300,0]) par1.set_xlim([0,6]) # # # Ploting horizontal lines for derivative depths # host.axhline(xmin=0.0, xmax=0.16,y=thedepthofmax1[i],linewidth=1, color='0.75') # host.axhline(xmin=0.0, xmax=0.14,y=thedepthofmax3[i],linewidth=1, color='m') # host.axhline(xmin=0.0, xmax=0.11,y=thedepthofmax5[i],linewidth=1, color='b') # host.axhline(xmin=0.0, xmax=0.10,y=thedepthofmax7[i],linewidth=1, color='k') # # # host.axhline(xmin=0.84, xmax=1,y=thedepthofmax21[i],linewidth=1, color='0.75') # host.axhline(xmin=0.86, xmax=1,y=thedepthofmax23[i],linewidth=1, color='m') # host.axhline(xmin=0.88, xmax=1,y=thedepthofmax25[i],linewidth=1, color='b') # host.axhline(xmin=0.9, xmax=1,y=thedepthofmax27[i],linewidth=1, color='k') # # axis labels host.set_xlabel(r'drho/dz (kgm$^{-4}$)') host.set_ylabel("Depth") par1.set_ylabel("Chlorophyll Fluorescence (micro g/l)") # # fig.suptitle("Date: " + str(floatdate[i])) fig.savefig(pp1, format = "pdf") pp1.close() dsgfds fdshfjksdah gridWdat = 1 if gridWdat == 1: # # =========================================================== # Read the Float Data and Interpolate for Gridding # =========================================================== # os.chdir('/data/orbprocess_mail/alex/Jan01_Jun04_2013/data/type/noday71/') rowslen = np.zeros(arraylen) for i in range(0,arraylen): dataout, rows, cols = csvread(fullnames[i]) Ygrid = np.linspace(-1800,-10,1791*2) interplength = len(Ygrid) GridData = np.zeros([interplength,7,arraylen]) ChlyDumby = np.zeros([interplength,arraylen]) for i in range(0,arraylen): dataout, rows, cols = csvread(fullnames[i]) # TemInterp = interpolate.interp1d(dataout[0,:], dataout[1,:],kind='linear') DenInterp = interpolate.interp1d(dataout[0,:], dataout[2,:],kind='linear') SalInterp = interpolate.interp1d(dataout[0,:], dataout[3,:],kind='linear') PreInterp = interpolate.interp1d(dataout[0,:], dataout[4,:],kind='linear') OxyInterp = interpolate.interp1d(dataout[0,:], dataout[7,:],kind='linear') ChlInterp = interpolate.interp1d(dataout[0,:], dataout[11,:],kind='linear') BskInterp = interpolate.interp1d(dataout[0,:], dataout[12,:],kind='linear') CDMInterp = interpolate.interp1d(dataout[0,:], dataout[13,:],kind='linear') # for j in range(0,interplength): GridData[j,0,i] = TemInterp(Ygrid[j]) GridData[j,1,i] = DenInterp(Ygrid[j]) GridData[j,2,i] = SalInterp(Ygrid[j]) GridData[j,3,i] = OxyInterp(Ygrid[j]) ChlyDumby[j,i] = ChlInterp(Ygrid[j]) GridData[j,5,i] = BskInterp(Ygrid[j]) GridData[j,6,i] = CDMInterp(Ygrid[j]) if (ChlInterp(Ygrid[j]) <= 0.0): GridData[j,4,i] = GridData[j-1,4,i] else: GridData[j,4,i] = ChlyDumby[j,i] # # Average density along single depth avgdenatdepth = np.zeros(interplength) DemeanDensity = np.zeros([interplength,arraylen]) for j in range(0,interplength): avgdenatdepth[j] = np.mean(GridData[j,1,:]) for i in range(0,arraylen): for j in range(0,interplength): DemeanDensity[j,i] = GridData[j,1,i] - avgdenatdepth[j] # # =========================================================== # # SUB-SETTING PROFILES by 400 CM2/S2 # # =========================================================== # floatMKEle400 = np.array([value for value in floatMKE if value <= 400.0]) floatMKEg400 = np.array([value for value in floatMKE if value > 400.0]) lenfloatMKEle400 = float(len(floatMKEle400)) lenfloatMKEg400 = float(len(floatMKEg400)) griddedprofiles_g400 = np.zeros([interplength,8,lenfloatMKEg400]) griddedprofiles_le400 = np.zeros([interplength,8,lenfloatMKEle400]) floatdate_le400 = np.zeros(lenfloatMKEle400) floatdate_g400 = np.zeros(lenfloatMKEg400) ii = 0 iii = 0 for i in range(0,arraylen): if floatMKE[i] <= 400.0: griddedprofiles_le400[:,0,ii] = GridData[:,0,i] griddedprofiles_le400[:,1,ii] = GridData[:,1,i] griddedprofiles_le400[:,2,ii] = GridData[:,2,i] griddedprofiles_le400[:,3,ii] = GridData[:,3,i] griddedprofiles_le400[:,4,ii] = GridData[:,4,i] griddedprofiles_le400[:,5,ii] = GridData[:,5,i] griddedprofiles_le400[:,6,ii] = GridData[:,6,i] griddedprofiles_le400[:,7,ii] = GridData[:,1,i] - GridData[interplength-1,1,i] floatdate_le400[ii] = floatdate[i] ii = ii + 1 else: griddedprofiles_g400[:,0,iii] = GridData[:,0,i] griddedprofiles_g400[:,1,iii] = GridData[:,1,i] griddedprofiles_g400[:,2,iii] = GridData[:,2,i] griddedprofiles_g400[:,3,iii] = GridData[:,3,i] griddedprofiles_g400[:,4,iii] = GridData[:,4,i] griddedprofiles_g400[:,5,iii] = GridData[:,5,i] griddedprofiles_g400[:,6,iii] = GridData[:,6,i] griddedprofiles_g400[:,7,iii] = GridData[:,1,i] - GridData[interplength-1,1,i] floatdate_g400[iii] = floatdate[i] iii = iii + 1 # # =========================================================== # # MEAN and STDEV PROFILES # # =========================================================== # meanprofilele400 = np.zeros(interplength) stdprofilele400 = np.zeros(interplength) meanprofileg400 = np.zeros(interplength) stdprofileg400 = np.zeros(interplength) meanprofilele400_plusstd = np.zeros(interplength) meanprofilele400_minusstd = np.zeros(interplength) meanprofileg400_plusstd = np.zeros(interplength) meanprofileg400_minusstd = np.zeros(interplength) for j in range(0,interplength): meanprofilele400[j] = np.mean(griddedprofiles_le400[j,7,:]) stdprofilele400[j] = np.std(griddedprofiles_le400[j,7,:]) meanprofilele400_plusstd[j] = meanprofilele400[j] + stdprofilele400[j] meanprofilele400_minusstd[j] = meanprofilele400[j] - stdprofilele400[j] # meanprofileg400[j] = np.mean(griddedprofiles_g400[j,7,:]) stdprofileg400[j] = np.std(griddedprofiles_g400[j,7,:]) meanprofileg400_plusstd[j] = meanprofileg400[j] + stdprofileg400[j] meanprofileg400_minusstd[j] = meanprofileg400[j] - stdprofileg400[j] # =========================================================== # # DrhoDz # # =========================================================== # # # > 400 drhodz_meang400 = np.zeros(interplength-1) drhodz_meang400_plusstd = np.zeros(interplength-1) drhodz_meang400_minusstd = np.zeros(interplength-1) drhodz_meang400_2 = np.zeros(interplength-1) drhodz_stdg400 = np.zeros(interplength-1) drhodz_g400 = np.zeros([interplength-1,lenfloatMKEg400]) drhodz_depth = np.zeros(interplength-1) for i in range(0,int(lenfloatMKEg400)): for j in range(0,interplength-1): drhodz_g400[j,i] = (griddedprofiles_g400[j+1,7,i]-griddedprofiles_g400[j,7,i])/(Ygrid[j+1]-Ygrid[j]) for j in range(0,interplength-1): drhodz_meang400[j] = (meanprofileg400[j+1]-meanprofileg400[j])/(Ygrid[j+1]-Ygrid[j]) drhodz_depth[j] = (Ygrid[j+1]+Ygrid[j])/2 drhodz_meang400_2[j] = np.mean(drhodz_g400[j,:]) drhodz_stdg400[j] = np.std(drhodz_g400[j,:]) for j in range(0,interplength-1): drhodz_meang400_plusstd[j] = abs(drhodz_meang400_2[j])+abs(drhodz_stdg400[j]) drhodz_meang400_minusstd[j] = abs(drhodz_meang400_2[j])-abs(drhodz_stdg400[j]) print drhodz_meang400 print drhodz_depth # # <= 400 drhodz_meanle400 = np.zeros(interplength-1) drhodz_meanle400_plusstd = np.zeros(interplength-1) drhodz_meanle400_minusstd = np.zeros(interplength-1) drhodz_meanle400_2 = np.zeros(interplength-1) drhodz_stdle400 = np.zeros(interplength-1) drhodz_le400 = np.zeros([interplength-1,lenfloatMKEle400]) for i in range(0,int(lenfloatMKEle400)): for j in range(0,interplength-1): drhodz_le400[j,i] = (griddedprofiles_le400[j+1,7,i]-griddedprofiles_le400[j,7,i])/(Ygrid[j+1]-Ygrid[j]) for j in range(0,interplength-1): drhodz_meanle400[j] = (meanprofilele400[j+1]-meanprofilele400[j])/(Ygrid[j+1]-Ygrid[j]) drhodz_meanle400_2[j] = np.mean(drhodz_le400[j,:]) drhodz_stdle400[j] = np.std(drhodz_le400[j,:]) for j in range(0,interplength-1): drhodz_meanle400_plusstd[j] = abs(drhodz_meanle400_2[j])+abs(drhodz_stdle400[j]) drhodz_meanle400_minusstd[j] = abs(drhodz_meanle400_2[j])-abs(drhodz_stdle400[j]) # # =========================================================== # # Plotting # # =========================================================== # # from matplotlib.font_manager import FontProperties legendfont = FontProperties() legendfont.set_name('Computer Modern Roman') legendfont.set_size('x-small') rcParams['axes.labelsize'] = 18 rcParams['xtick.labelsize'] = 18 rcParams['ytick.labelsize'] = 18 rcParams['legend.fontsize'] = 14 # from matplotlib import rcParams rcParams['font.family'] = 'serif' rcParams['font.serif'] = ['Computer Modern Roman'] rcParams['text.usetex'] = True # fig = plt.figure() # ax4 = fig.add_axes([0.12,0.12,0.8,0.8]) from pylab import * # # Plot Set-up import math as ma from mpl_toolkits.axes_grid1 import make_axes_locatable fig = plt.figure() ax = fig.add_axes([0.15,0.1,0.68,0.85]) ax.plot(abs(drhodz_meang400),drhodz_depth,'k', linewidth = 3, zorder = 1) ax.plot(abs(drhodz_meanle400),drhodz_depth,'g', linewidth = 3, zorder = 2) ax.set_xlabel(r'drho/dz (kgm$^{-4}$)') ax.set_ylabel(r'Depth (m)') ax.set_ylim([-300,0]) #ax.set_xlim([-.1,1.]) #axvline((1027.45-1026.), linewidth=3, color='b',linestyle=':', zorder = 2) # # Save and spc to storm plt.savefig('/home/ardavies/satdata/OSCAR/floatMKE_gandlemeandrhodz_profile.png') from matplotlib.font_manager import FontProperties legendfont = FontProperties() legendfont.set_name('Computer Modern Roman') legendfont.set_size('x-small') rcParams['axes.labelsize'] = 18 rcParams['xtick.labelsize'] = 18 rcParams['ytick.labelsize'] = 18 rcParams['legend.fontsize'] = 14 # from matplotlib import rcParams rcParams['font.family'] = 'serif' rcParams['font.serif'] = ['Computer Modern Roman'] rcParams['text.usetex'] = True # fig = plt.figure() # ax4 = fig.add_axes([0.12,0.12,0.8,0.8]) from pylab import * # # Plot Set-up import math as ma from mpl_toolkits.axes_grid1 import make_axes_locatable fig = plt.figure() ax = fig.add_axes([0.15,0.1,0.68,0.85]) print len(drhodz_depth) print len(drhodz_meang400_plusstd) print len(drhodz_meang400_minusstd) ax.plot(abs(drhodz_meang400_2),drhodz_depth,'k', linewidth = 3, zorder = 12) # ax.plot(drhodz_meang400_minusstd,drhodz_depth,'k', linewidth = .5, zorder = 10) # ax.plot(drhodz_meang400_plusstd,drhodz_depth,'k', linewidth = .5, zorder = 9) ax.plot(abs(drhodz_meanle400_2),drhodz_depth,'g', linewidth = 3, zorder = 11) # ax.plot(drhodz_meanle400_minusstd,drhodz_depth,'g', linewidth = .5, zorder = 7) # ax.plot(drhodz_meanle400_plusstd,drhodz_depth,'g', linewidth = .5, zorder = 6) ax.fill_betweenx(drhodz_depth,drhodz_meang400_minusstd,drhodz_meang400_plusstd, facecolor = 'k', alpha=0.2,zorder = 1) ax.fill_betweenx(drhodz_depth,drhodz_meanle400_minusstd,drhodz_meanle400_plusstd, facecolor = 'g', alpha=0.2,zorder = 2) ax.axvline(x = 0.001, ymin=0, ymax=1, linewidth=1, color='b') ax.axvline(x = 0.002, ymin=0, ymax=1, linewidth=1, color='b') ax.axvline(x = 0.003, ymin=0, ymax=1, linewidth=1, color='b') ax.axvline(x = 0.004, ymin=0, ymax=1, linewidth=1, color='b') ax.axvline(x = 0.005, ymin=0, ymax=1, linewidth=1, color='b') #ax.fill_between(abs(drhodz_meanle400_2)-abs(drhodz_stdle400), abs(drhodz_meanle400_2)+abs(drhodz_stdle400), drhodz_depth,'g', alpha=0.2,zorder = 2) #ax2.set_ylabel('between y1 and 1') ax.set_xlabel(r'drho/dz (kgm$^{-4}$)') ax.set_ylabel(r'Depth (m)') ax.set_ylim([-300,0]) #ax.set_xlim([-.1,1.]) #axvline((1027.45-1026.), linewidth=3, color='b',linestyle=':', zorder = 2) # # Save and spc to storm plt.savefig('/home/ardavies/satdata/OSCAR/floatMKE_gandlemeandrhodz_profile_withstd.png') #pp = PdfPages('floatMKE_le400_profiles.pdf') from matplotlib.font_manager import FontProperties legendfont = FontProperties() legendfont.set_name('Computer Modern Roman') legendfont.set_size('x-small') rcParams['axes.labelsize'] = 18 rcParams['xtick.labelsize'] = 18 rcParams['ytick.labelsize'] = 18 rcParams['legend.fontsize'] = 14 # from matplotlib import rcParams rcParams['font.family'] = 'serif' rcParams['font.serif'] = ['Computer Modern Roman'] rcParams['text.usetex'] = True # fig = plt.figure() # ax4 = fig.add_axes([0.12,0.12,0.8,0.8]) from pylab import * contourdat= np.zeros([interplength,lenfloatMKEle400]) for ci in range(0,interplength): for cj in range(0,int(lenfloatMKEle400)): contourdat[ci,cj] = griddedprofiles_le400[ci,7,cj] # # Plot Set-up import math as ma from mpl_toolkits.axes_grid1 import make_axes_locatable fig = plt.figure() ax = fig.add_axes([0.15,0.1,0.68,0.85]) for cj in range(0,int(lenfloatMKEle400)): ax.plot(contourdat[:,cj],Ygrid,'0.75',linewidth = 1, zorder = 0) #ax.set_yticks([0,-100,-200,-300,-400, -500]) #ax.set_yticklabels([r'0', r'100',r'200',r'300',r'400',r'500']) # ax.set_xticks([40,60,80,100,120]) # ax.set_xticklabels([ r'40',r'60',r'80',r'100',r'120']) ax.set_xlabel(r'Density (kgm$^{-3}$)') ax.set_ylabel(r'Depth (m)') ax.set_ylim([-300,0]) ax.set_xlim([-.1,1.]) ax.plot(meanprofilele400,Ygrid,'k', linewidth = 3, zorder = 1) #axvline((1027.45-1026.), linewidth=3, color='b',linestyle=':', zorder = 2) # # Save and spc to storm plt.savefig('/home/ardavies/satdata/OSCAR/floatMKE_le400_profile.png') #pp.close() #os.system('scp /home/ardavies/satdata/OSCAR/pdfoutput/multiprofiles.pdf ardavies@storm.ceoe.udel.edu:/dev/ardavies/grlpaperplots/') from matplotlib.font_manager import FontProperties legendfont = FontProperties() legendfont.set_name('Computer Modern Roman') legendfont.set_size('x-small') rcParams['axes.labelsize'] = 18 rcParams['xtick.labelsize'] = 18 rcParams['ytick.labelsize'] = 18 rcParams['legend.fontsize'] = 14 # from matplotlib import rcParams rcParams['font.family'] = 'serif' rcParams['font.serif'] = ['Computer Modern Roman'] rcParams['text.usetex'] = True # fig = plt.figure() # ax4 = fig.add_axes([0.12,0.12,0.8,0.8]) from pylab import * contourdat= np.zeros([interplength,lenfloatMKEg400]) for ci in range(0,interplength): for cj in range(0,int(lenfloatMKEg400)): contourdat[ci,cj] = griddedprofiles_g400[ci,7,cj] # # Plot Set-up import math as ma from mpl_toolkits.axes_grid1 import make_axes_locatable fig = plt.figure() ax = fig.add_axes([0.15,0.1,0.68,0.85]) for cj in range(0,int(lenfloatMKEg400)): ax.plot(contourdat[:,cj],Ygrid,'0.75',linewidth = 1, zorder = 0) #ax.set_yticks([0,-100,-200,-300,-400, -500]) #ax.set_yticklabels([r'0', r'100',r'200',r'300',r'400',r'500']) # ax.set_xticks([40,60,80,100,120]) # ax.set_xticklabels([ r'40',r'60',r'80',r'100',r'120']) ax.set_xlabel(r'Density (kgm$^{-3}$)') ax.set_ylabel(r'Depth (m)') ax.set_ylim([-300,0]) ax.set_xlim([-.1,1.]) ax.plot(meanprofileg400,Ygrid,'k', linewidth = 3, zorder = 1) #axvline((1027.45-1026.), linewidth=3, color='b',linestyle=':', zorder = 2) # # Save and spc to storm plt.savefig('/home/ardavies/satdata/OSCAR/floatMKE_g400_profile.png') from matplotlib.font_manager import FontProperties legendfont = FontProperties() legendfont.set_name('Computer Modern Roman') legendfont.set_size('x-small') rcParams['axes.labelsize'] = 18 rcParams['xtick.labelsize'] = 18 rcParams['ytick.labelsize'] = 18 rcParams['legend.fontsize'] = 14 # from matplotlib import rcParams rcParams['font.family'] = 'serif' rcParams['font.serif'] = ['Computer Modern Roman'] rcParams['text.usetex'] = True # fig = plt.figure() # ax4 = fig.add_axes([0.12,0.12,0.8,0.8]) from pylab import * contourdat= np.zeros([interplength,lenfloatMKEg400]) for ci in range(0,interplength): for cj in range(0,int(lenfloatMKEg400)): contourdat[ci,cj] = griddedprofiles_g400[ci,7,cj] # # Plot Set-up import math as ma from mpl_toolkits.axes_grid1 import make_axes_locatable fig = plt.figure() ax = fig.add_axes([0.15,0.1,0.68,0.85]) ax.plot(meanprofileg400,Ygrid,'k', linewidth = 3, zorder = 1) ax.plot(meanprofilele400,Ygrid,'g', linewidth = 3, zorder = 1) ax.fill_betweenx(Ygrid,meanprofileg400_minusstd,meanprofileg400_plusstd, facecolor = 'k', alpha=0.2,zorder = 1) ax.fill_betweenx(Ygrid,meanprofilele400_minusstd,meanprofilele400_plusstd, facecolor = 'g', alpha=0.2,zorder = 2) ax.set_xlabel(r'Density (kgm$^{-3}$)') ax.set_ylabel(r'Depth (m)') ax.set_ylim([-300,0]) ax.set_xlim([-.1,1.]) #axvline((1027.45-1026.), linewidth=3, color='b',linestyle=':', zorder = 2) # # Save and spc to storm plt.savefig('/home/ardavies/satdata/OSCAR/floatMKE_gandlemean_profile.png') # # Do you want to run this code? gridWplt = 1 if gridWplt == 1: # # Opening pdf for printing at the end of loop pp1 = PdfPages('/home/ardavies/satdata/OSCAR/linear_density_chl_profiles_300m.pdf') # for i in range(0,arraylen): # # =========================================================== # Plotting Density Profiles and Filters (full profiiles) # =========================================================== # # Font Set-up import matplotlib.pyplot as plt import numpy as np import math as ma from mpl_toolkits.axes_grid1 import make_axes_locatable from matplotlib.font_manager import FontProperties legendfont = FontProperties() legendfont.set_name('Computer Modern Roman') legendfont.set_size('x-small') rcParams['axes.labelsize'] = 18 rcParams['xtick.labelsize'] = 18 rcParams['ytick.labelsize'] = 18 rcParams['legend.fontsize'] = 14 # from matplotlib import rcParams rcParams['font.family'] = 'serif' rcParams['font.serif'] = ['Computer Modern Roman'] rcParams['text.usetex'] = True # # fig = plt.figure() # ax = fig.add_subplot(1, 1, 1) # fig = plt.figure() host = fig.add_subplot(111) from pylab import * # # Sharing y axis par1 = host.twiny() # # Plotting p1 = host.plot(GridData[:,1,i],Ygrid,'k',linewidth = 3) # host.scatter(dataout[2,:],dataout[0,:],c='0.75') # host.plot(denavg3,depthavg3,'m',linewidth = 1) # host.plot(denavg5,depthavg5,'b',linewidth = 1) # host.plot(denavg7,depthavg7,'k',linewidth = 1) p2 = par1.plot(GridData[:,4,i],Ygrid,'g',linewidth = 1) # # Setting-up axis host.xaxis.set_major_locator(MaxNLocator(5)) host.set_ylim([-300,0]) par1.set_xlim([0,6]) # # # Ploting horizontal lines for derivative depths # host.axhline(xmin=0.0, xmax=0.16,y=thedepthofmax1[i],linewidth=1, color='0.75') # host.axhline(xmin=0.0, xmax=0.14,y=thedepthofmax3[i],linewidth=1, color='m') # host.axhline(xmin=0.0, xmax=0.11,y=thedepthofmax5[i],linewidth=1, color='b') # host.axhline(xmin=0.0, xmax=0.10,y=thedepthofmax7[i],linewidth=1, color='k') # # # host.axhline(xmin=0.84, xmax=1,y=thedepthofmax21[i],linewidth=1, color='0.75') # host.axhline(xmin=0.86, xmax=1,y=thedepthofmax23[i],linewidth=1, color='m') # host.axhline(xmin=0.88, xmax=1,y=thedepthofmax25[i],linewidth=1, color='b') # host.axhline(xmin=0.9, xmax=1,y=thedepthofmax27[i],linewidth=1, color='k') # # axis labels host.set_xlabel("Density") host.set_ylabel("Depth") par1.set_ylabel("Chlorophyll Fluorescence (micro g/l)") # # fig.suptitle("Date: " + str(floatdate[i])) fig.savefig(pp1, format = "pdf") pp1.close() #os.system('scp /home/ardavies/satdata/OSCAR/pdfoutput/linear_density_JanJunMission_full.pdf ardavies@storm.ceoe.udel.edu:/dev/ardavies/') # # Do you want to run this code? gridWplt = 1 if gridWplt == 1: # # Opening pdf for printing at the end of loop pp1 = PdfPages('/home/ardavies/satdata/OSCAR/linear_density_chl_profilesg400_300m.pdf') # for i in range(0,int(lenfloatMKEg400)): # # =========================================================== # Plotting Density Profiles and Filters (full profiiles) # =========================================================== # # Font Set-up import matplotlib.pyplot as plt import numpy as np import math as ma from mpl_toolkits.axes_grid1 import make_axes_locatable from matplotlib.font_manager import FontProperties legendfont = FontProperties() legendfont.set_name('Computer Modern Roman') legendfont.set_size('x-small') rcParams['axes.labelsize'] = 18 rcParams['xtick.labelsize'] = 18 rcParams['ytick.labelsize'] = 18 rcParams['legend.fontsize'] = 14 # from matplotlib import rcParams rcParams['font.family'] = 'serif' rcParams['font.serif'] = ['Computer Modern Roman'] rcParams['text.usetex'] = True # # fig = plt.figure() # ax = fig.add_subplot(1, 1, 1) # fig = plt.figure() host = fig.add_subplot(111) from pylab import * # # Sharing y axis par1 = host.twiny() # # Plotting p1 = host.plot(griddedprofiles_g400[:,7,i],Ygrid,'k',linewidth = 3) # host.scatter(dataout[2,:],dataout[0,:],c='0.75') # host.plot(denavg3,depthavg3,'m',linewidth = 1) # host.plot(denavg5,depthavg5,'b',linewidth = 1) # host.plot(denavg7,depthavg7,'k',linewidth = 1) p2 = par1.plot(griddedprofiles_g400[:,4,i],Ygrid,'g',linewidth = 1) # # Setting-up axis host.set_xlim([-.1,1]) host.set_ylim([-300,0]) par1.set_xlim([0,6]) # # # Ploting horizontal lines for derivative depths # host.axhline(xmin=0.0, xmax=0.16,y=thedepthofmax1[i],linewidth=1, color='0.75') # host.axhline(xmin=0.0, xmax=0.14,y=thedepthofmax3[i],linewidth=1, color='m') # host.axhline(xmin=0.0, xmax=0.11,y=thedepthofmax5[i],linewidth=1, color='b') # host.axhline(xmin=0.0, xmax=0.10,y=thedepthofmax7[i],linewidth=1, color='k') # # # host.axhline(xmin=0.84, xmax=1,y=thedepthofmax21[i],linewidth=1, color='0.75') # host.axhline(xmin=0.86, xmax=1,y=thedepthofmax23[i],linewidth=1, color='m') # host.axhline(xmin=0.88, xmax=1,y=thedepthofmax25[i],linewidth=1, color='b') # host.axhline(xmin=0.9, xmax=1,y=thedepthofmax27[i],linewidth=1, color='k') # # axis labels host.set_xlabel("Density") host.set_ylabel("Depth") par1.set_ylabel("Chlorophyll Fluorescence (micro g/l)") # # fig.suptitle("Date: " + str(floatdate_g400[i])) fig.savefig(pp1, format = "pdf") pp1.close() #os.system('scp /home/ardavies/satdata/OSCAR/pdfoutput/linear_density_JanJunMission_full.pdf ardavies@storm.ceoe.udel.edu:/dev/ardavies/') # # Do you want to run this code? gridWplt = 1 if gridWplt == 1: # # Opening pdf for printing at the end of loop pp1 = PdfPages('/home/ardavies/satdata/OSCAR/linear_density_chl_profilesle400_300m.pdf') # for i in range(0,int(lenfloatMKEle400)): # # =========================================================== # Plotting Density Profiles and Filters (full profiiles) # =========================================================== # # Font Set-up import matplotlib.pyplot as plt import numpy as np import math as ma from mpl_toolkits.axes_grid1 import make_axes_locatable from matplotlib.font_manager import FontProperties legendfont = FontProperties() legendfont.set_name('Computer Modern Roman') legendfont.set_size('x-small') rcParams['axes.labelsize'] = 18 rcParams['xtick.labelsize'] = 18 rcParams['ytick.labelsize'] = 18 rcParams['legend.fontsize'] = 14 # from matplotlib import rcParams rcParams['font.family'] = 'serif' rcParams['font.serif'] = ['Computer Modern Roman'] rcParams['text.usetex'] = True # # fig = plt.figure() # ax = fig.add_subplot(1, 1, 1) # fig = plt.figure() host = fig.add_subplot(111) from pylab import * # # Sharing y axis par1 = host.twiny() # # Plotting p1 = host.plot(griddedprofiles_le400[:,7,i],Ygrid,'k',linewidth = 3) # host.scatter(dataout[2,:],dataout[0,:],c='0.75') # host.plot(denavg3,depthavg3,'m',linewidth = 1) # host.plot(denavg5,depthavg5,'b',linewidth = 1) # host.plot(denavg7,depthavg7,'k',linewidth = 1) p2 = par1.plot(griddedprofiles_le400[:,4,i],Ygrid,'g',linewidth = 1) # # Setting-up axis host.set_xlim([-.1,1]) host.set_ylim([-300,0]) par1.set_xlim([0,6]) # # # Ploting horizontal lines for derivative depths # host.axhline(xmin=0.0, xmax=0.16,y=thedepthofmax1[i],linewidth=1, color='0.75') # host.axhline(xmin=0.0, xmax=0.14,y=thedepthofmax3[i],linewidth=1, color='m') # host.axhline(xmin=0.0, xmax=0.11,y=thedepthofmax5[i],linewidth=1, color='b') # host.axhline(xmin=0.0, xmax=0.10,y=thedepthofmax7[i],linewidth=1, color='k') # # # host.axhline(xmin=0.84, xmax=1,y=thedepthofmax21[i],linewidth=1, color='0.75') # host.axhline(xmin=0.86, xmax=1,y=thedepthofmax23[i],linewidth=1, color='m') # host.axhline(xmin=0.88, xmax=1,y=thedepthofmax25[i],linewidth=1, color='b') # host.axhline(xmin=0.9, xmax=1,y=thedepthofmax27[i],linewidth=1, color='k') # # axis labels host.set_xlabel("Density") host.set_ylabel("Depth") par1.set_ylabel("Chlorophyll Fluorescence (micro g/l)") # # fig.suptitle("Date: " + str(floatdate_le400[i])) fig.savefig(pp1, format = "pdf") pp1.close() #os.system('scp /home/ardavies/satdata/OSCAR/pdfoutput/linear_density_JanJunMission_full.pdf ardavies@storm.ceoe.udel.edu:/dev/ardavies/') # # Do you want to run this code? gridWplt = 1 if gridWplt == 1: # # Opening pdf for printing at the end of loop pp1 = PdfPages('/home/ardavies/satdata/OSCAR/linear_density_chl_drhodzprofilesg400_300m.pdf') # for i in range(0,int(lenfloatMKEg400)): # # =========================================================== # Plotting Density Profiles and Filters (full profiiles) # =========================================================== # # Font Set-up import matplotlib.pyplot as plt import numpy as np import math as ma from mpl_toolkits.axes_grid1 import make_axes_locatable from matplotlib.font_manager import FontProperties legendfont = FontProperties() legendfont.set_name('Computer Modern Roman') legendfont.set_size('x-small') rcParams['axes.labelsize'] = 18 rcParams['xtick.labelsize'] = 18 rcParams['ytick.labelsize'] = 18 rcParams['legend.fontsize'] = 14 # from matplotlib import rcParams rcParams['font.family'] = 'serif' rcParams['font.serif'] = ['Computer Modern Roman'] rcParams['text.usetex'] = True # # fig = plt.figure() # ax = fig.add_subplot(1, 1, 1) # fig = plt.figure() host = fig.add_subplot(111) from pylab import * # # Sharing y axis par1 = host.twiny() # # Plotting p1 = host.plot(abs(drhodz_g400[:,i]),drhodz_depth,'k',linewidth = 3) # host.scatter(dataout[2,:],dataout[0,:],c='0.75') # host.plot(denavg3,depthavg3,'m',linewidth = 1) # host.plot(denavg5,depthavg5,'b',linewidth = 1) # host.plot(denavg7,depthavg7,'k',linewidth = 1) p2 = par1.plot(griddedprofiles_g400[:,4,i],Ygrid,'g',linewidth = 1) # # Setting-up axis host.set_xlim([-.005,.02]) host.set_ylim([-300,0]) par1.set_xlim([0,6]) # # # Ploting horizontal lines for derivative depths # host.axhline(xmin=0.0, xmax=0.16,y=thedepthofmax1[i],linewidth=1, color='0.75') # host.axhline(xmin=0.0, xmax=0.14,y=thedepthofmax3[i],linewidth=1, color='m') # host.axhline(xmin=0.0, xmax=0.11,y=thedepthofmax5[i],linewidth=1, color='b') # host.axhline(xmin=0.0, xmax=0.10,y=thedepthofmax7[i],linewidth=1, color='k') # # # host.axhline(xmin=0.84, xmax=1,y=thedepthofmax21[i],linewidth=1, color='0.75') # host.axhline(xmin=0.86, xmax=1,y=thedepthofmax23[i],linewidth=1, color='m') # host.axhline(xmin=0.88, xmax=1,y=thedepthofmax25[i],linewidth=1, color='b') # host.axhline(xmin=0.9, xmax=1,y=thedepthofmax27[i],linewidth=1, color='k') host.axvline(x = 0.001, ymin=0, ymax=1, linewidth=1, color='b') host.axvline(x = 0.002, ymin=0, ymax=1, linewidth=1, color='b') host.axvline(x = 0.003, ymin=0, ymax=1, linewidth=1, color='b') host.axvline(x = 0.004, ymin=0, ymax=1, linewidth=1, color='b') host.axvline(x = 0.005, ymin=0, ymax=1, linewidth=1, color='b') # # axis labels host.set_xlabel("Density") host.set_ylabel("Depth") par1.set_ylabel("Chlorophyll Fluorescence (micro g/l)") # # fig.suptitle("Date: " + str(floatdate_g400[i])) fig.savefig(pp1, format = "pdf") pp1.close() #os.system('scp /home/ardavies/satdata/OSCAR/pdfoutput/linear_density_JanJunMission_full.pdf ardavies@storm.ceoe.udel.edu:/dev/ardavies/') # # Do you want to run this code? gridWplt = 1 if gridWplt == 1: # # Opening pdf for printing at the end of loop pp1 = PdfPages('/home/ardavies/satdata/OSCAR/linear_density_chl_drhodzprofilesle400_300m.pdf') # for i in range(0,int(lenfloatMKEle400)): # # =========================================================== # Plotting Density Profiles and Filters (full profiiles) # =========================================================== # # Font Set-up import matplotlib.pyplot as plt import numpy as np import math as ma from mpl_toolkits.axes_grid1 import make_axes_locatable from matplotlib.font_manager import FontProperties legendfont = FontProperties() legendfont.set_name('Computer Modern Roman') legendfont.set_size('x-small') rcParams['axes.labelsize'] = 18 rcParams['xtick.labelsize'] = 18 rcParams['ytick.labelsize'] = 18 rcParams['legend.fontsize'] = 14 # from matplotlib import rcParams rcParams['font.family'] = 'serif' rcParams['font.serif'] = ['Computer Modern Roman'] rcParams['text.usetex'] = True # # fig = plt.figure() # ax = fig.add_subplot(1, 1, 1) # fig = plt.figure() host = fig.add_subplot(111) from pylab import * # # Sharing y axis par1 = host.twiny() # # Plotting p1 = host.plot(abs(drhodz_le400[:,i]),drhodz_depth,'k',linewidth = 3) # host.scatter(dataout[2,:],dataout[0,:],c='0.75') # host.plot(denavg3,depthavg3,'m',linewidth = 1) # host.plot(denavg5,depthavg5,'b',linewidth = 1) # host.plot(denavg7,depthavg7,'k',linewidth = 1) p2 = par1.plot(griddedprofiles_le400[:,4,i],Ygrid,'g',linewidth = 1) # # Setting-up axis host.set_xlim([-.005,.02]) host.set_ylim([-300,0]) par1.set_xlim([0,6]) # # # Ploting horizontal lines for derivative depths # host.axhline(xmin=0.0, xmax=0.16,y=thedepthofmax1[i],linewidth=1, color='0.75') # host.axhline(xmin=0.0, xmax=0.14,y=thedepthofmax3[i],linewidth=1, color='m') # host.axhline(xmin=0.0, xmax=0.11,y=thedepthofmax5[i],linewidth=1, color='b') # host.axhline(xmin=0.0, xmax=0.10,y=thedepthofmax7[i],linewidth=1, color='k') # # # host.axhline(xmin=0.84, xmax=1,y=thedepthofmax21[i],linewidth=1, color='0.75') # host.axhline(xmin=0.86, xmax=1,y=thedepthofmax23[i],linewidth=1, color='m') host.axvline(x = 0.001, ymin=0, ymax=1, linewidth=1, color='b') host.axvline(x = 0.002, ymin=0, ymax=1, linewidth=1, color='b') host.axvline(x = 0.003, ymin=0, ymax=1, linewidth=1, color='b') host.axvline(x = 0.004, ymin=0, ymax=1, linewidth=1, color='b') host.axvline(x = 0.005, ymin=0, ymax=1, linewidth=1, color='b') # # axis labels host.set_xlabel("Density") host.set_ylabel("Depth") par1.set_ylabel("Chlorophyll Fluorescence (micro g/l)") # # fig.suptitle("Date: " + str(floatdate_le400[i])) fig.savefig(pp1, format = "pdf") pp1.close() #os.system('scp /home/ardavies/satdata/OSCAR/pdfoutput/linear_density_JanJunMission_full.pdf ardavies@storm.ceoe.udel.edu:/dev/ardavies/') # drhodz_g400 print lenfloatMKEg400 print lenfloatMKEle400 jhjkhjkh # # =========================================================== # # THORPE DISPLACEMENT # # =========================================================== # # Do you want to plot this coderunner = 1 if coderunner == 1: # # Array Initialization Density_newlist = np.zeros([interplength,arraylen]) Density_new = np.zeros([interplength,arraylen]) sortindex = np.zeros([interplength,arraylen]) z_new = np.zeros([interplength,arraylen]) thorpedisplacement = np.zeros([interplength,arraylen]) # # Re-sort densities so they are stratifies for i in range(0,arraylen): Density_newlist[:,i] = sorted(GridData[:,1,i]) # # Re-arrange so the ording is same (bottom to top) as grided profiles for i in range(0,arraylen): jj = interplength -1 for j in range(0,interplength): Density_new[j,i] = np.float(Density_newlist[jj,i]) jj = jj - 1 # # Calculate Thorpe Displacement Length for i in range(0,arraylen): for j in range(0,interplength): value, sortindex[j,i] = find_nearest(Density_new[:,i],GridData[j,1,i]) for i in range(0,arraylen): for j in range(0,interplength): z_new[j,i] = Ygrid[np.int(sortindex[j,i])] for i in range(0,arraylen): for j in range(0,interplength): thorpedisplacement[j,i] = z_new[j,i] - Ygrid[j] # # =========================================================== # # THORPE LENGTH SCALE PLOTTING # # =========================================================== # # Do you want to plot this? gridWplt = 2 if gridWplt == 1: # # =========================================================== # Plotting Individual Density and Thorpe Scale Profiles # =========================================================== # # Changing Directory to plotting output os.chdir('/home/ardavies/satdata/OSCAR/pdfoutput') pp = PdfPages('thorpedisplacementprofiles.pdf') for i in range(0,arraylen): # # Fonts from matplotlib.colors import LogNorm from mpl_toolkits.basemap import Basemap import matplotlib.pyplot as plt import numpy as np # fig = plt.figure() # # ax = fig.add_axes() # ax = fig.add_subplot(1, 1, 1) from matplotlib.font_manager import FontProperties legendfont = FontProperties() legendfont.set_name('Computer Modern Roman') legendfont.set_size('x-small') rcParams['axes.labelsize'] = 18 rcParams['xtick.labelsize'] = 18 rcParams['ytick.labelsize'] = 18 rcParams['legend.fontsize'] = 14 # from pylab import * from matplotlib import rcParams rcParams['font.family'] = 'serif' rcParams['font.serif'] = ['Computer Modern Roman'] rcParams['text.usetex'] = True fig = plt.figure() # ax = fig.add_axes() ax = fig.add_subplot(111) from matplotlib.font_manager import FontProperties legendfont = FontProperties() legendfont.set_name('Computer Modern Roman') legendfont.set_size('x-small') rcParams['axes.labelsize'] = 18 rcParams['xtick.labelsize'] = 18 rcParams['ytick.labelsize'] = 18 rcParams['legend.fontsize'] = 14 # # Plot set-up f = plt.figure() ax3 = f.add_subplot(121) ax = f.add_subplot(122) # # Plotting ax3.plot((GridData[:,1,i]-1026),Ygrid,'k',linewidth = 2) ax.plot(thorpedisplacement[:,i],Ygrid,'k',linewidth = 2) # # Axis and Label Set-up ax3.set_ylim([-2000,0]) ax3.set_xlim([0.5,2.0]) ax3.set_ylabel(r'Depth (m)') ax3.set_xlabel(r'Density + 1026 (kg m$^{-3}$)') # ax.set_ylim([-2000,0]) ax.set_xlim([-120,120]) ax.set_xlabel(r'Thorpe Displacement Length (m)') plt.setp(ax.get_yticklabels(), visible=False) # # Title f.suptitle('DOY ' + str(int(floatdate[i]))) # # Saving figure plt.savefig(pp, format = "pdf") pp.close() os.system('scp /home/ardavies/satdata/OSCAR/pdfoutput/thorpedisplacementprofiles.pdf ardavies@storm.ceoe.udel.edu:/dev/ardavies/grlpaperplots/') # # Do you want to plot this? gridWplt = 2 if gridWplt == 1: # # =========================================================== # Contour Plotting Thorpe Scale # =========================================================== # # Plot Set-up pp = PdfPages('thorpedisplacement.pdf') from matplotlib.colors import LogNorm from mpl_toolkits.basemap import Basemap import matplotlib.pyplot as plt import numpy as np import math as ma from mpl_toolkits.axes_grid1 import make_axes_locatable fig = plt.figure() ax = fig.add_axes([0.15,0.1,0.68,0.85]) from matplotlib.font_manager import FontProperties legendfont = FontProperties() legendfont.set_name('Computer Modern Roman') legendfont.set_size('x-small') rcParams['axes.labelsize'] = 18 rcParams['xtick.labelsize'] = 18 rcParams['ytick.labelsize'] = 18 rcParams['legend.fontsize'] = 14 # # Contour Plotting & Color Bar from pylab import * # # Contouring cs = plt.contourf(floatdate,Ygrid, thorpedisplacement, levels = np.linspace(-100,100,201)) plt.set_cmap('seismic') cs = plt.contourf(floatdate,Ygrid, thorpedisplacement, levels = np.linspace(-100,100,201)) plt.set_cmap('seismic') # # Building Colorbar divider = make_axes_locatable(ax) cax = divider.append_axes("right", size="5%", pad=0.05) # cbar = plt.colorbar(cs,cax=cax) # cbar.set_label(r"Thorpe Displacement Length (m)") cbar.set_ticks([-100,-50, 0, 50, 100]) cbar.set_ticklabels([r'-100.0',r'-50.0', r'0.0',r'50.0', r'100.0']) # ax.set_yticks([0,-100,-200,-300,-400, -500]) ax.set_yticklabels([r'0', r'100',r'200',r'300',r'400',r'500']) # ax.set_xticks([40,60,80,100,120]) ax.set_xticklabels([ r'40',r'60',r'80',r'100',r'120']) # ax.set_xlabel(r'Day of Year') ax.set_ylabel(r'Depth (m)') ax.set_ylim([-500,0]) ax.set_xlim([31,139]) # # Saving plt.savefig(pp, format = "pdf") pp.close() os.system('scp /home/ardavies/satdata/OSCAR/pdfoutput/thorpedisplacement.pdf ardavies@storm.ceoe.udel.edu:/dev/ardavies/grlpaperplots/') # # =========================================================== # # FINDING ISOPYCNAL DEPTHS & PLOTTING # # =========================================================== # # Do you want to run this code? gridWplt = 1 if gridWplt == 1: # # Initializing arrays depths thedepthofmax1 = np.zeros(arraylen) thedepthofmax3 = np.zeros(arraylen) thedepthofmax5 = np.zeros(arraylen) thedepthofmax7 = np.zeros(arraylen) thedepthofmax21 = np.zeros(arraylen) thedepthofmax23 = np.zeros(arraylen) thedepthofmax25 = np.zeros(arraylen) thedepthofmax27 = np.zeros(arraylen) # depthof1026= np.zeros(arraylen) depthof10272= np.zeros(arraylen) depthof10273= np.zeros(arraylen) depthof10274= np.zeros(arraylen) depthof102745= np.zeros(arraylen) depthof10275= np.zeros(arraylen) depthof102755= np.zeros(arraylen) # # for i in range(0,arraylen): # # =========================================================== # Get Data # =========================================================== # dataout, rows, cols = csvread('/data/orbprocess_mail/alex/Jan01_Jun04_2013/data/type/test_Mar2014/' + fullnames[i]) # # =========================================================== # Filtering Profiles: filter in wave num = avg in depth domain # =========================================================== # # Initializing arrays denavg3 = np.zeros(rows-2) depthavg3 = np.zeros(rows-2) denavg5 = np.zeros(rows-4) depthavg5 = np.zeros(rows-4) denavg7 = np.zeros(rows-6) depthavg7 = np.zeros(rows-6) # # Filtering for ii in range(0,rows-2): denavg3[ii] = (dataout[2,ii] + dataout[2,ii + 1] + dataout[2,ii + 2])/3 depthavg3[ii] = (dataout[0,ii] + dataout[0,ii + 1] + dataout[0,ii + 2])/3 for iii in range(0,rows-4): denavg5[iii] = (dataout[2,iii] + dataout[2,iii + 1] + dataout[2,iii + 2]+ dataout[2,iii + 3]+ dataout[2,iii + 4])/5 depthavg5[iii] = (dataout[0,iii] + dataout[0,iii + 1] + dataout[0,iii + 2]+ dataout[0,iii + 3]+ dataout[0,iii + 4])/5 for iiii in range(0,rows-6): denavg7[iiii] = (dataout[2,iiii] + dataout[2,iiii + 1] + dataout[2,iiii + 2]+ dataout[2,iiii + 3] + dataout[2,iiii + 4] + dataout[2,iiii + 5] + dataout[2,iiii + 6])/7 depthavg7[iiii] = (dataout[0,iiii] + dataout[0,iiii + 1] + dataout[0,iiii + 2]+ dataout[0,iiii + 3] + dataout[0,iiii + 4] + dataout[0,iiii + 5] + dataout[0,iiii + 6])/7 # # =========================================================== # d(rho)/dz of Filtered Profiles # =========================================================== # # Initialization for first derivative drhodz = np.zeros(rows-1) drhodz_depth = np.zeros(rows-1) drhodz3 = np.zeros(rows-3) drhodz3_depth = np.zeros(rows-3) drhodz5 = np.zeros(rows-5) drhodz5_depth = np.zeros(rows-5) drhodz7 = np.zeros(rows-7) drhodz7_depth = np.zeros(rows-7) # # Initialization for second derivative d2rhodz = np.zeros(rows-2) d2rhodz_depth = np.zeros(rows-2) d2rhodz3 = np.zeros(rows-4) d2rhodz3_depth = np.zeros(rows-4) d2rhodz5 = np.zeros(rows-6) d2rhodz5_depth = np.zeros(rows-6) d2rhodz7 = np.zeros(rows-8) d2rhodz7_depth = np.zeros(rows-8) # # Derivative for ii in range(0,rows-1): drhodz[ii] = (dataout[2,ii+1]-dataout[2,ii])/(dataout[0,ii+1]-dataout[0,ii]) drhodz_depth[ii] = (dataout[0,ii]+dataout[0,ii+1])/2 for ii in range(0,rows-3): drhodz3[ii] = (denavg3[ii+1]-denavg3[ii])/(depthavg3[ii+1] -depthavg3[ii]) drhodz3_depth[ii] = (depthavg3[ii] + depthavg3[ii+1])/2 for ii in range(0,rows-5): drhodz5[ii] = (denavg5[ii+1]-denavg5[ii])/(depthavg5[ii+1] -depthavg5[ii]) drhodz5_depth[ii] = (depthavg5[ii] + depthavg5[ii+1])/2 for ii in range(0,rows-7): drhodz7[ii] = (denavg7[ii+1]-denavg7[ii])/(depthavg7[ii+1] -depthavg7[ii]) drhodz7_depth[ii] = (depthavg7[ii] + depthavg7[ii+1])/2 # # Find the max d(rho)/dz (right hand corrdinates) mx1 = np.amin(drhodz) mx3 = np.amin(drhodz3) mx5 = np.amin(drhodz5) mx7 = np.amin(drhodz7) # # Find index of max d(rho)/dz idx1 = (np.abs(drhodz-mx1)).argmin() idx3 = (np.abs(drhodz3-mx3)).argmin() idx5 = (np.abs(drhodz5-mx5)).argmin() idx7 = (np.abs(drhodz7-mx7)).argmin() # # Depth of max d(rho)/dz thedepthofmax1[i] = drhodz_depth[idx1] thedepthofmax3[i] = drhodz3_depth[idx3] thedepthofmax5[i] = drhodz5_depth[idx5] thedepthofmax7[i] = drhodz7_depth[idx7] # # Second Derivative for ii in range(0,rows-2): d2rhodz[ii] = (drhodz[ii+1]-drhodz[ii])/(drhodz_depth[ii+1]-drhodz_depth[ii]) d2rhodz_depth[ii] = (drhodz_depth[ii+1]+drhodz_depth[ii])/2 for ii in range(0,rows-4): d2rhodz3[ii] = (drhodz3[ii+1]-drhodz3[ii])/(drhodz3_depth[ii+1] -drhodz3_depth[ii]) d2rhodz3_depth[ii] = (drhodz3_depth[ii] + drhodz3_depth[ii+1])/2 for ii in range(0,rows-6): d2rhodz5[ii] = (drhodz5[ii+1]-drhodz5[ii])/(drhodz5_depth[ii+1] -drhodz5_depth[ii]) d2rhodz5_depth[ii] = (drhodz5_depth[ii] + drhodz5_depth[ii+1])/2 for ii in range(0,rows-8): d2rhodz7[ii] = (drhodz7[ii+1]-drhodz7[ii])/(drhodz7_depth[ii+1] -drhodz7_depth[ii]) d2rhodz7_depth[ii] = (drhodz7_depth[ii] + drhodz7_depth[ii+1])/2 # # Find the index where d2(rho)/dz2 = 0 (right hand corrdinates) idx21 = (np.abs(d2rhodz-0.0)).argmin() idx23 = (np.abs(d2rhodz3-0.0)).argmin() idx25 = (np.abs(d2rhodz5-0.0)).argmin() idx27 = (np.abs(d2rhodz7-0.0)).argmin() # # Depth of inflection point thedepthofmax21[i] = d2rhodz_depth[idx1] thedepthofmax23[i] = d2rhodz3_depth[idx3] thedepthofmax25[i] = d2rhodz5_depth[idx5] thedepthofmax27[i] = d2rhodz7_depth[idx7] # # =========================================================== # Filtering Isopychnal Depths # =========================================================== # # Finding each isopycnals in each profile value1026, index1026 = find_nearest(dataout[2,:],1026) value10272, index10272 = find_nearest(dataout[2,:],1027.2) value10273, index10273 = find_nearest(dataout[2,:],1027.3) value10274, index10274 = find_nearest(dataout[2,:],1027.4) value102745, index102745 = find_nearest(dataout[2,:],1027.45) value10275, index10275 = find_nearest(dataout[2,:],1027.5) value102755, index102755 = find_nearest(dataout[2,:],1027.55) # # Finding depths of isopycnals in each profile depthof1026[i] = dataout[0,index1026] depthof10272[i] = dataout[0,index10272] depthof10273[i] = dataout[0,index10273] depthof10274[i] = dataout[0,index10274] depthof102745[i] = dataout[0,index102745] depthof10275[i] = dataout[0,index10275] depthof102755[i] = dataout[0,index102755]