import matplotlib.pyplot as plt import pandas as pd from salishsea_tools import stormtools import datetime from dateutil import tz import numpy as np %matplotlib inline def date_parser(s): unaware =datetime.datetime.strptime(s, '%Y-%m-%d %H:%M') aware = unaware.replace(tzinfo=tz.tzutc()) return aware DFO_file='/data/nsoontie/MEOPAR/analysis/Nancy/tides/Observed-DFO.csv' DFO=pd.read_csv(DFO_file,parse_dates=[1],header=0,names=['station','time','tides','e1','e2'],date_parser=date_parser) my_file1='/data/nsoontie/MEOPAR/analysis/Nancy/tides/Point Atkinson_t_tide_compare8_31-Dec-2013_02-Jan-2015.csv' ttide1,msl1=stormtools.load_tidal_predictions(my_file1) my_file1='/data/nsoontie/anne_tides/Point Atkinson_atide_compare8_31-Dec-2013_02-Dec-2015.csv' ttide2,msl2=stormtools.load_tidal_predictions(my_file1) fig,ax=plt.subplots(1,1,figsize=(10,5)) ax.plot(DFO.time,DFO.tides,'-o',label='DFO') ax.plot(ttide1.time,ttide1.pred_all +msl1,'-o',label='pred_all') ax.plot(ttide1.time,ttide1.pred_8+msl1,label='pred_8') ax.plot(ttide2.time,ttide2.pred_all +msl2,'-o',label='pred_all annes') ax.plot(ttide2.time,ttide2.pred_8+msl2,label='pred_8 annes') ax.set_xlim(datetime.datetime(2014,12,7),datetime.datetime(2014,12,12)) ax.legend(loc=0) fig,ax=plt.subplots(1,1,figsize=(10,5)) ax.plot(DFO.time,DFO.tides,'-o',label='DFO') ax.plot(ttide1.time,ttide1.pred_all +msl1,'-o',label='pred_all') ax.plot(ttide1.time,ttide1.pred_8+msl1,label='pred_8') ax.plot(ttide2.time,ttide2.pred_all +msl2,'-o',label='pred_all annes') ax.plot(ttide2.time,ttide2.pred_8+msl2,label='pred_8 annes') ax.set_xlim(datetime.datetime(2014,12,8),datetime.datetime(2014,12,9)) ax.set_ylim([0,1]) ax.legend(loc=0) ax.grid() fig,ax=plt.subplots(1,1,figsize=(10,5)) ax.plot(DFO.time,DFO.tides,'-o',label='DFO') ax.plot(ttide1.time,ttide1.pred_all +msl1,'-o',label='pred_all') ax.plot(ttide1.time,ttide1.pred_8+msl1,label='pred_8') ax.plot(ttide2.time,ttide2.pred_all +msl2,'-o',label='pred_all annes') ax.plot(ttide2.time,ttide2.pred_8+msl2,label='pred_8 annes') ax.set_xlim(datetime.datetime(2014,12,8),datetime.datetime(2014,12,9)) ax.set_ylim([4,5]) ax.legend(loc=0) ax.grid() fig,ax=plt.subplots(1,1,figsize=(10,5)) ax.plot(DFO.time,DFO.tides,'-o',label='DFO') ax.plot(ttide1.time,ttide1.pred_all +msl1,'-o',label='pred_all') ax.plot(ttide1.time,ttide1.pred_8+msl1,label='pred_8') ax.plot(ttide2.time,ttide2.pred_all +msl2,'-o',label='pred_all annes') ax.plot(ttide2.time,ttide2.pred_8+msl2,label='pred_8 annes') ax.set_xlim(datetime.datetime(2014,12,8),datetime.datetime(2014,12,9)) ax.set_ylim([3,4]) ax.legend(loc=0) ax.grid() t1=datetime.datetime(2014,12,8); t2=datetime.datetime(2014,12,12) t1 = t1.replace(tzinfo=tz.tzutc()) t2 = t2.replace(tzinfo=tz.tzutc()) subDFO=np.array(DFO.tides) timDFO=np.array(DFO.time) indices = np.where(np.logical_and(timDFO >= t1, timDFO <= t2)) subDFO=subDFO[indices]; timDFO=timDFO[indices]; meanDFO=np.nanmean(subDFO); maxDFO=np.nanmax(subDFO); minDFO=np.nanmin(subDFO) submine=np.array(ttide2.pred_all); timmine=np.array(ttide2.time); indices = np.where(np.logical_and(timmine >= t1, timmine <= t2)) submine=submine[indices] +msl2; timmine=timmine[indices]; meanmine=np.nanmean(submine); maxmine=np.nanmax(submine); minmine=np.nanmin(submine) print 'DFO: Mean {}, Max {}, Min {}'.format(meanDFO,maxDFO,minDFO) print 'Mine: Mean {}, Max {}, Min {}'.format(meanmine,maxmine,minmine) print 'Differences: Mean {}, Max {}, Min {}'.format(meanDFO-meanmine,maxDFO-maxmine,minDFO-minmine) from salishsea_tools.nowcast import figures import datetime from glob import glob import os import netCDF4 as nc import scipy.io as sio %matplotlib inline def results_dataset(period, grid, results_dir): """Return the results dataset for period (e.g. 1h or 1d) and grid (e.g. grid_T, grid_U) from results_dir. """ filename_pattern = 'SalishSea_{period}_*_{grid}.nc' filepaths = glob(os.path.join(results_dir, filename_pattern.format(period=period, grid=grid))) return nc.Dataset(filepaths[0]) run_date = datetime.datetime(2015,2,8) # Results dataset location results_home = '/data/dlatorne/MEOPAR/SalishSea/nowcast/' results_dir = os.path.join(results_home, run_date.strftime('%d%b%y').lower()) # model winds model_path = '/ocean/sallen/allen/research/MEOPAR/Operational/' coastline = sio.loadmat('/ocean/rich/more/mmapbase/bcgeo/PNW.mat') grid_T_hr = results_dataset('1h', 'grid_T', results_dir) bathy = nc.Dataset('/data/nsoontie/MEOPAR/NEMO-forcing/grid/bathy_meter_SalishSea2.nc') from salishsea_tools.nowcast import analyze from salishsea_tools import tidetools, stormtools import matplotlib.pyplot as plt from dateutil import tz t_orig=datetime.datetime(2014,12,1); t_final=datetime.datetime(2015,1,31) path='/data/dlatorne/MEOPAR/SalishSea/nowcast/' files_all=analyze.get_filenames(t_orig,t_final,'1h','grid_T','nowcast') def get_new_tides(name): if name == 'Point Atkinson': path = ('/data/nsoontie/anne_tides/') filename ='{}_atide_compare8_31-Dec-2013_02-Dec-2015.csv'.format(name) else: # Tide file covers 2014 and 2015. Harmonics from a 2013 time series. path = ( '/data/nsoontie/MEOPAR/tools/SalishSeaTools/salishsea_tools/' 'nowcast/tidal_predictions/') filename = '{}_t_tide_compare8_31-Dec-2013_02-Dec-2015.csv'.format(name) tfile = os.path.join(path, filename) ttide, msl = stormtools.load_tidal_predictions(tfile) return ttide lats,lons=figures.station_coords() b,X,Y= tidetools.get_bathy_data(bathy) name='Point Atkinson' j,i=tidetools.find_closest_model_point(lons[name],lats[name],X,Y,b,) ssh, time=analyze.combine_files(files_all,'sossheig','None',j,i) ttide_new=get_new_tides(name) ttide_old = figures.get_tides(name) msl=3.09 sdt = t_orig.replace(tzinfo=tz.tzutc()) edt = t_final +datetime.timedelta(days=1) edt=edt.replace(tzinfo=tz.tzutc()) ssh_corr_new = stormtools.correct_model(ssh,ttide_new,sdt,edt) ssh_corr_old = stormtools.correct_model(ssh,ttide_old,sdt,edt) obs=stormtools.load_observations((t_orig-datetime.timedelta(days=1)).strftime('%d-%b-%Y'), (t_final+datetime.timedelta(days=1)).strftime('%d-%b-%Y'),'PointAtkinson') #obs=figures.load_PA_observations() fig,ax=plt.subplots(1,1, figsize=(15,5)) ax.plot(obs.time,obs.slev,'g',lw=1.5,label='observations') ax.plot(time,ssh+msl,'k--',lw=1.5,label='model') ax.plot(time,ssh_corr_new+msl,'b',lw=1.5,label='corrected model - anns tides') ax.plot(time,ssh_corr_old+msl,'r',lw=1.5,label='corrected model - my tides') plt.legend(loc=0) td=datetime.timedelta(days=17) win=datetime.timedelta(days=5) ax.set_xlim([t_orig+td,t_orig+td+win]) ax.grid() def interp_to_model_time(time_model,varp,tp): """ interpolate to model time """ # Strategy: convert times to seconds past a reference value. # Use this as the independent variable in interpolation. #Set epoc (reference) time. epoc=time_model[0]; # Determine tp times wrt epc tp_wrt_epoc=[] for t in tp: tp_wrt_epoc.append((t-epoc).total_seconds()) # Interpolate observations to model times varp_interp=[] for t in time_model: mod_wrt_epoc=(t-epoc).total_seconds(); varp_interp.append(np.interp(mod_wrt_epoc,tp_wrt_epoc,varp)) return varp_interp obs_interp=interp_to_model_time(time,obs.slev,obs.time) fig,ax=plt.subplots(1,1, figsize=(15,5)) ax.plot(time,obs_interp,'g',lw=1.5,label='observations') ax.plot(time,ssh+msl,'k--',lw=1.5,label='model') ax.plot(time,ssh_corr_new+msl,'b',lw=1.5,label='corrected model - anns tides') ax.plot(time,ssh_corr_old+msl,'r',lw=1.5,label='corrected model - my tides') plt.legend(loc=0) td=datetime.timedelta(days=17) win=datetime.timedelta(days=5) #ax.set_xlim([t_orig+td,t_orig+td+win]) ax.grid() errors={} errors['new']=ssh_corr_new+msl-obs_interp errors['old'] = ssh_corr_old +msl - obs_interp errors['none'] = ssh +msl -obs_interp for n in ['new','old','none']: print n, 'mean error:', np.mean(errors[n]), 'max error:', np.max(errors[n]), 'min error:', np.min(errors[n]) 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'}) filename = '/data/nsoontie/MEOPAR/analysis/Idalia/other_constituents.csv' harm_other = pd.read_csv(filename,sep=',',header=0) harm_other = harm_other.rename(columns={'Site': 'site', 'Lat': 'lat', 'Lon': 'lon', 'O1 amp': 'O1_amp', 'O1 phase (deg UT)': 'O1_pha', 'P1 amp': 'P1_amp', 'P1 phase (deg UT)': 'P1_pha', 'Q1 amp': 'Q1_amp', 'Q1 phase (deg UT)': 'Q1_pha', 'S2 amp': 'S2_amp', 'S2 phase (deg UT)': 'S2_pha', 'N2 amp': 'N2_amp', 'N2 phase (deg UT)': 'N2_pha', 'K2 amp': 'K2_amp', 'K2 phase (deg UT)': 'K2_pha'}) s='Point Atkinson' harm_other[harm_other.site==s] ampsF={}; phasF={} for t in ['K1','M2']: astr=t + '_amp' pstr=t + '_pha' ampsF[t]= harm_obs[harm_obs.site==s][astr].values phasF[t]= harm_obs[harm_obs.site==s][pstr].values tides=['O1','P1','Q1','S2','N2','K2'] for t in tides: astr=t + '_amp' pstr=t + '_pha' ampsF[t]= harm_other[harm_other.site==s][astr].values phasF[t]= harm_other[harm_other.site==s][pstr].values annes='/data/nsoontie/anne_tides/07795const.wlev' atides=pd.read_table(annes,skiprows=24,names=['per','amp','pha','a','b','c','d','e','f','g'],delim_whitespace=True) ampsA={};phasA={} tides=['K1','O1','P1','Q1','M2','S2','N2','K2'] for t in tides: pha=atides.ix[t].pha +8.0*360/atides.ix[t].per #convert to UTC ampsA[t]=atides.ix[t].amp*100 phasA[t]= pha -360*(pha>360) fig,axs=plt.subplots(2,1) for t in tides: axs[0].plot(ampsF[t],ampsA[t],'o',label=t) axs[1].plot(phasF[t],phasA[t],'o',label=t) axs[0].plot([0,100],[0,100],'r') axs[1].plot([0,360],[0,360],'r') axs[0].set_xlabel('Foreman Amp') axs[0].set_ylabel('Annes Amp') axs[1].set_xlabel('Foreman Phase') axs[1].set_ylabel('Annes Phase') axs[0].legend(loc=0) print 'Constitents, ampA, ampF, phaA, phaF' for t in tides: print t, ampsA[t], ampsF[t], phasA[t], phasF[t]