import netCDF4 as nc import numpy as np import matplotlib.pyplot as plt from scipy.optimize import curve_fit import pandas as pd import csv %matplotlib inline path = '/data/nsoontie/MEOPAR/SalishSea/results/stratification/no_strat/2d_grid_T.nc' gridT=nc.Dataset(path) temp=gridT.variables['votemper'] sal=gridT.variables['vosaline'] sal.shape fig,axs=plt.subplots(1,2,figsize=(10,10)) ax=axs[0] mesh=ax.pcolormesh(sal[4,0,:,:]) fig.colorbar(mesh,ax=ax) ax.set_title('Salinity') ax=axs[1] mesh=ax.pcolormesh(temp[4,0,:,:]) fig.colorbar(mesh,ax=ax) ax.set_title('Temperature') dep_d = gridT.variables['deptht'] # Call thalweg lines = np.loadtxt('/data/nsoontie/MEOPAR/tools/bathymetry/thalweg_working.txt', delimiter=" ", unpack=False) lines = lines.astype(int) ds=np.arange(0,lines.shape[0],1); XX,ZZ = np.meshgrid(ds,-dep_d[:]) # Salinity along thalweg salP=sal[-1,:,lines[:,0],lines[:,1]] salP= np.ma.masked_values(salP,0) tempP=temp[-1,:,lines[:,0],lines[:,1]] tempP= np.ma.masked_values(tempP,0) fig,axs=plt.subplots(2,1,figsize=(20,5)) mesh=axs[0].pcolormesh(XX,ZZ,salP) axs[0].set_title('salinity') fig.colorbar(mesh,ax=axs[0]) mesh=axs[1].pcolormesh(XX,ZZ,tempP) axs[1].set_title('temperature') fig.colorbar(mesh,ax=axs[1]) # pathname for data - all of the tide runs are stored in this directory path = '/data/nsoontie/MEOPAR/SalishSea/results/stratification/no_strat/' #path = '../../myResults/' #joining the two string together name = path print name #observations filename = '/data/nsoontie/MEOPAR/analysis/compare_tides/obs_tidal_wlev_const_all.csv' harm_obs = pd.read_csv(filename,sep=';',header=0) harm_obs = harm_obs.rename(columns={'Site': 'site', 'Lat': 'lat', 'Lon': 'lon', 'M2 amp': 'M2_amp', 'M2 phase (deg UT)': 'M2_pha', 'K1 amp': 'K1_amp', 'K1 phase (deg UT)': 'K1_pha'}) print harm_obs stations = ['PortRenfrew','SheringhamPoint','PedderBay', 'Esquimalt', 'Victoria','CloverPoint','FinnertyCove', 'FulfordHarbour', 'TumboChannel','PatosIsland','WhalerBay', 'Tsawwassen', 'Sandheads', 'PointGrey','PointAtkinson','GibsonsLanding', 'WinchelseaIs', 'HalfmoonBay','IrvinesLanding','PowellRiver', 'LittleRiver', 'Lund', 'TwinIslets','CampbellRiver','MaudeIslandE', 'NympheCove', 'SeymourNarrows','BrownBay','ChathamPoint','KelseyBay','YorkeIsland'] numsta=len(stations) #again with spaces because the text file likes that stations_obs = ['Port Renfrew','Sheringham Point','Pedder Bay', 'Esquimalt', 'Victoria','Clover Point','Finnerty Cove', 'Fulford Harbour', 'Tumbo Channel','Patos Island','Whaler Bay', 'Tsawwassen', 'Sandheads', 'Point Grey','Point Atkinson','Gibsons Landing', 'Winchelsea', 'Halfmoon Bay','Irvines Landing','Powell River', 'Little River', 'Lund', 'Twin Islets','Campbell River','Maude Island E', 'Nymphe Cove', 'Seymour Narrows','Brown Bay','Chatham Point','Kelsey Bay','Yorke Island'] stn_lon=np.zeros(numsta); stn_lat=np.zeros(numsta) #ax.pcolormesh(X,Y,bathy,cmap='winter_r') for stn in range(numsta): location = stations_obs[stn] lon=-harm_obs.lon[harm_obs.site==location] lat=harm_obs.lat[harm_obs.site==location] stn_lon[stn]=lon stn_lat[stn]=lat #constants and fitting # M2 M2freq = 28.984106 # degrees per hour M2freq = M2freq*np.pi/180. # radians per hour #K1 K1freq = 15.041069*np.pi/180. #O1 O1freq = 13.943036*np.pi/180. #S2 S2freq = 30.000002*np.pi/180. #P1 P1freq = 14.958932*np.pi/180. #N2 N2freq = 28.439730*np.pi/180. #Q1 Q1freq = 13.398661*np.pi/180. #K2 K2freq = 30.082138*np.pi/180. # initial phase calculation # our start is currently April 21 2003 # data for phase output from bdytides.F90; found in ocean.output K1ft = 1.065505 K1uvt = 111.481741 M2ft = 0.982328 M2uvt = 250.506179 O1ft = 1.105495 O1uvt = 142.040782 S2ft = 1.0 S2uvt = 0.0 P1ft = 1.0 P1uvt =241.335269 N2ft = 0.982328 N2uvt =205.684028 Q1ft = 1.105495 Q1uvt =97.218631 K2ft = 1.159036 K2uvt = 42.361669 # function for fitting 8 frequencies def octuple(x, M2amp, M2pha, K1amp, K1pha, O1amp, O1pha, S2amp, S2pha, P1amp, P1pha, N2amp, N2pha, Q1amp, Q1pha, K2amp, K2pha): return (M2amp*np.cos(M2freq*x-M2pha*np.pi/180.)+ K1amp*np.cos(K1freq*x-K1pha*np.pi/180.)+ O1amp*np.cos(O1freq*x-O1pha*np.pi/180.)+ S2amp*np.cos(S2freq*x-S2pha*np.pi/180.)+ P1amp*np.cos(P1freq*x-P1pha*np.pi/180.)+ N2amp*np.cos(N2freq*x-N2pha*np.pi/180.)+ Q1amp*np.cos(Q1freq*x-Q1pha*np.pi/180.)+ K2amp*np.cos(K2freq*x-K2pha*np.pi/180.)) fig, ax = plt.subplots(1,1,figsize=(12,5)) for stn in (0,4,14,23): print stations[stn] fT1 = nc.Dataset(name+stations[stn]+'.nc','r') time = fT1.variables["time_counter"][:]/3600. # want hours not seconds ssh = fT1.variables["sossheig"][:,0,0] ax.plot(time,ssh) #allocate space for our arrays M2_amp=[]; M2_pha=[]; K1_amp=[]; K1_pha=[] M2_amp_obs=np.zeros(numsta); M2_pha_obs=np.zeros(numsta) K1_amp_obs=np.zeros(numsta); K1_pha_obs=np.zeros(numsta) ts = 240 te = ssh.shape[0] for stn in range(numsta): fT1 = nc.Dataset(name+stations[stn]+'.nc','r') time = fT1.variables["time_counter"][:]/3600. # want hours not seconds ssh = fT1.variables["sossheig"][:,0,0] fitted, cov = curve_fit(octuple,time[ts:te],ssh[ts:te]) if fitted[0] < 0: fitted[0] = -fitted[0] fitted[1] = fitted[1]+180 if fitted[2] < 0: fitted[2] = -fitted[2] fitted[3] = fitted[3]+180 M2_amp.append(fitted[0]*M2ft) pha = fitted[1]+M2uvt if pha > 360: pha=pha-360 elif pha < 0: pha = pha+360 if stn == 6: print pha M2_pha.append(pha) K1_amp.append(fitted[2]*K1ft) pha= fitted[3]+K1uvt if pha > 360: pha=pha-360 K1_pha.append(pha) #now the observations location=stations_obs[stn] M2_amp_obs[stn]=harm_obs.M2_amp[harm_obs.site==location]/100 M2_pha_obs[stn]=harm_obs.M2_pha[harm_obs.site==location] K1_amp_obs[stn]=harm_obs.K1_amp[harm_obs.site==location]/100 K1_pha_obs[stn]=harm_obs.K1_pha[harm_obs.site==location] def mean(diff): return np.mean(abs(diff)) def rms(diff): return np.sqrt(np.mean(diff**2)) def complex_diff(Ao,go,Am,gm): #calculates complex differences between observations and model #Ao, go - amplitude and phase from observations #Am, gm - amplitude and phase from model D = np.sqrt((Ao*np.cos(np.pi*go/180)-Am*np.cos(np.pi*gm/180))**2 + (Ao*np.sin(np.pi*go/180)-Am*np.sin(np.pi*gm/180))**2) return D #R R_M2 = M2_amp/M2_amp_obs R_K1 = K1_amp/K1_amp_obs #delta phi (adjust so between -180, 180) Dphi_M2=M2_pha-M2_pha_obs; Dphi_M2 = Dphi_M2 -360*(Dphi_M2>180) + 360*(Dphi_M2<-180) Dphi_K1=K1_pha-K1_pha_obs Dphi_K1 = Dphi_K1 -360*(Dphi_K1>180) + 360*(Dphi_K1<-180) #Complex differences D_M2= complex_diff(np.array(M2_amp_obs),np.array(M2_pha_obs), np.array(M2_amp),np.array(M2_pha)) D_K1= complex_diff(np.array(K1_amp_obs),np.array(K1_pha_obs), np.array(K1_amp),np.array(K1_pha)) outfile = 'tides_nostrat.csv' split1=8; split2=20 with open(outfile, 'wb') as csvfile: writer = csv.writer(csvfile, delimiter=',') writer.writerow([ 'Station Number', 'Station Name', 'Longitude', 'Latitude', 'R (M2)', 'Delta phi (M2)', 'D (M2)', 'R (K1)', 'Delta phi (K1)', 'D (K1)' ]) for stn in range(numsta): location = stations_obs[stn] writer.writerow( [str(stn +1),stations_obs[stn], str(stn_lon[stn]), str(stn_lat[stn]), R_M2[stn], Dphi_M2[stn], D_M2[stn], R_K1[stn], Dphi_K1[stn], D_K1[stn]]) #write averages and rms writer.writerow(['','','','Mean Difference', mean(M2_amp-M2_amp_obs),mean(Dphi_M2),mean(D_M2), mean(K1_amp-K1_amp_obs),mean(Dphi_K1),mean(D_K1)]) writer.writerow(['','','','RMS Difference', rms(M2_amp-M2_amp_obs),rms(Dphi_M2),rms(D_M2), rms(K1_amp-K1_amp_obs),rms(Dphi_K1),rms(D_K1)]) #without the north writer.writerow(['','','','Mean Difference no North no PR', mean(M2_amp[1:split2]-M2_amp_obs[1:split2]),mean(Dphi_M2[1:split2]),mean(D_M2[1:split2]), mean(K1_amp[1:split2]-K1_amp_obs[1:split2]),mean(Dphi_K1[1:split2]),mean(D_K1[1:split2])]) writer.writerow(['','','','RMS Difference no North no PR', rms(M2_amp[1:split2]-M2_amp_obs[1:split2]),rms(Dphi_M2[1:split2]),rms(D_M2[1:split2]), rms(K1_amp[1:split2]-K1_amp_obs[1:split2]),rms(Dphi_K1[1:split2]),rms(D_K1[1:split2])])