%-------------------------------------------------------------------------- % %-------------------------------------------------------------------------- %Satellite Data %ftp://data.nodc.noaa.gov/pub/data.nodc/ghrsst/L4/GLOB/JPL/MUR/2004/064/FileSpecificMetadata/ filename=dir(); for f=3:96 com = sprintf('bunzip2 %s',filename(f).name); [stat,res] = system(com); ncid=netcdf.open(filename(f).name(1:end-4),'nowrite'); varid = netcdf.inqVarID(ncid,'lon'); lon = netcdf.getVar(ncid,varid); varid = netcdf.inqVarID(ncid,'lat'); lat = netcdf.getVar(ncid,varid); varid=netcdf.inqVarID(ncid,'analysed_sst'); sst = netcdf.getVar(ncid,varid); lon(lon<0)=lon(lon<0)+360; i=find(lon>240&lon<250); j=find(lat>23&lat<28); sub_sst=sst(i,j);clear sst sub_sst=double(sub_sst)*.001+298-273.15; eval(['SST',int2str(f),'.sst=sub_sst;']) eval(['SST',int2str(f),'.lat=lat(j);']) eval(['SST',int2str(f),'.lon=lon(i);']) eval(['SST',int2str(f),'.datetime=filename(f).name;']) [x,y]=meshgrid(lat(j),lon(i)); clf end %-------------------------------------------------------------------------- % Front Characterization %Front profiles twtz=[reshape(Front.Tw,[numel(Front.Tw) 1]) reshape(Front.Depth,[numel(Front.Depth) 1])]; tbins=5:.1:30; z=0:10:400; cycle_cutoff=.01; split_tw=nan(length(z),1);w_p=split_tw;c_p=w_p; for i=1:length(z)-1 ind=find(twtz(:,2)>=z(i)&twtz(:,2)10 %Fit histogram to distribution hdata=hist(twtz(ind,1),tbins);hdata=hdata/nansum(hdata); lp=[.001 max(twtz(ind,1))-0.5*(max(twtz(ind,1))-min(twtz(ind,1))) .01 .001 min(twtz(ind,1)) .01]; up=[10 max(twtz(ind,1)) 10 10 min(twtz(ind,1))+0.5*(max(twtz(ind,1))-min(twtz(ind,1))) 10]; foptions=fitoptions('Method','NonLinearLeastSquares','MaxIter',10000,'Robust','LAR','Lower',lp,'Upper',up); ff=fit(tbins',hdata','gauss2',foptions); cold_pdf=ff.a2*exp(-((3:.05:30-ff.b2)/ff.c2).^2); warm_pdf=ff.a1*exp(-((3:.05:30-ff.b1)/ff.c1).^2); if warm_pdf(tbins==ceil(ff.b2))<=cycle_cutoff && cold_pdf(tbins==floor(ff.b1))<=cycle_cutoff && ff.b2<=ff.b1 w_p(i)=ff.b1; c_p(i)=ff.b2; split_tw(i)=ff.b2+(ff.b1-ff.b2)/2; else w_p(i)=nan; c_p(i)=nan; split_tw(i)=nan; end end end hires_z=0:.1:400; ind=~isnan(w_p); w_p=smooth(w_p(ind)); w_p=interp1(z(ind),w_p,hires_z); ind=~isnan(c_p); c_p=smooth(c_p(ind)); c_p=interp1(z(ind),c_p,hires_z); ind=~isnan(split_tw); split_tw=smooth(split_tw(ind)); split_tw=interp1(z(ind),split_tw,hires_z); warm_mld=hires_z(find(w_pcold_mld+5 Front.box(i,f)=3; %if it was in between, then not in a region else Front.box(i,f)=0; end end end %determine if the temperature is warmer than split temperature if Front.Tw(i,f)>=split_tw(ind)+.25 %if it is, and above the mixed layer depth if Front.Depth(i,f)warm_mld+5 Front.box(i,f)=4; %if it was in between then not in a region else Front.box(i,f)=0; end end end %if the temperature was in the split temperature, then not in region if Front.Tw(i,f)>split_tw(ind)-.25 & Front.Tw(i,f)=244.6&sst.lon<=246.6); j=find(sst.lat>=25&sst.lat<=27); [x,y]=meshgrid(sst.lat(j),sst.lon(i)); cold_ind=find(roundn(sst.sst(i,j),-1)==roundn(max(c_p)+(max(w_p)-max(c_p))/2-.25,-1)); coldloc=double(deg2rad([x(cold_ind) y(cold_ind)])); warm_ind=find(roundn(sst.sst(i,j),-1)==roundn(max(c_p)+(max(w_p)-max(c_p))/2+.25,-1)); warmloc=double(deg2rad([x(warm_ind) y(warm_ind)])); for c=1:length(coldloc) dist=acos(sin(coldloc(c,1)).*sin(warmloc(:,1)) + cos(coldloc(c,1)).*cos(warmloc(:,1)).*cos(coldloc(c,2)-warmloc(:,2)))*6371; nnD(c)=min(dist); end dayDistStats(d,:)=[nanmin(nnD) prctile(nnD,5) prctile(nnD,25) nanmedian(nnD) prctile(nnD,75) prctile(nnD,95) nanmax(nnD)]; clear nnD end Front.CrossDistance=dayDistStats(:,2); % %-------------------------------------------------------------------------- % % Behavior at Front % % Distribution Front.DayDist.CS=nan(length(dates),4);Front.NightDist.CS=nan(length(dates),4); Front.DayDist.CD=nan(length(dates),4);Front.NightDist.CD=nan(length(dates),4); Front.DayDist.WS=nan(length(dates),4);Front.NightDist.WS=nan(length(dates),4); Front.DayDist.WD=nan(length(dates),4);Front.NightDist.WD=nan(length(dates),4); Front.DayDist.Transit=nan(length(dates),4);Front.NightDist.Transit=nan(length(dates),4); for f=1:4 for d=1:length(dates) dtotal=length(find(Front.dp(:,f)==1&round(Front.Time(:,f))==dates(d))); ntotal=length(find(Front.dp(:,f)==-1&round(Front.Time(:,f))==dates(d))); Front.DayDist.CS(d,f)=length(find(Front.box(:,f)==1&Front.dp(:,f)==1&round(Front.Time(:,f))==dates(d)))/dtotal; Front.NightDist.CS(d,f)=length(find(Front.box(:,f)==1&Front.dp(:,f)==-1&round(Front.Time(:,f))==dates(d)))/ntotal; Front.DayDist.CD(d,f)=length(find(Front.box(:,f)==3&Front.dp(:,f)==1&round(Front.Time(:,f))==dates(d)))/dtotal; Front.NightDist.CD(d,f)=length(find(Front.box(:,f)==3&Front.dp(:,f)==-1&round(Front.Time(:,f))==dates(d)))/ntotal; Front.DayDist.WS(d,f)=length(find(Front.box(:,f)==2&Front.dp(:,f)==1&round(Front.Time(:,f))==dates(d)))/dtotal; Front.NightDist.WS(d,f)=length(find(Front.box(:,f)==2&Front.dp(:,f)==-1&round(Front.Time(:,f))==dates(d)))/ntotal; Front.DayDist.WD(d,f)=length(find(Front.box(:,f)==4&Front.dp(:,f)==1&round(Front.Time(:,f))==dates(d)))/dtotal; Front.NightDist.WD(d,f)=length(find(Front.box(:,f)==4&Front.dp(:,f)==-1&round(Front.Time(:,f))==dates(d)))/ntotal; Front.DayDist.Transit(d,f)=length(find(Front.box(:,f)==0&Front.dp(:,f)==1&round(Front.Time(:,f))==dates(d)))/dtotal; Front.NightDist.Transit(d,f)=length(find(Front.box(:,f)==0&Front.dp(:,f)==-1&round(Front.Time(:,f))==dates(d)))/ntotal; end end %Movement Front.TransitStats=nan(5000,10);t_count=1; for f=1:4 %Calculate Switches res_loc=Front.box(1,f);t_time=1;reside_time=1;reside_ts=[Front.Tb(1,f) Front.Tb(1,f)-Front.Tw(1,f)];old_loc=nan; for i=2:length(Front.box(:,f)) switch Front.box(i,f)-res_loc case -res_loc %if transitioning t_time=t_time+1; case 0 %if staying reside_time=reside_time+1; t_time=1; reside_ts=[reside_ts;Front.Tb(i,f) Front.Tb(i,f)-Front.Tw(i,f)]; otherwise %if moving Front.TransitStats(t_count,:)=[Front.Time(i,f) f old_loc res_loc Front.box(i,f) Front.dp(i,f) median(reside_ts(:,1)) median(reside_ts(:,2)) t_time reside_time]; t_time=1; t_count=t_count+1; old_loc=res_loc; res_loc=Front.box(i,f); reside_time=1;reside_ts=[Front.Tb(i,f) Front.Tb(i,f)-Front.Tw(i,f)]; end end end %Paths Taken for f=1:4 for d=1:length(dates) Front.DayMove.TotalTransits(d,f)=length(find(round(Front.TransitStats(:,1))==dates(d)&Front.TransitStats(:,2)==f&Front.TransitStats(:,6)==1)); Front.NightMove.TotalTransits(d,f)=length(find(round(Front.TransitStats(:,1))==dates(d)&Front.TransitStats(:,2)==f&Front.TransitStats(:,6)==-1)); Front.DayMove.AB(d,f)=length(find(round(Front.TransitStats(:,1))==dates(d)&Front.TransitStats(:,2)==f&Front.TransitStats(:,4)==1&Front.TransitStats(:,5)==2&Front.TransitStats(:,6)==1))/Front.DayMove.TotalTransits(d,f); Front.DayMove.AC(d,f)=length(find(round(Front.TransitStats(:,1))==dates(d)&Front.TransitStats(:,2)==f&Front.TransitStats(:,4)==1&Front.TransitStats(:,5)==3&Front.TransitStats(:,6)==1))/Front.DayMove.TotalTransits(d,f); Front.DayMove.AD(d,f)=length(find(round(Front.TransitStats(:,1))==dates(d)&Front.TransitStats(:,2)==f&Front.TransitStats(:,4)==1&Front.TransitStats(:,5)==4&Front.TransitStats(:,6)==1))/Front.DayMove.TotalTransits(d,f); Front.NightMove.AB(d,f)=length(find(round(Front.TransitStats(:,1))==dates(d)&Front.TransitStats(:,2)==f&Front.TransitStats(:,4)==1&Front.TransitStats(:,5)==2&Front.TransitStats(:,6)==-1))/Front.NightMove.TotalTransits(d,f); Front.NightMove.AC(d,f)=length(find(round(Front.TransitStats(:,1))==dates(d)&Front.TransitStats(:,2)==f&Front.TransitStats(:,4)==1&Front.TransitStats(:,5)==3&Front.TransitStats(:,6)==-1))/Front.NightMove.TotalTransits(d,f); Front.NightMove.AD(d,f)=length(find(round(Front.TransitStats(:,1))==dates(d)&Front.TransitStats(:,2)==f&Front.TransitStats(:,4)==1&Front.TransitStats(:,5)==4&Front.TransitStats(:,6)==-1))/Front.NightMove.TotalTransits(d,f); Front.DayMove.BA(d,f)=length(find(round(Front.TransitStats(:,1))==dates(d)&Front.TransitStats(:,2)==f&Front.TransitStats(:,4)==2&Front.TransitStats(:,5)==1&Front.TransitStats(:,6)==1))/Front.DayMove.TotalTransits(d,f); Front.DayMove.BC(d,f)=length(find(round(Front.TransitStats(:,1))==dates(d)&Front.TransitStats(:,2)==f&Front.TransitStats(:,4)==2&Front.TransitStats(:,5)==3&Front.TransitStats(:,6)==1))/Front.DayMove.TotalTransits(d,f); Front.DayMove.BD(d,f)=length(find(round(Front.TransitStats(:,1))==dates(d)&Front.TransitStats(:,2)==f&Front.TransitStats(:,4)==2&Front.TransitStats(:,5)==4&Front.TransitStats(:,6)==1))/Front.DayMove.TotalTransits(d,f); Front.NightMove.BA(d,f)=length(find(round(Front.TransitStats(:,1))==dates(d)&Front.TransitStats(:,2)==f&Front.TransitStats(:,4)==2&Front.TransitStats(:,5)==1&Front.TransitStats(:,6)==-1))/Front.NightMove.TotalTransits(d,f); Front.NightMove.BC(d,f)=length(find(round(Front.TransitStats(:,1))==dates(d)&Front.TransitStats(:,2)==f&Front.TransitStats(:,4)==2&Front.TransitStats(:,5)==3&Front.TransitStats(:,6)==-1))/Front.NightMove.TotalTransits(d,f); Front.NightMove.BD(d,f)=length(find(round(Front.TransitStats(:,1))==dates(d)&Front.TransitStats(:,2)==f&Front.TransitStats(:,4)==2&Front.TransitStats(:,5)==4&Front.TransitStats(:,6)==-1))/Front.NightMove.TotalTransits(d,f); Front.DayMove.CA(d,f)=length(find(round(Front.TransitStats(:,1))==dates(d)&Front.TransitStats(:,2)==f&Front.TransitStats(:,4)==3&Front.TransitStats(:,5)==1&Front.TransitStats(:,6)==1))/Front.DayMove.TotalTransits(d,f); Front.DayMove.CB(d,f)=length(find(round(Front.TransitStats(:,1))==dates(d)&Front.TransitStats(:,2)==f&Front.TransitStats(:,4)==3&Front.TransitStats(:,5)==2&Front.TransitStats(:,6)==1))/Front.DayMove.TotalTransits(d,f); Front.DayMove.CD(d,f)=length(find(round(Front.TransitStats(:,1))==dates(d)&Front.TransitStats(:,2)==f&Front.TransitStats(:,4)==3&Front.TransitStats(:,5)==4&Front.TransitStats(:,6)==1))/Front.DayMove.TotalTransits(d,f); Front.NightMove.CA(d,f)=length(find(round(Front.TransitStats(:,1))==dates(d)&Front.TransitStats(:,2)==f&Front.TransitStats(:,4)==3&Front.TransitStats(:,5)==1&Front.TransitStats(:,6)==-1))/Front.NightMove.TotalTransits(d,f); Front.NightMove.CB(d,f)=length(find(round(Front.TransitStats(:,1))==dates(d)&Front.TransitStats(:,2)==f&Front.TransitStats(:,4)==3&Front.TransitStats(:,5)==2&Front.TransitStats(:,6)==-1))/Front.NightMove.TotalTransits(d,f); Front.NightMove.CD(d,f)=length(find(round(Front.TransitStats(:,1))==dates(d)&Front.TransitStats(:,2)==f&Front.TransitStats(:,4)==3&Front.TransitStats(:,5)==4&Front.TransitStats(:,6)==-1))/Front.NightMove.TotalTransits(d,f); Front.DayMove.DA(d,f)=length(find(round(Front.TransitStats(:,1))==dates(d)&Front.TransitStats(:,2)==f&Front.TransitStats(:,4)==4&Front.TransitStats(:,5)==1&Front.TransitStats(:,6)==1))/Front.DayMove.TotalTransits(d,f); Front.DayMove.DB(d,f)=length(find(round(Front.TransitStats(:,1))==dates(d)&Front.TransitStats(:,2)==f&Front.TransitStats(:,4)==4&Front.TransitStats(:,5)==2&Front.TransitStats(:,6)==1))/Front.DayMove.TotalTransits(d,f); Front.DayMove.DC(d,f)=length(find(round(Front.TransitStats(:,1))==dates(d)&Front.TransitStats(:,2)==f&Front.TransitStats(:,4)==4&Front.TransitStats(:,5)==3&Front.TransitStats(:,6)==1))/Front.DayMove.TotalTransits(d,f); Front.NightMove.DA(d,f)=length(find(round(Front.TransitStats(:,1))==dates(d)&Front.TransitStats(:,2)==f&Front.TransitStats(:,4)==4&Front.TransitStats(:,5)==1&Front.TransitStats(:,6)==-1))/Front.NightMove.TotalTransits(d,f); Front.NightMove.DB(d,f)=length(find(round(Front.TransitStats(:,1))==dates(d)&Front.TransitStats(:,2)==f&Front.TransitStats(:,4)==4&Front.TransitStats(:,5)==2&Front.TransitStats(:,6)==-1))/Front.NightMove.TotalTransits(d,f); Front.NightMove.DC(d,f)=length(find(round(Front.TransitStats(:,1))==dates(d)&Front.TransitStats(:,2)==f&Front.TransitStats(:,4)==4&Front.TransitStats(:,5)==3&Front.TransitStats(:,6)==-1))/Front.NightMove.TotalTransits(d,f); end end %-------------------------------------------------------------------------- %Figure 2 for f=1:4 dtotal=length(find(Front.dp(:,f)==1)); zbins=[0:2:225]; d_chist(f,:)=hist(Front.Depth(Front.dp(:,f)==1&Front.box(:,f)==1|Front.box(:,f)==3&Front.dp(:,f)==1,f),zbins)./dtotal; d_whist(f,:)=hist(Front.Depth(Front.dp(:,f)==1&Front.box(:,f)==2|Front.box(:,f)==4&Front.dp(:,f)==1,f),zbins)./dtotal; d_thist(f,:)=hist(Front.Depth(Front.dp(:,f)==1&Front.box(:,f)==0,f),zbins)./dtotal; ntotal=length(find(Front.dp(:,f)==-1)); n_chist(f,:)=hist(Front.Depth(Front.dp(:,f)==-1&Front.box(:,f)==1|Front.box(:,f)==3&Front.dp(:,f)==-1,f),zbins)./ntotal; n_whist(f,:)=hist(Front.Depth(Front.dp(:,f)==-1&Front.box(:,f)==2|Front.box(:,f)==4&Front.dp(:,f)==-1,f),zbins)./ntotal; n_thist(f,:)=hist(Front.Depth(Front.dp(:,f)==-1&Front.box(:,f)==0,f),zbins)./ntotal; end %2c Movement paths path_stats=[mean(mean(Front.DayMove.AB)) sqrt(0.25*(sum([mean(Front.DayMove.AB)-mean(mean(Front.DayMove.AB))].^2)))];%Results in excel file. %-------------------------------------------------------------------------- %Stats for text trips_CrossFront=[mean(mean(Front.DayMove.TotalTransits+Front.NightMove.TotalTransits)) sqrt(0.25*(sum([mean(Front.DayMove.TotalTransits+Front.NightMove.TotalTransits)-mean(mean(Front.DayMove.TotalTransits+Front.NightMove.TotalTransits))].^2)))] prcDay_trips_CrossFront=[mean(mean(Front.DayMove.TotalTransits./(Front.DayMove.TotalTransits+Front.NightMove.TotalTransits))) sqrt(0.25*(sum([mean(Front.DayMove.TotalTransits./(Front.DayMove.TotalTransits+Front.NightMove.TotalTransits))-mean(mean(Front.DayMove.TotalTransits./(Front.DayMove.TotalTransits+Front.NightMove.TotalTransits)))].^2)))] Front.swim=[nan(1,4);abs(diff(Front.Depth))]; for f=1:4 for d=1:length(dates) ind=find(round(Front.Time(:,f))==dates(d)&Front.dp(:,f)==1); day_distance(d,f)=nansum(abs(Front.swim(ind,f))); ind=find(round(Front.Time(:,f))==dates(d)&Front.dp(:,f)==-1); night_distance(d,f)=nansum(abs(Front.swim(ind,f))); end end day_distStats=[mean(mean(day_distance)) sqrt(0.25*(sum([mean(day_distance)-mean(mean(day_distance))].^2)))]; night_distStats=[mean(mean(night_distance)) sqrt(0.25*(sum([mean(night_distance)-mean(mean(night_distance))].^2)))]; ws_day_perctime=[mean(mean(Front.DayDist.WS)) sqrt(0.25*(sum([mean(Front.DayDist.WS)-mean(mean(Front.DayDist.WS))].^2)))]; size(Front.TransitStats(Front.TransitStats(:,6)==-1,:)) night_distStats=[mean(mean(night_distance)) sqrt(0.25*(sum([mean(night_distance)-mean(mean(night_distance))].^2)))]; ws_night=[mean(mean(Front.NightDist.WS+Front.NightDist.WD)) sqrt(0.25*(sum([mean(Front.NightDist.WS+Front.NightDist.WD)-mean(mean(Front.NightDist.WS+Front.NightDist.WD))].^2)))]; for f=1:4 tb_CS2CD(f)=mean(Front.TransitStats(Front.TransitStats(:,3)==1&Front.TransitStats(:,6)==1&Front.TransitStats(:,4)==3&Front.TransitStats(:,2)==f,7)); tb_WS2CD(f)=mean(Front.TransitStats(Front.TransitStats(:,3)==2&Front.TransitStats(:,6)==1&Front.TransitStats(:,4)==3&Front.TransitStats(:,2)==f,7)); tx_CS2CD(f)=mean(Front.TransitStats(Front.TransitStats(:,3)==1&Front.TransitStats(:,6)==1&Front.TransitStats(:,4)==3&Front.TransitStats(:,2)==f,8)); tx_WS2CD(f)=mean(Front.TransitStats(Front.TransitStats(:,3)==2&Front.TransitStats(:,6)==1&Front.TransitStats(:,4)==3&Front.TransitStats(:,2)==f,8)); end tb_stats_CS=[mean(tb_CS2CD) sqrt(0.25*(sum([tb_CS2CD-mean(tb_CS2CD)].^2)))]; %-------------------------------------------------------------------------- %Figure 3 tagnum=[1246;1967;1973;390167]; for f=1:length(tagnum) eval(['tag=tag',int2str(tagnum(f)),';']) ss_index=find(tag.insitu.daypart==2); sr_index=find(tag.insitu.daypart==0); if ss_index(1)>sr_index(1) sr_index=sr_index(2:end); end figure;count=1; for n=1:length(sr_index) i=ss_index(n):sr_index(n); if round(mean(tag.insitu.date(i)))>=dates(1)&round(mean(tag.insitu.date(i)))<=dates(end) nightlyDate(count)=mean(tag.insitu.date(i)); ts=smooth(tag.insitu.tb(i),120); dts=diff(ts); nightlyTb(count)=mode(roundn(ts(ts<=mean(ts)+.05),-1)); tstart=find(ts(1:end-1)>nightlyTb(count)+.05&dts>0,1,'first'); tend=tstart+find(ts(tstart:end-1)<=nightlyTb(count)+.05&dts(tstart:end)<0,1,'first'); if tend-tstart>120 nightlyAUC(count)=sum(ts(tstart:tend)-nightlyTb(count)); else tstart=tend+find(ts(tend:end-1)>nightlyTb(count)+.05&dts(tend:end)>0,1,'first'); tend=tstart+find(ts(tstart:end-1)<=nightlyTb(count)+.05&dts(tstart:end)<0,1,'first'); if tend-tstart>120 nightlyAUC(count)=sum(ts(tstart:tend)-nightlyTb(count)); nightlyDate(count)=mean(tag.insitu.date(i)); else tstart=nan; tend=nan; end end plot(ts,'-r.') hold on plot(nightlyTb(count)*ones(length(ts),1),'-k','LineWidth',2) plot(tstart,nightlyTb(count),'ko','MarkerFace','w') plot(tend,nightlyTb(count),'ko','MarkerFace','w') ylim([16 25]) title(int2str(datevec(nightlyDate(count)))) pause clf count=count+1; end end eval(['Forage',int2str(f),'.tb=nightlyTb;']) eval(['Forage',int2str(f),'.auc=nightlyAUC;']) eval(['Forage',int2str(f),'.date=nightlyDate;']) clear nightlyTb nightlyAUC nightlyDate end