%matplotlib inline from matplotlib import pylab import matplotlib.pyplot as plt import matplotlib as mpl import netCDF4 as NC import numpy as np import scipy.interpolate as sp import math def combine_data(data_list): #combines ouput files so that they are easier to analyze #data_list is a dict object that contains the netcdf handles for the files of importance. us={}; vs={}; lats={}; lons={}; sals={}; tmps={}; sshs={}; for k in data_list: net=data_list.get(k) us[k]=net.variables['vozocrtx'] vs[k]=net.variables['vomecrty'] lats[k]=net.variables['nav_lat'] lons[k]=net.variables['nav_lon'] tmps[k]=net.variables['votemper'] sals[k]=net.variables['vosaline'] sshs[k]=net.variables['sossheig'] return us, vs, lats, lons, tmps, sals, sshs #load up the bathymetry. grid = NC.Dataset('/data/nsoontie/MEOPAR/NEMO-forcing/grid/bathy_meter_SalishSea2.nc','r') bathy = grid.variables['Bathymetry'][:,:] X = grid.variables['nav_lon'][:,:] Y = grid.variables['nav_lat'][:,:] #load in the datas. Trying something new with dict objects... r1 = 'r3942-4510_long'; runname1 = '/data/nsoontie/MEOPAR/SalishSea/results/merge-tests/' + r1 print runname1 run1_thal={} for k in range(1,7): string = 'thal' + str(k) run1_thal[string] = NC.Dataset(runname1 +'/1h_Thalweg' +str(k)+ '.nc','r'); #second dataset r2 = 'r3942-4510Jasper_long' runname2 = '/data/nsoontie/MEOPAR/SalishSea/results/merge-tests/' +r2 print runname2 run2_thal={} for k in range(1,7): string = 'thal' + str(k) run2_thal[string] = NC.Dataset(runname2 +'/1h_Thalweg' +str(k)+ '.nc','r'); [us1, vs1, lats1, lons1, tmps1, sals1, sshs1] = combine_data(run1_thal) [us2, vs2, lats2, lons2, tmps2, sals2, sshs2] = combine_data(run2_thal) #finally grab a times placer and depths net=run1_thal.get('thal1') times = net.variables['time_counter'] depths = net.variables['deptht'] #Plot bathy ax=[-126.3,-122.1,47,51] istorms=np.array([329, 196, 214, 124]); jstorms=np.array([469, 299, 352, 748]); iplume = 290 jplume = 440 fig =plt.figure(figsize=(6,8)) pylab.pcolor(X,Y,bathy) pylab.colorbar() for key in lats1: pylab.plot(lons1[key][0,0], lats1[key][0,0],'o',label=key) pylab.plot(X[jstorms-1,istorms-1],Y[jstorms-1,istorms-1], 'ko') pylab.plot(X[jplume-1,iplume-1],Y[jplume-1,iplume-1], 'ko') pylab.axis(ax) pylab.legend(loc=0) plt.figure(figsize=(20,8)) panel=1 ax=[0,7,-0.002,0.002]; for key in run1_thal: pylab.subplot(3,2,panel) ssh_diff = sshs1[key][:,0,0] - sshs2[key][:,0,0]; #pylab.plot(times[:]/86400,sshs1[key][:,0,0],label=r1) #pylab.plot(times[:]/86400,sshs2[key][:,0,0],label=r2) pylab.plot(times[:]/86400,ssh_diff,label=r1 +' - ' +r2) pylab.title(key +' (' + str(lats1[key][0,0]) +', '+str(lons1[key][0,0]) +')') pylab.ylabel('ssh (m)') pylab.axis(ax) panel = panel+1 pylab.legend(loc=2) plt.figure(figsize=(20,8)) panel=1 ax=[0,7,-0.05,0.05]; depthlevel=0; print 'Plotting diffefrence in salinity: ' +r1 +' - ' +r2 for key in run1_thal: pylab.subplot(3,2,panel) sal_diff = sals1[key][:,depthlevel,0,0] - sals2[key][:,depthlevel,0,0]; #pylab.plot(times[:]/86400,sals1[key][:,depthlevel,0,0],label=r1) #pylab.plot(times[:]/86400,sals2[key][:,depthlevel,0,0],label=r2) pylab.plot(times[:]/86400,sal_diff,label = 'surface') pylab.title(key +' (' + str(lats1[key][0,0]) +', '+str(lons1[key][0,0]) +')') pylab.ylabel('Salinity') pylab.axis(ax) panel = panel+1 # do it again for 25m depth panel = 1 depthlevel=25 for key in run1_thal: pylab.subplot(3,2,panel) sals1P=np.zeros(times.shape[0]); sals2P=np.zeros(times.shape[0]); # interpolate to 25m first for t in range(0,times.shape[0]): f = sp.interp1d(depths[:],sals1[key][t,:,0,0]); sals1P[t] = f(depthlevel) f = sp.interp1d(depths[:],sals2[key][t,:,0,0]); sals2P[t] = f(depthlevel) sal_diff = sals1P-sals2P #pylab.plot(times[:]/86400,sals1P,'b*',label=depthlevel) #pylab.plot(times[:]/86400,sals2P,'g*') pylab.plot(times[:]/86400,sal_diff, label = depthlevel) panel = panel+1 # do it again for 100m depth panel = 1 depthlevel=100 for key in run1_thal: pylab.subplot(3,2,panel) sals1P=np.zeros(times.shape[0]); sals2P=np.zeros(times.shape[0]); # interpolate to 25m first for t in range(0,times.shape[0]): f = sp.interp1d(depths[:],sals1[key][t,:,0,0]); sals1P[t] = f(depthlevel) f = sp.interp1d(depths[:],sals2[key][t,:,0,0]); sals2P[t] = f(depthlevel) sal_diff = sals1P-sals2P #pylab.plot(times[:]/86400,sals1P,'bs',label=depthlevel) #pylab.plot(times[:]/86400,sals2P,'gs') pylab.plot(times[:]/86400,sal_diff,label = depthlevel) panel = panel+1 pylab.legend(loc=2) plt.figure(figsize=(20,8)) panel=1 ax=[0,7,-.005,.005]; depthlevel=0; for key in run1_thal: pylab.subplot(3,2,panel) u_diff = us1[key][:,depthlevel,0,0] - us2[key][:,depthlevel,0,0]; #pylab.plot(times[:]/86400,us1[key][:,depthlevel,0,0],label=r1) #pylab.plot(times[:]/86400,us2[key][:,depthlevel,0,0],label=r2) pylab.plot(times[:]/86400,u_diff,label='surface') pylab.title(key +' (' + str(lats1[key][0,0]) +', '+str(lons1[key][0,0]) +')') pylab.ylabel('u (m/s)') pylab.axis(ax) panel = panel+1 # do it again for 25m depth panel = 1 depthlevel=25 for key in run1_thal: pylab.subplot(3,2,panel) us1P=np.zeros(times.shape[0]); us2P=np.zeros(times.shape[0]); # interpolate to 25m first for t in range(0,times.shape[0]): f = sp.interp1d(depths[:],us1[key][t,:,0,0]); us1P[t] = f(depthlevel) f = sp.interp1d(depths[:],us2[key][t,:,0,0]); us2P[t] = f(depthlevel) u_diff = us1P-us2P #pylab.plot(times[:]/86400,us1P,'b*',label=depthlevel) #pylab.plot(times[:]/86400,us2P,'g*') pylab.plot(times[:]/86400,u_diff,label=depthlevel) panel = panel+1 # do it again for 100m depth panel = 1 depthlevel=100 for key in run1_thal: pylab.subplot(3,2,panel) us1P=np.zeros(times.shape[0]); us2P=np.zeros(times.shape[0]); # interpolate to 25m first for t in range(0,times.shape[0]): f = sp.interp1d(depths[:],us1[key][t,:,0,0]); us1P[t] = f(depthlevel) f = sp.interp1d(depths[:],us2[key][t,:,0,0]); us2P[t] = f(depthlevel) u_diff = us1P-us2P #pylab.plot(times[:]/86400,us1P,'bs',label=depthlevel) #pylab.plot(times[:]/86400,us2P,'gs') pylab.plot(times[:]/86400,u_diff,label=depthlevel) panel = panel+1 pylab.legend(loc=2) plt.figure(figsize=(20,8)) panel=1 ax=[0,7,-0.005,0.005]; depthlevel=0; for key in run1_thal: pylab.subplot(3,2,panel) v_diff = vs1[key][:,depthlevel,0,0] - vs2[key][:,depthlevel,0,0]; #pylab.plot(times[:]/86400,vs1[key][:,depthlevel,0,0],label=r1) #pylab.plot(times[:]/86400,vs2[key][:,depthlevel,0,0],label=r2) pylab.plot(times[:]/86400,v_diff,label='surface') pylab.title(key +' (' + str(lats1[key][0,0]) +', '+str(lons1[key][0,0]) +')') pylab.ylabel('v (m/s)') pylab.axis(ax) panel = panel+1 # do it again for 25m depth panel = 1 depthlevel=25 for key in run1_thal: pylab.subplot(3,2,panel) vs1P=np.zeros(times.shape[0]); vs2P=np.zeros(times.shape[0]); # interpolate to 25m first for t in range(0,times.shape[0]): f = sp.interp1d(depths[:],vs1[key][t,:,0,0]); vs1P[t] = f(depthlevel) f = sp.interp1d(depths[:],vs2[key][t,:,0,0]); vs2P[t] = f(depthlevel) v_diff = vs1P-vs2P #pylab.plot(times[:]/86400,vs1P,'b*',label=depthlevel) #pylab.plot(times[:]/86400,vs2P,'g*') pylab.plot(times[:]/86400,v_diff,label=depthlevel) panel = panel+1 # do it again for 100m depth panel = 1 depthlevel=100 for key in run1_thal: pylab.subplot(3,2,panel) vs1P=np.zeros(times.shape[0]); vs2P=np.zeros(times.shape[0]); # interpolate to 100m first for t in range(0,times.shape[0]): f = sp.interp1d(depths[:],vs1[key][t,:,0,0]); vs1P[t] = f(depthlevel) f = sp.interp1d(depths[:],vs2[key][t,:,0,0]); vs2P[t] = f(depthlevel) v_diff = vs1P-vs2P #pylab.plot(times[:]/86400,vs1P,'bs',label=depthlevel) #pylab.plot(times[:]/86400,vs2P,'gs') pylab.plot(times[:]/86400,v_diff,label=depthlevel) panel = panel+1 pylab.legend(loc=2) plt.figure(figsize=(20,8)) panel=1 ax=[0,7,-.01,0.01]; depthlevel=0; for key in run1_thal: pylab.subplot(3,2,panel) T_diff = tmps1[key][:,depthlevel,0,0] - tmps2[key][:,depthlevel,0,0]; pylab.plot(times[:]/86400,T_diff,label='surface') pylab.title(key +' (' + str(lats1[key][0,0]) +', '+str(lons1[key][0,0]) +')') pylab.ylabel('Temperature') pylab.axis(ax) panel = panel+1 # do it again for 25m depth panel = 1 depthlevel=25 for key in run1_thal: pylab.subplot(3,2,panel) tmps1P=np.zeros(times.shape[0]); tmps2P=np.zeros(times.shape[0]); # interpolate to 25m first for t in range(0,times.shape[0]): f = sp.interp1d(depths[:],tmps1[key][t,:,0,0]); tmps1P[t] = f(depthlevel) f = sp.interp1d(depths[:],tmps2[key][t,:,0,0]); tmps2P[t] = f(depthlevel) T_diff = tmps1P-tmps2P pylab.plot(times[:]/86400,T_diff,label=depthlevel) panel = panel+1 # do it again for 100m depth panel = 1 depthlevel=100 for key in run1_thal: pylab.subplot(3,2,panel) tmps1P=np.zeros(times.shape[0]); tmps2P=np.zeros(times.shape[0]); # interpolate to 25m first for t in range(0,times.shape[0]): f = sp.interp1d(depths[:],tmps1[key][t,:,0,0]); tmps1P[t] = f(depthlevel) f = sp.interp1d(depths[:],tmps2[key][t,:,0,0]); tmps2P[t] = f(depthlevel) T_diff = tmps1P-tmps2P pylab.plot(times[:]/86400,T_diff,label=depthlevel) panel = panel+1 pylab.legend(loc=2) stations = {'PointAtkinson': 1, 'Victoria': 1, 'PatriciaBay': 1, 'CampbellRiver': 1, 'Plume': 1} run1_stations={} run2_stations ={} for key in stations: string = runname1 + '/1h_' + key + '.nc' run1_stations[key] = NC.Dataset(string,'r'); string = runname2 + '/1h_' + key + '.nc' run2_stations[key] = NC.Dataset(string,'r'); [us1, vs1, lats1, lons1, tmps1, sals1, sshs1] = combine_data(run1_stations) [us2, vs2, lats2, lons2, tmps2, sals2, sshs2] = combine_data(run2_stations) print 'plotting ' + r1 +' - ' + r2 fig, ((ax_ssh, ax_sal), (ax_u, ax_v)) = plt.subplots(2, 2, figsize=(18, 10)) for locn in stations: ssh_diff=sshs1[locn][:,0,0] -sshs2[locn][:,0,0] u_diff=us1[locn][:,0,0,0] -us2[locn][:,0,0,0] v_diff=vs1[locn][:,0,0,0] -vs2[locn][:,0,0,0] sals_diff=sals1[locn][:,0,0,0] -sals2[locn][:,0,0,0] ax_ssh.plot(times[:]/86400, ssh_diff, marker='+', label=locn) ax_sal.plot(times[:]/86400,sals_diff, marker='+', label=locn) ax_u.plot(times[:]/86400,u_diff, marker='+', label=locn) ax_v.plot(times[:]/86400,v_diff, marker='+', label=locn) for ax in (ax_ssh, ax_sal, ax_u, ax_v): start, end = ax.set_xlim(0, 7) ax.set_xlabel('time') vmin=-0.005; vmax=0.005 ax_ssh.set_ylim(vmin,vmax); ax_ssh.set_ylabel('SSH (m)'); ax_u.set_ylim(-0.001,0.001); ax_u.set_ylabel('Surface u (m/s)'); ax_v.set_ylim(-0.001,0.001); ax_v.set_ylabel('Surface v (m/s)'); ax_sal.set_ylim(-.05,.05); ax_sal.set_ylabel('Surface Salinity'); ax_sal.legend(loc='best') fT1= NC.Dataset(runname1 + '/SalishSea_4h_20020915_20020921_grid_T.nc','r'); fU1= NC.Dataset(runname1 + '/SalishSea_4h_20020915_20020921_grid_U.nc','r'); fV1= NC.Dataset(runname1 + '/SalishSea_4h_20020915_20020921_grid_V.nc','r'); fT2= NC.Dataset(runname2 + '/SalishSea_4h_20020915_20020921_grid_T.nc','r'); fU2= NC.Dataset(runname2 + '/SalishSea_4h_20020915_20020921_grid_U.nc','r'); fV2= NC.Dataset(runname2 + '/SalishSea_4h_20020915_20020921_grid_V.nc','r'); u1 = fU1.variables['vozocrtx']; v1 = fV1.variables['vomecrty']; ssh1 = fT1.variables['sossheig']; u2 = fU2.variables['vozocrtx']; v2 = fV2.variables['vomecrty']; ssh2 = fT2.variables['sossheig']; times_full = fT2.variables['time_counter']; #Calculate max difference over the field for all times. fig =plt.figure(figsize=(12,12)) print ' Plotting difference in ssh over time: ' + r1 + ' - ' + r2 sshmaxdiff=np.zeros((42,1)); for ii in range(0,42): sshdiff =ssh1[ii,:,:]-ssh2[ii,:,:]; sshmaxdiff[ii] = abs(sshdiff[:]).max() fig =plt.figure(figsize=(12,6)) pylab.plot(times_full[:]/86400,sshmaxdiff) pylab.xlabel('time (d)') pylab.ylabel('difference in ssh (m)') pylab.title('Maximum difference in ssh') #Calculate max difference over the field for all times. fig =plt.figure(figsize=(12,12)) print ' Plotting difference in u over time: ' + r1 + ' - ' + r2 umaxdiff=np.zeros((42,1)); for ii in range(0,42): udiff =u1[ii,:,:,:]-u2[ii,:,:,:]; umaxdiff[ii] = abs(udiff[:]).max() fig =plt.figure(figsize=(12,6)) pylab.plot(times_full[:]/86400,umaxdiff) pylab.xlabel('time (d)') pylab.ylabel('difference in u (m/s)') pylab.title('Maximum difference in u') #Calculate max difference over the field for all times. vmaxdiff=np.zeros((42,1)); for ii in range(0,42): vdiff =v1[ii,:,:,:]-v2[ii,:,:,:]; vmaxdiff[ii] = abs(vdiff[:]).max() print ' Plotting difference in v over time: ' + r1 + ' - ' + r2 fig =plt.figure(figsize=(12,6)) pylab.plot(times_full[:]/86400,vmaxdiff) pylab.xlabel('time (d)') pylab.ylabel('difference in v (m/s)') pylab.title('Maximum difference in v')