% % airtrend.m Computes trend from MLOST v4 data done June 2015 % tried to update for latest data July 4, 2022 % compute trend of giseof.m % program to use sst.mnmean_ersst.nc Smith NOAA web site %data starts in 1880 montly means 2-degree resolution 88N-88S 0-358 %close all clear all % PLAN % make a data set of just North Pacific, or Atlantic data % a two dimensional array with basin points in one direction and time in the other % Compute EOFs and PCs of this % Regress PCs on the whole map of SST to show patterns Would help to have land mask I think % Suggest loading topo map, interpolating to SST grid, and using elevation aas land-sea mask ilag=0 % We're going to look at the SST regression before the anomaly develops icor= 0 % this is to make correlation maps instead of regression maps if = 1 jfilt = 0 detrend = 1 gmremove = 1 montrend = 0 ivmx=0 fcut = 72.; yrbeg= 1950; yrbeg= 1900; yrbeg= 1998; yrbeg= 1880; yrbeg=1901; yrbeg= 1979; yskip= (yrbeg-1880)*12 it1=1+yskip % sub select season you want % dmo1= 5; dm02 = 9 % dmo1= 1; dm02 = 12 dmo1= 10; dm02=3 dmo2=dm02; if dm02 < dmo1 dmo2=dm02+12; end im2=dmo2; im1=dmo1; % OK we are ready to make our array, using some prescribed lat long zones and our land mask % below for equatorial ENSO xlat1 = 10.; xlat2 = -10.; xlon1 = 120; xlon2 = 280.; % below for all Tropics xlat1 = 10.; xlat2 = -10.; xlon1 = 45; xlon2 = 280.; % Below for full pacific N of 30S xlat1 = 65.; xlat2 = -30.; xlon1 = 120; xlon2 = 265.; % Below for Global Pacific ST xlat1 = 87.; xlat2 = -87.; xlon1 = 1; xlon2 = 285.0; % Below for Global SST xlat1 = 70.; xlat2 = -70.; xlon1 = 1; xlon2 = 357.0; % Below for Atlantic SST xlat1 = 80.; xlat2 = -20.; xlon1 = 285; xlon2 = 357.0; % Below for PDO region N of 20N xlat1 = 20.; xlat2 = 75.; xlon1 = -179; xlon2 = 179.; iymx=160; yint=32 delda = yint * 12 %iymx=1 coast = load( 'coast.mat'); %filenme =['sst.mnmean_2.nc'] % Through January of 2015 %filenme =['gistemp1200_ERSST.nc'] % Through January of 2015 %filenme =['sst.mnmean.v4.nc']; % Through May of 2015 filenme =['air.mon.anomv3.5.4.nc']; % Through May of 2015 version='v4' if version=='v5' filenme =['air.mon.anom.v5.nc']; % Through May of 2022 elseif version=='v4' filenme =['air.mon.anom.nc']; % Through May of 2022 end fillmo=7; % months added to make integer years % Lats go down from 88N at 2.0 degrees Longs are 0-358 %xlat1=5.0; xlat2=-5.0; %ila1=int32(1+(88-xlat1)/2.); ila2=int32(1+(88-xlat2)/2.); % Specifies indices time = ncread(filenme,'time'); timxo=length(time) timx=length(time)-(yskip); %year = time/365.242199 + 1800; %time=(time-time(1)); %now in months lat = ncread(filenme,'lat'); lon = ncread(filenme,'lon'); ssto = ncread(filenme,'air'); latx = int32(length(lat)); lonx = int32(length(lon)); lat = double(lat); dlon=lon(2)-lon(1) lonxp = lonx+1 lonp=cat(1,lon,lon(lonx)+dlon); lonp= double(lonp); ssto2=ssto(:,:,it1:timxo); timx=size(ssto2,3); % make even number of years by adding missing values to end ssta=NaN(lonx,latx,timx+fillmo); %modified to stop in April 2015, so need extra 10 months for even years. ssta(:,:,1:timx)=ssto2; sst=ssta; timx=timx+fillmo % modified to include April 2015 %year = [1:timx]/12.+yrbeg; month = mod(1:timx,12); year = [1:timx]/12.+yrbeg; month(month<1)=12; clear ssto: clear ssto2: clear ssta; maxsst=max(max(max(sst))) timx=int32(size(sst,3)) % replace missing values with NaN for i = 1:lonx for j = 1:latx for k = 1:timx if sst(i,j,k) < -100 sst(i,j,k)=NaN; end end end end sst=double(sst); maxsst=max(max(max(sst))) % here we develop a mask for topography load('topo.mat','topo') topo(topo<0.)=0.; lont=[0:359]; latt=[-89:90]; lontx=length(lont); lattx=length(latt); % need to make two-D arrays of lats and longs %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% [lontb, lattb] = ndgrid(lont, latt); [lonb, latb] = ndgrid(lon, lat); lontb=double(lontb); lonb=double(lonb); lattb=double(lattb); latb=double(latb); topoz=interp2(lontb',lattb',topo,lonb,latb); %figure %contour(lon,lat,topoz',[0 0]) % title('Interpolated Topography') % OK we are ready to make our array, using some prescribed lat long zones and our land mask %xlat1 = 60.; xlat2 = 20.; %xlon1 = 130; xlon2 = 255.; dlon=lon(2)-lon(1); dlat=lat(2)-lat(1); %ilat1= 1 + (xlat1-88.0)/dlat; ilat2= 1 + (xlat2-88.0)/dlat; ilat1= 1 + (xlat1+89.0)/dlat; ilat2= 1 + (xlat2+89.0)/dlat; %GISTEMP %ilon1= 1 + xlon1/dlon; ilon2= 1 + xlon2/dlon; ilon1= 1 + (xlon1+179.)/dlon; ilon2= 1 + (xlon2+179.)/dlon; ilat1=int32(ilat1); ilat2=int32(ilat2); ilon1=int32(ilon1); ilon2=int32(ilon2); % remove the global mean at each time step lat=double(lat); cs=cos(lat*pi/180.); ls=ones(lonx,1,'double'); [CS,LS]=meshgrid(cs,ls); CSr=reshape(CS,lonx*latx,1); sstr=reshape(sst,lonx*latx,timx); for im=1:timx weit(:,im)=isfinite(sstr(:,im)); mnwt(im)=squeeze(nansum(weit(:,im).*CSr,1)); mnsst(im)=squeeze(nansum(sstr(:,im).*CSr,1)); mnsst(im)=mnsst(im)/mnwt(im); ssta(:,im)=sstr(:,im)-mnsst(im); end sstdif=sstr-ssta; maxdif=squeeze(nanmax(nanmax(nanmax((sstdif))))); maxdif for im=1:timx mnsst2(im)=squeeze(nansum(ssta(:,im).*CSr,1)); mnsst2(im)=mnsst2(im)/mnwt(im); %sstb(:,im)=ssta(:,im)-mnsst2(im); end if gmremove == 1 sstr=ssta; end figure % Fig. 1 line plot 1 plot(year,mnsst) title ('Global mean at each month') maxx=max(mnsst)*1.001; minn=min(mnsst)*0.999; axis([yrbeg 2020 minn maxx]); figure % Fig. 2 line plot plot(year,mnsst2) title ('Global mean at each month second time') maxx=max(mnsst2)*1.001; minn=min(mnsst2)*1.001; axis([yrbeg 2020 minn maxx]); % Before we do anything else we want to remove the climatology, errors and weight anomalies by sqrt(cos(lat) % compute and remove climatology whos sst yrmx=int32(timx/12) sstr = reshape(sst,lonx*latx,12,yrmx); sstc = nanmean(sstr,3); for j=1:12 for k=1:yrmx sstr(:,j,k)=sstr(:,j,k)-sstc(:,j); end end maxsstr=max(max(max(sstr))) % Here I should plot the standard deviation after the climatological annual cycle is removed. % * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * % First compute for all months and all years sstrg=reshape(sstr,lonx*latx,12*yrmx); sstst1=nanstd(sstrg,0,2); % Let's plot the variance of SST, Just for fun % sstv=nanstd(ssts,0,2); maxx=squeeze(max(max(abs(sstst1)))) % maxx=2.0 sstd=reshape(sstst1, lonx,latx); sstdp=cat(1,sstd,sstd(1,:)); % cint=0.15 cint = maxx/10. maxx = 4.0 cint=0.5 conts=[cint:cint:maxx]; maxe=maxx; pslat1=-90.; pslat2=90.0; FH =figure; % Try doing a polar stereographic projection Figure 3 ax = axesm('MapProjection','robinson','Grid','on','Frame','on','FLineWidth',1.0,'MapLatLimit',[pslat1 pslat2],'MapLonLimit',[180 540]); %plotm(coast.lat,coast.long,'Color','k','LineWidth',0.4); %contourm(lat,lonp,regp',contm,'linecolor','magenta','LineWidth',0.50); [cmap]=buildcmap('wr'); colormap(cmap); [C, H] = contourfm(lat,lonp,sstdp',conts,'linewidth',0.5); caxis([-0.01 maxe]); %,'LineWidth',0.5,'Fill','on','LineColor','none'); %h5 = clabelm(C,H); %,'Rotation',0); %,'LabelSpacing',10); hold on %plotm(coast.lat,coast.long,'Color',[0.5 0.5 0.5],'LineWidth',0.3); %geoshow(coast.lat,coast.long,'DisplayType','polygon','FaceColor','none') geoshow(coast.lat,coast.long,'DisplayType','polygon','FaceColor','none') set(gca,'FontSize',12); colorbar('SouthOutside') ax = axesm('MapProjection','robinson','Grid','on','Frame','on','FLineWidth',1.0,'MapLatLimit',[pslat1 pslat2],'MapLonLimit',[180 540]); title(['STD Surf Temp all Monthly ',num2str(yrbeg),', Season ',num2str(dmo1),'-',num2str(dm02),', ',num2str(year(1),4),'Filt=',num2str(jfilt),num2str(icor),' Lag=',num2str(ilag)],'FontSize',9); set(gca,'Color','none') hold off saveas(gca,['SSTSIG_Monthly4_',num2str(yrbeg),'_All',num2str(jfilt),'_',num2str(year(1),4),'_',num2str(dmo1),'-',num2str(dm02),'-',num2str(icor),num2str(ilag),'.epsc']); export_fig ['SSTSIG_Monthly4_',num2str(yrbeg),'_All',num2str(jfilt),'_',num2str(year(1),4),'_',num2str(dmo1),'-',num2str(dm02),'-',num2str(icor),num2str(ilag),'.epsc'] -transparent % Next just for the NH summer JJA sstrg=reshape(sstr(:,6:8,:),lonx*latx,3*yrmx); sstst1=nanstd(sstrg,0,2); % Let's plot the variance of SST, Just for fun % sstv=nanstd(ssts,0,2); maxx=squeeze(max(max(abs(sstst1)))) % maxx=2.0 sstd=reshape(sstst1, lonx,latx); sstdp=cat(1,sstd,sstd(1,:)); % cint=0.15 cint = maxx/10. maxx = 4.0 cint=0.5 conts=[cint:cint:maxx]; maxe=maxx*0.8; pslat1=-90.; pslat2=90.0; FH =figure; % Try doing a polar stereographic projection Fig. 4 ax = axesm('MapProjection','robinson','Grid','on','Frame','on','FLineWidth',1.0,'MapLatLimit',[pslat1 pslat2],'MapLonLimit',[180 540]); %plotm(coast.lat,coast.long,'Color','k','LineWidth',0.4); %contourm(lat,lonp,regp',contm,'linecolor','magenta','LineWidth',0.50); [cmap]=buildcmap('wr'); colormap(cmap); [C, H] = contourfm(lat,lonp,sstdp',conts,'linewidth',0.5); caxis([-0.01 maxe]); %,'LineWidth',0.5,'Fill','on','LineColor','none'); %h5 = clabelm(C,H); %,'Rotation',0); %,'LabelSpacing',10); hold on %plotm(coast.lat,coast.long,'Color',[0.5 0.5 0.5],'LineWidth',0.3); %geoshow(coast.lat,coast.long,'DisplayType','polygon','FaceColor','none') geoshow(coast.lat,coast.long,'DisplayType','polygon','FaceColor','none') set(gca,'FontSize',12); colorbar('SouthOutside') ax = axesm('MapProjection','robinson','Grid','on','Frame','on','FLineWidth',1.0,'MapLatLimit',[pslat1 pslat2],'MapLonLimit',[180 540]); title(['STD Surf Temp JJA ',num2str(yrbeg),' Season ',num2str(dmo1),'-',num2str(dm02),', ',num2str(year(1),4),'Filt=',num2str(jfilt),num2str(icor),' Lag=',num2str(ilag)],'FontSize',9); set(gca,'Color','none') hold off saveas(gca,['STSIG_1901_JJA',num2str(jfilt),'_',num2str(year(1),4),'_',num2str(dmo1),'-',num2str(dm02),'-',num2str(icor),num2str(ilag),'.epsc']); export_fig ['STSIG_1901_JJA',num2str(jfilt),'_',num2str(year(1),4),'_',num2str(dmo1),'-',num2str(dm02),'-',num2str(icor),num2str(ilag),'.epsc'] -transparent % Completely Separate section to plot trend of annual means % First compute annual means I think my basic data is in sstr(lonx*latx,12,yrmx) % modified to do seasonal subsets July 4, 2022 %sstam=squeeze(nanmean(sstr,2)); sstr=reshape(sstr,lonx*latx,12*yrmx); yrs=yrmx; if im2>12 yrs=yrs-1; end im=im1 dmo=im2-im1; for i=1:yrs sstam(:,i)=nanmean(sstr(:,im:im+dmo),2); im=im+12; end % end of years cycle yrmx=yrs; % OK now compute standard deviation of annual means stdam=nanstd(sstam,0,2); sstst1=stdam; sstd=reshape(sstst1, lonx,latx); sstdp=cat(1,sstd,sstd(1,:)); % cint=0.15 cint = maxx/10. annmin = min(min(sstd)) maxx = 2.0 cint=0.2 conts=[-cint:cint:maxx]; maxe=maxx; pslat1=-90.; pslat2=90.0; FH =figure; % Try doing a polar stereographic projection fig. 5 Standard Deviation of Annual Means ax = axesm('MapProjection','robinson','Grid','on','Frame','on','FLineWidth',1.0,'MapLatLimit',[pslat1 pslat2],'MapLonLimit',[180 540]); %plotm(coast.lat,coast.long,'Color','k','LineWidth',0.4); %contourm(lat,lonp,regp',contm,'linecolor','magenta','LineWidth',0.50); [cmap]=buildcmap('wr'); colormap(cmap); [C, H] = contourfm(lat,lonp,sstdp',conts,'linewidth',0.5); caxis([-0.01 maxe]); %,'LineWidth',0.5,'Fill','on','LineColor','none'); %h5 = clabelm(C,H); %,'Rotation',0); %,'LabelSpacing',10); hold on %plotm(coast.lat,coast.long,'Color',[0.5 0.5 0.5],'LineWidth',0.3); %geoshow(coast.lat,coast.long,'DisplayType','polygon','FaceColor','none') geoshow(coast.lat,coast.long,'DisplayType','polygon','FaceColor','none') set(gca,'FontSize',12); colorbar('SouthOutside') ax = axesm('MapProjection','robinson','Grid','on','Frame','on','FLineWidth',1.0,'MapLatLimit',[pslat1 pslat2],'MapLonLimit',[180 540]); title(['STD Surf Temp Annual, Season ',num2str(dmo1),'-',num2str(dm02),', ',num2str(yrbeg),'Filt=',num2str(jfilt),num2str(icor),' Lag=',num2str(ilag)],'FontSize',9); set(gca,'Color','none') hold off saveas(gca,['SSTSIG_ANN_V4_',num2str(jfilt),'_',num2str(yrbeg),'_',num2str(dmo1),'-',num2str(dm02),'-',num2str(icor),num2str(ilag),'.epsc']); export_fig ['SSTSIG_ANN_V4_',num2str(jfilt),'_',num2str(yrbeg),'_',num2str(dmo1),'-',num2str(dm02),'-',num2str(icor),num2str(ilag),'.epsc'] -transparent t=double([1:yrmx]); t=(t-mean(t))/std(t); sstamm=nanmean(sstam,2); for it = 1:yrmx sstam(:,it)=sstam(:,it)-sstamm; end sstam(isnan(sstam))=0.0; % need to set anomalies to zero before computing trend, idiot %maxssts=max(max(max(ssts))) corr = sstam*t'/double(yrmx); sstrnd=corr*(t(yrmx)-t(1)); % This is the linear change over time. %% OK now to plot this corr=reshape(sstrnd,lonx,latx); %% can I smooth this a little % 1-2-1 for ismooth=1:1 corrs=corr; for j=1:latx for i=2:lonx-1 corrs(i,j)=(corr(i-1,j)+2*corr(i,j)+corr(i+1,j))*0.25; end corrs(1,j)=(corr(lonx,j)+2*corr(1,j)+corr(2,j))*0.25; corrs(lonx,j)=(corr(lonx-1,j)+2*corr(lonx,j)+corr(1,j))*0.25; end % zonal smooth completed %next meridional smooth for i=1:lonx for j=2,latx-1 corrs(i,j)=(corr(i,j-1)+2*corr(1,j)+corr(i,j+1))*0.25; end end % meridional smooth completed corr=corrs; end % end of multiple smooths corrp=cat(1,corr,corr(1,:)); maxx = max(max(corrp)) cint = maxx/10; maxx=3.0 cint=0.3 maxe=maxx; conts = [-maxx:cint:maxx]; FH =figure; % Try doing a polar stereographic projection of trend Fig 6 pslat1=-90. ; pslat2=90. ; ax = axesm('MapProjection','robinson','Grid','on','Frame','on','FLineWidth',1.0,'MapLatLimit',[pslat1 pslat2],'MapLonLimit',[180 540]); %plotm(coast.lat,coast.long,'Color','k','LineWidth',0.4); %contourm(lat,lonp,regp',contm,'linecolor','magenta','LineWidth',0.50); [cmap]=buildcmap('bwr'); colormap(cmap); [C, H] = contourfm(lat,lonp,corrp',conts,'linewidth',0.5); caxis([-maxe*1.02 maxe]); %,'LineWidth',0.5,'Fill','on','LineColor','none'); %h5 = clabelm(C,H); %,'Rotation',0); %,'LabelSpacing',10); hold on %plotm(coast.lat,coast.long,'Color',[0.5 0.5 0.5],'LineWidth',0.3); %geoshow(coast.lat,coast.long,'DisplayType','polygon','FaceColor','none') geoshow(coast.lat,coast.long,'DisplayType','polygon','FaceColor','none') set(gca,'FontSize',12); colorbar('SouthOutside') ax = axesm('MapProjection','robinson','Grid','on','Frame','on','FLineWidth',1.0,'MapLatLimit',[pslat1 pslat2],'MapLonLimit',[180 540]); title(['Surf Temp Trend Annual ',num2str(yrbeg),' Annual Means ',num2str(dmo1),'-',num2str(dm02),', ',num2str(year(1),4),'Filt=',num2str(jfilt),num2str(icor),' Lag=',num2str(ilag),version],'FontSize',9); set(gca,'Color','none') hold off saveas(gca,['ST_TREND_ANN_V4_',num2str(yrbeg),'_',num2str(jfilt),'_',num2str(year(1),4),'_',num2str(dmo1),'-',num2str(dm02),'-',num2str(icor),num2str(ilag),version,'.epsc']); export_fig ['ST_TREND_ANN_V4_,'num2str(yrbeg),'_',num2str(jfilt),'_',num2str(year(1),4),'_',num2str(dmo1),'-',num2str(dm02),'-',num2str(icor),num2str(ilag),version.'.epsc'] -transparent % Also should remove trend % OK, now climatology has been removed, Next remove trend if montrend ==1 %% OK I am going to remove the trend from each month separately t=double([1:yrmx]); t=(t-mean(t))/std(t); ssts=reshape(sstr, lonx*latx*12,yrmx); sstm=nanmean(ssts,2); for it = 1:yrmx ssts(:,it)=ssts(:,it)-sstm; end ssts(isnan(ssts))=0.0; % need to set anomalies to zero before computing trend, idiot maxssts=max(max(max(ssts))) corr = ssts*t'/double(yrmx); sstrnd=corr*t; if detrend == 1 ssts=ssts-sstrnd; % HERE is where TREND is removed end maxssts=max(max(max(ssts))) sst=reshape(ssts,lonx*latx,12*yrmx); ssts=sst; else t=double([1:timx]); t=(t-mean(t))/std(t); ssts=reshape(sstr,lonx*latx,timx); sstm=nanmean(ssts,2); for it = 1:timx ssts(:,it)=ssts(:,it)-sstm; end ssts(isnan(ssts))=0.0; % need to set anomalies to zero before computing trend, idiot maxssts=max(max(max(ssts))) corr = ssts*t'/double(timx); corp=corr*(t(timx)-t(1)); corll = reshape(corp,lonx,latx); corrp=cat(1,corll,corll(1,:)); fskip =0 ; if fskip == 1 % plot trend figure contourf(lon,lat,corll') title(['map of trend ',num2str(yrbeg)]) colorbar('SouthOutside') end %fskip maxx = max(max(corrp)) cint = maxx/10 maxe=maxx*0.8; maxx=3.0 cint=0.3 maxe=maxx; conts = [-maxx:cint:maxx]; FH =figure; % Try doing a polar stereographic projection of trend Fig. 6 pslat1=-90. ; pslat2=90. ; ax = axesm('MapProjection','robinson','Grid','on','Frame','on','FLineWidth',1.0,'MapLatLimit',[pslat1 pslat2],'MapLonLimit',[180 540]); %plotm(coast.lat,coast.long,'Color','k','LineWidth',0.4); %contourm(lat,lonp,regp',contm,'linecolor','magenta','LineWidth',0.50); [cmap]=buildcmap('bwr'); colormap(cmap); [C, H] = contourfm(lat,lonp,corrp',conts,'linewidth',0.5); caxis([-maxe*1.02 maxe]); %,'LineWidth',0.5,'Fill','on','LineColor','none'); %h5 = clabelm(C,H); %,'Rotation',0); %,'LabelSpacing',10); hold on %plotm(coast.lat,coast.long,'Color',[0.5 0.5 0.5],'LineWidth',0.3); %geoshow(coast.lat,coast.long,'DisplayType','polygon','FaceColor','none') geoshow(coast.lat,coast.long,'DisplayType','polygon','FaceColor','none') set(gca,'FontSize',12); colorbar('SouthOutside') ax = axesm('MapProjection','robinson','Grid','on','Frame','on','FLineWidth',1.0,'MapLatLimit',[pslat1 pslat2],'MapLonLimit',[180 540]); title(['Surf Temp Trend Monthly ',num2str(yrbeg),' Season ',num2str(dmo1),'-',num2str(dm02),', ',num2str(year(1),4),'Filt=',num2str(jfilt),num2str(icor),' Lag=',num2str(ilag)],'FontSize',9); set(gca,'Color','none') hold off saveas(gca,['ST_TREND_MONTHLY_',num2str(yrbeg),'_',num2str(jfilt),'_',num2str(year(1),4),'_',num2str(dmo1),'-',num2str(dm02),'-',num2str(icor),num2str(ilag),'.epsc']); export_fig ['ST_TREND_MONTHLY_',num2str(yrbeg),'_',num2str(jfilt),'_',num2str(year(1),4),'_',num2str(dmo1),'-',num2str(dm02),'-',num2str(icor),num2str(ilag),'.epsc'] -transparent maxcorr=max(max(abs(corr))) sstrnd=corr*t; if detrend == 1 ssts=ssts-sstrnd; % HERE is where TREND is removed end maxssts=max(max(max(ssts))) sst=ssts; end % end of montrend if % OK if this has worked we have ssts(lonx*latx,timx) with seasons and trend removed and ready for something % Next thing would be filterin, then seasonal stratification std2plot=0 % do not plot variance with trend removed if std2plot == 1 % Let's plot the variance of SST, Just for fun sstv=nanstd(ssts,0,2); maxx=squeeze(max(max(abs(sstv)))) % maxx=2.0 sstd=reshape(sstv, lonx,latx); sstdp=cat(1,sstd,sstd(1,:)); % cint=0.15 cint = maxx/10. maxx = 4.0 cint=0.5 conts=[cint:cint:maxx]; maxe=maxx*0.8; pslat1=-90.; pslat2=90.0; FH =figure; % Try doing a polar stereographic projection ax = axesm('MapProjection','robinson','Grid','on','Frame','on','FLineWidth',1.0,'MapLatLimit',[pslat1 pslat2],'MapLonLimit',[180 540]); %plotm(coast.lat,coast.long,'Color','k','LineWidth',0.4); %contourm(lat,lonp,regp',contm,'linecolor','magenta','LineWidth',0.50); [cmap]=buildcmap('wr'); colormap(cmap); [C, H] = contourfm(lat,lonp,sstdp',conts,'linewidth',0.5); caxis([0 maxe]); %,'LineWidth',0.5,'Fill','on','LineColor','none'); %h5 = clabelm(C,H); %,'Rotation',0); %,'LabelSpacing',10); hold on %plotm(coast.lat,coast.long,'Color',[0.5 0.5 0.5],'LineWidth',0.3); %geoshow(coast.lat,coast.long,'DisplayType','polygon','FaceColor','none') geoshow(coast.lat,coast.long,'DisplayType','polygon','FaceColor','none') set(gca,'FontSize',12); colorbar('SouthOutside') ax = axesm('MapProjection','robinson','Grid','on','Frame','on','FLineWidth',1.0,'MapLatLimit',[pslat1 pslat2],'MapLonLimit',[180 540]); title(['STD Surf Temp -20A_1901_, Season ',num2str(dmo1),'-',num2str(dm02),', ',num2str(year(1),4),'Filt=',num2str(jfilt),num2str(icor),' Lag=',num2str(ilag)],'FontSize',9); set(gca,'Color','none') hold off saveas(gca,['SSTSIG_20P_1901_',num2str(jfilt),'_',num2str(year(1),4),'_',num2str(dmo1),'-',num2str(dm02),'-',num2str(icor),num2str(ilag),'.epsc']); export_fig ['SSTSIG_20P_1901_',num2str(jfilt),'_',num2str(year(1),4),'_',num2str(dmo1),'-',num2str(dm02),'-',num2str(icor),num2str(ilag),'.epsc'] -transparent end % End of if to skip plotting STD with trend removed % ******* ******** ************* ********* ********* if ivmx > 0 % if to skip all eigen analysis if jfilt == 1 % Next let's low-pass filter the data N=8 cut=2.0/fcut; % make cut at 2/Tscale , so that cut is at Tscale 30 day cut here [b,a]=butter(N,cut); [h,w]=freqz(b,a,128); r=abs(h); rlow=r.*r; f=w/pi/2.; %figure %plot(f,r,':r',f,rlow) %title(['Butterworth Response: N= ',num2str(N),' cut= ',num2str(cut)]) %xlabel('Cycles per time step') %ylabel('Response Function') %drawnow % OK, now we will try to filter the data sst for i = 1:lonx*latx sst90(i,:)=filtfilt(b,a,ssts(i,:)); end sst=sst90; clear sst90; clear ssts; else % this else is applied for jfilt not= 1 i.e.0 sst=ssts; clear ssts; end % end of filter if loop ******************************* % NOW seasonal subsetting delmo=dmo2-dmo1; if dmo2 > 12 yrmx=yrmx-1; end im1=1; im2=dmo2-dmo1+1; imo1=dmo1; imo2=dmo2; spcdim=int32(lonx*latx); timdim=int32(yrmx*(delmo+1)) z=zeros(spcdim,timdim,'double'); years=zeros(timdim,1,'double'); for iy = 1:yrmx z(:,im1:im2)=sst(:,imo1:imo2); years(im1:im2)=year(imo1:imo2); im1=im2+1; im2=im1+delmo; imo1=imo1+12; imo2=imo2+12; end % OK seasonal subset is in z(lonx*latx,??) % Next do sqrt(cos) weighting of anomalies sqc=sqrt(cos(lat*(pi/180.))); timxz=size(z,2) sst=reshape(z,lonx,latx,timxz); for il = 1:latx sst(:,il,:)=sst(:,il,:)*sqc(il); end % OK NOW I think we are ready to select regional subsets for EOF Analysis mask = zeros(lonx,latx); ind=1 for i=ilon1:ilon2 for j=ilat1:ilat2 ss(ind,:)=sst(i,j,:); mask(i,j)=1.0; ind=ind+1; end end ind if ilag == 0 figure contour(lon,lat,mask') hold on contour(lon,lat,topoz',[0 0],'Color','b') title('selection mask') maxss=max(max(max(ss))) ss(isnan(ss))=0.0; maxsss=max(max(max(ss))) end % now ss includes the region I want % Hey, before we do SVD, let's estimate the autocorrelation of the data NT = size(ss,2) aa = diag(ss*ss')/NT; % this is variance at each point ab = diag(ss(:,1:NT-1)*ss(:,2:NT)')/(NT-1); % one lag covariance ac=ab./aa; autt = nanmean(ac.*aa)/nanmean(aa); dof=NT*(1-autt*autt)/(1+autt*autt) [U S V] = svd(ss); sd=diag(S); % Get the sign right for 1950 to be consisten as El Nino if yrbeg == 1950 V(:,1)=-V(:,1); end sd=sd.*sd/double(ind); sd=sd/squeeze(sum(sd))*100; factor=sqrt(2./dof); sde=sd*factor; neofs = size(ss,1); indx=[1:neofs]; figure kk=12; errorbar(indx(1:kk),sd(1:kk),sde(1:kk),'LineWidth',1.5) % plot(sd(1:10)) title('Eigenvalue Spectrum','FontSize',16) ylabel('Percent of Variance','FontSize',16) xlabel('Mode Number','FontSize',16) set(gca,'FontSize',16) set(gca,'Color','none') saveas(gca,['EigenSpec_1515_00_',num2str(xlat1),'_',num2str(jfilt),'_',num2str(year(1),4),'_',num2str(dmo1),'-',num2str(dm02),'-',num2str(icor),num2str(ilag),'.epsc']); export_fig ['EigenSpec_1515_00_',num2str(xlat1),'_',num2str(jfilt),'_',num2str(year(1),4),'_',num2str(dmo1),'-',num2str(dm02),'-',num2str(icor),num2str(ilag),'.epsc'] -transparent % OK, time for regression % loop on eigenvalues %year=([1:timxz])/12. + 1854.; for iv = 1:ivmx % loop on eigenmodes *** ***** **** *** *** *** pc=V(:,iv); pc=(pc-mean(pc))/std(pc); if yrbeg == 1901 if xlat2 == -20 if iv == 1 pc=-pc; end end end if yrbeg == 1901 if xlat2 == 20 if iv == 2 pc=-pc; end end end if yrbeg == 1901 if detrend == 0 if iv < 3 pc=-pc; end end end if xlat2 == -30 if yrbeg == 1900 if iv <3 pc=-pc; end else if iv == 1 pc=-pc; end end end if xlat2 == -70. if iv < 1 pc=-pc; end if iv > 2 pc=-pc; end end if iv == 1 if yrbeg == 1900 pc=pc; end if yrbeg == 1950 pc=pc; end end if iv == 2 if yrbeg == 1900 pc=-pc; end end if yrbeg == 1950 pc=-pc; end if xlat1 == 30. if iv < 3 pc=-pc; end end if xlat2 == -70. if iv == 2 pc=-pc; end end if ilag == 0 % Let's add some code to filter the pc time series Wn=1.0/36.; nf = 7; [b,a] = butter(nf,Wn); [H,W] = freqz(b,a,100); %figure %plot(W,H); pcf = filtfilt(b,a,pc); figure ykmx=length(years); ykmx5=ykmx-11; % how do I knnow how many to cut off?? I am assuming data through Jan. 2015 plot(years(1:ykmx5),pc(1:ykmx5),'b',years(1:ykmx5),pcf(1:ykmx5),'r'); hold on plot(years(1:ykmx5),pcf(1:ykmx5),'r','LineWidth',1.5); title(['20P_1901_V(:,1)',' PC#',num2str(iv)]) set(gca,'FontSize',14); ylabel('PC Amplitude'); xlabel('Year') axis([yrbeg 2020 -3 4]); daspect([35,3.5,1]); daspect set(gca,'Color','none') saveas(gca,['PC_20P_1901_TimeSeries_',num2str(iv),'_',num2str(year(1)),'.epsc']); export_fig ['PC_20P_1901_Timeseries_',num2str(iv),'_',num2str(year(1)),'.epsc'] -transparent figure ykmxm=ykmx-72; plot(years(ykmxm:ykmx5),pc(ykmxm:ykmx5)); title(['20P_1901_V(:,1)',' PC#',num2str(iv)]) set(gca,'FontSize',20); ylabel('PC Amplitude'); xlabel('Year') axis([2010 2015 -3 4]); end % end of plot if ilag if ilag == 0 pc=pc'; fileout=['PDO_20P_1901_',num2str(yrbeg),'_NoTrend',num2str(iv),'.mat'] save (fileout,'pc','year'); % OK we can now compute the power spectra of the principle components, if we want m=256; n=m/2; m2m=int32(n/3); m2=m/2; h = spectrum.welch('hamming',m,50); P1 = pwelch(pc,m,n,m); f=(0:m2)/m; f=f*12; figure plot(f(1:m2m),P1(1:m2m),'Color','r','LineWidth',1.5) legend('Power Spectrum') title(['Power Spectrum','-',num2str(m),'-',num2str(year(1),4),' ',num2str(iv),]) xlabel('Frequency in Cycles per Year','FontSize',14) ylabel('Power','FontSize',14) set(gca,'FontSize',14) %xtick=[0:0.2:2]; %ytick=[0:1:3]; %set(gca,'Xtick',xtick); %set(gca,'Ytick',ytick); set(gca,'Color','none') saveas(gca,['PDO20P_1901_Power_',num2str(iv),'_',num2str(m),'_',num2str(year(1)),'.epsc']) export_fig ['PDO20P_1901_Power_',num2str(iv),'_',num2str(m),'_',num2str(year(1)),'.epsc'] -transparent end % of plot ilag if % regression sstr = reshape(sst,lonx*latx,timxz); if icor == 1 % make regression maps instead of regression maps % remove standard deviation from sstr sstrm=nanmean(sstr,2); sstrd=nanstd(sstr,1,2); for i = 1:timxz sstra(:,i)=(sstr(:,i)-sstrm)./sstrd; end sstr=sstra; end % end of correlation loop % OK, here we have an option to do lagged regression of SST on iTself NT=double(timxz); if ilag == 0 regr=sstr*pc'/NT; elseif ilag > 0 clear sstrs; clear pcs; NTM = NT -ilag; sstrs(:,1:NTM)=sstr(:,ilag+1:NT); pcs(1:NTM)=pc(1:NTM); regr=sstrs*pcs'/NTM; elseif ilag < 0 clear sstrs; clear pcs; jlag = abs(ilag) NTM = NT -jlag; sstrs(:,1:NTM)=sstr(:,1:NTM); pcs(1:NTM)=pc(jlag+1:NT); regr=sstrs*pcs'/NTM; end % end of LAG decision loop reg = reshape(regr,lonx,latx); if iv == 1 maxx=max(max(abs(reg)))*1.05 cint=maxx/10. %maxx=1.3 %cint=0.1 contm=[-maxx:cint:-cint,cint:cint:maxx]; % contm=[-maxx:cint:maxx]; maxe=maxx*0.7; conts=contm; if icor == 1 conts = [-1:0.1:1.0]; maxe=0.7; end end % end of icor==1 loop % maxx=max(max(abs(reg))) [cmap]=buildcmap('bwr'); colormap(cmap) %need to expand plotted arrays regp=cat(1,reg,reg(1,:)); pslat1=-90.; pslat2=90.0; FH =figure; % Try doing a polar stereographic projection ax = axesm('MapProjection','robinson','Grid','on','Frame','on','FLineWidth',1.0,'MapLatLimit',[pslat1 pslat2],'MapLonLimit',[180 540]); %plotm(coast.lat,coast.long,'Color','k','LineWidth',0.4); %contourm(lat,lonp,regp',contm,'linecolor','magenta','LineWidth',0.50); [cmap]=buildcmap('bwr'); colormap(cmap); [C, H] = contourfm(lat,lonp,regp',conts,'linewidth',0.5); [C, H] = contourm(lat,lonp,regp',[0. 0.],'linewidth',2.0); caxis([-maxe*1.04 maxe]); %,'LineWidth',0.5,'Fill','on','LineColor','none'); %h5 = clabelm(C,H); %,'Rotation',0); %,'LabelSpacing',10); hold on %plotm(coast.lat,coast.long,'Color',[0.5 0.5 0.5],'LineWidth',0.3); geoshow(coast.lat,coast.long,'DisplayType','polygon','FaceColor','none') set(gca,'FontSize',12); ax = axesm('MapProjection','robinson','Grid','on','Frame','on','FLineWidth',1.0,'MapLatLimit',[pslat1 pslat2],'MapLonLimit',[180 540]); title(['PDO Regression, Season ',num2str(dmo1),'-',num2str(dm02),' mode=',num2str(iv),', ',num2str(year(1),4),'Filt=',num2str(jfilt),num2str(icor),' Lag=',num2str(ilag)]); geoshow(coast.lat,coast.long,'DisplayType','polygon','FaceColor','none') set(gca,'Color','none') hold off saveas(gca,['PDO20P_1901_',num2str(iv),'_',num2str(jfilt),'_',num2str(year(1),4),'_',num2str(dmo1),'-',num2str(dm02),'-',num2str(icor),num2str(ilag),'.epsc']); export_fig ['PDO20P_1901_',num2str(iv),'_',num2str(jfilt),'_',num2str(year(1),4),'_',num2str(dmo1),'-',num2str(dm02),'-',num2str(icor),num2str(ilag),'.epsc'] -transparent end % end of eigenvalue loop end % end of if to skip all eigenanalysis *******************