% sstsam13.m Program to composite sst on SAM based on monthly data % Modified from sstsam9.m to use data from 1959 ERA5 %major rehack to remove ENSO from subsets consistently % compoosites ERA5 on ERA5-derived index % this one uses ERA52m Temperature % April 2020 got latest data from psl.noaa.gov/datagriddeddata.noaa.ersse.v5.html See references there % based on ssteof.m Make eofs from SST in various basins % program to use sst.mnmean_ersst.nc Smith NOAA web site %data starts in 1854 montly means 2-degree resolution 88N-88S 0-358 % Dennis L. Hartmann January 2016 **1944 months in todays sst.mnmean.nc V3b %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; detrend1 = 1; gmremove = 1; montrend = 1; ivmx=3; fcut = 72.; ensoremove=1; % regress nino3.4 out of sst field if ensoremove==1 ensoremove2=1; var1='v10' var1='u10' var1='sst' var0='nino'; var0='ssts'; im1=1; im2=12; idof=24; im1=10; im2=15; idof=24; lag=1 yrstart=1979; yrstart=1959; % for extended ERA5 data set yrbeg= 2000; yrend= 1990; yrend= 1989; yrbeg= 1990; yrbeg= 1979; yrbeg= 1959; yrend= 2021; % 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 tropical Pacific xlat1 = 30.; xlat2 = -30.; xlon1 = 110; xlon2 = 280.; % Below for Global Pacific ST xlat1 = 87.; xlat2 = -87.; xlon1 = 1; xlon2 = 285.0; % below for all Tropics xlat1 = 20.; xlat2 = -20.; xlon1 = 1.; xlon2 = 357.; % Below for full pacific N of 65S xlat1 = 65.; xlat2 = -65.; xlon1 = 120; xlon2 = 255.; % Below for Rachel's Atlantic xlat1 = 65.; xlat2 = -65.; xlon1 = 230; xlon2 = 320.; % Below for Atlantic SST xlat1 = 80.; xlat2 = -65.; xlon1 = 285; xlon2 = 357.0; expname='Atlantic'; % Below for PDO region N of 20N xlat1 = 65.; xlat2 = 20.; xlon1 = 120; xlon2 = 260.; expname='PDORegion'; % Below for Global SST xlat1 = 70.; xlat2 = -70.; xlon1 = 0.; xlon2 = 360.0; expname='Global'; coast = load( 'coast.mat'); %filenme =['/home/disk/eos10/dennis/ERAUW.d/sst.mnmean.v4.nc'] % Through December of 2017 Version 4 %filenme =['/home/disk/eos10/dennis/ERAUW.d/sst.mnmeanV5.nc'] % Through January of 2021 Version 5 %filenme =['/home/disk/eos10/dennis/ERAUW.d/ERA5_2mTemp.nc'] % Through March of 2022 ERA5 monthly if var1=='u10' %filenme =['/home/disk/eos10/dennis/ERAUW.d/ERA5_10mU.nc'] % Through March of 2022 ERA5 monthly filenme =['/home/disk/eos10/dennis/ERAUW.d/ERA5_1959_u10.nc'] % Through March of 2022 ERA5 monthly elseif var1=='sst' %filenme =['/home/disk/eos10/dennis/ERAUW.d/ERA5_SST.nc'] % Through March of 2022 ERA5 monthly filenme =['/home/disk/eos10/dennis/ERAUW.d/ERA5_1959_sst.nc'] % Through March of 2022 ERA5 monthly elseif var1=='v10' filenme =['/home/disk/eos10/dennis/ERAUW.d/ERA5_1959_v10.nc'] % Through March of 2022 ERA5 monthly end % Lats go down from North Pole 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 yskip=0; 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,'latitude'); lon = ncread(filenme,'longitude'); % ssto = ncread(filenme,'t2m'); % ssto = ncread(filenme,'sst'); ssto = ncread(filenme,var1); expver=1; nsmooth=60; latx = int32(length(lat)); lonx = int32(length(lon)); % subsample quarter-degree to one degree ids=4; % pick every 4th data point %timx=size(ssto,3) %timx=timx-3 ; % data ends with march so subtract 3 to get on integer years timx=timx ; % 1959 data was cut off at integer 2021 data ends with march so subtract 3 to get on integer years ibeg=1; iend=timx; % here we are using all the data in integral years 1959-2021 in this case % OK, let's allow for dividing the sample sst=squeeze(ssto(1:ids:lonx,1:ids:latx,ibeg:iend)); % reduced volume of data note no expver SZ=size(sst); timx=SZ(3) lon=lon(1:ids:lonx); lat=lat(1:ids:latx); lonx=length(lon); latx=length(lat); lat = double(lat); lonxp = lonx+1; dlon=lon(2)-lon(1); lonp=cat(1,lon,lon(lonx)+dlon); lonp= double(lonp); %ssto2=ssto(:,:,it1:timxo); % make even number of years by adding missing values to end %ssta=NaN(lonx,latx,timx); %The plus zero here has to be changed depending on what month the series ends on %ssta(:,:,1:timx)=ssto2;sst=ssta; %sst(:,:,1:timx)=ssto(:,:,1:timx); 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 sst(sst< -100)= NaN; 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; ilon1= 1 + xlon1/dlon; ilon2= 1 + xlon2/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 sst=reshape(ssta,lonx,latx,timx); end iplothis=0; if iplothis==1 figure plot(year,mnsst) title ('Global mean at each month') maxx=max(mnsst)*1.001; minn=min(mnsst)*0.999; axis([yrbeg yrend minn maxx]); figure plot(year,mnsst2) title ('Global mean at each month second time') maxx=max(mnsst2)*1.001; minn=min(mnsst2)*1.001; axis([yrbeg yrend minn maxx]); end % iplothis % 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); % removes climatology end end maxsstr=max(max(max(sstr))) % 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 instead of NaN before computing trend, idiot maxssts=max(max(max(ssts))) corr = ssts*t'/double(yrmx); sstrend=corr*t; if detrend1 == 1 ssts=ssts-sstrend; % 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); plotrend=0; if plotrend==1 % plot trend figure corll = reshape(corr,lonx,latx); contourf(lon,lat,corll') title(['map of trend ',num2str(yrbeg)]) colorbar drawnow; end % of plotrend if maxcorr=max(max(abs(corr))) sstrend=corr*t; if detrend1 == 1 ssts=ssts-sstrend; % 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 % ****************************************************************************************************************** % read in raw SAM and ENSO data so that we can regress out ENSO from subsets of data filenso='/home/disk/eos10/dennis/nino/Nino_3-4.mat'; %filesam='/home/disk/eos10/dennis/ERAUW.d/SAM5-60u10m0-359.mat'; %filesam='/home/disk/eos10/dennis/ERAUW.d/SAM5-u10m.mat'; filesam='/home/disk/eos10/dennis/ERAUW.d/SAM5-59u10m0-359.mat'; load(filenso,'nino34','year'); load(filesam,'xc'); % ENSO indices start in 1950, so need to skip 29*12 timx=length(xc) % New Section to regress nino34 out of SST field prior to compositing on SAM % Next choose chunk of data we want to use ibeg=(yrbeg-yrstart)*12+1; iend=timx-(2021-yrend)*12; % Here we set indices for Nino3.4, which starts in 1950 ibeg2=(yrstart-1950)*12+1; iend2=length(nino34)-3; % nino34 starts in 1950 and goes through march of 2022 nino34=nino34(ibeg2:iend2); % lines nino34 up with entire xc sequence nino34=nino34(ibeg:iend); % picks subset of nino34 same as subset of xc xc=xc(ibeg:iend); % now need to also subset the input sst field ssts=ssts(:,ibeg:iend); year=year(ibeg2:iend2); year=year(ibeg:iend); mnx=iend-ibeg+1; yrmx=mnx/12; % Need to recompute yrmx here when we use less that whole ERA5 time series iplothis=0; if iplothis==1 figure plot(year,xc) title('wind index vs time') end % iplothis if % need to remove seasonal cycle from index for i=1:12 ann(i)=nanmean(xc(i:12:mnx)); ann2(i)=nanmean(nino34(i:12:mnx)); end xc=reshape(xc,12,yrmx); nino34=reshape(nino34,12,yrmx); for i=1:yrmx xc(:,i)=xc(:,i)-ann'; nino34(:,i)=nino34(:,i)-ann2'; end xc=reshape(xc,mnx,1); nino34=reshape(nino34,mnx,1); % need to subset what seasons we want here yrs=mnx/12 %xca=reshape(xc,12,yrs); %ssta=reshape(ssts,lonx*latx,12,yrs); %ninoa=reshape(nino34,12,yrmx); ima=im1+lag; imb=im2+lag; im=1; dm=im2-im1; if im2>12 yrs=yrs-1; end for iy=1:yrs xcs(im:im+dm)=xc((iy-1)*12+im1:(iy-1)*12+im2); years(im:im+dm)=year((iy-1)*12+im1:(iy-1)*12+im2); nino(im:im+dm)=nino34((iy-1)*12+im1:(iy-1)*12+im2); ninl(im:im+dm)=nino34((iy-1)*12+ima:(iy-1)*12+imb); sstc(:,im:im+dm)=ssts(:,(iy-1)*12+ima:(iy-1)*12+imb); im=im+dm; end % year loop % ok now we have subsetted to just the season we want xc=xcs; nino34=nino; nino3l=ninl; ssts=sstc; year=years; mnx=length(xc) if ensoremove2==1 % regress Nino3.4 out of SAM index code to regress out nino3.4 from SAM index %lets normalize both time series nino34=(nino34-mean(nino34))/std(nino34,1); xc=(xc-mean(xc))/std(xc,1); %compute correlation cor=xc*nino34'/length(xc) %whos cor xcn=xc-cor*nino34; iplothis=0; if iplothis==1 figure plot(year,xcn,year,xc,year,nino34) legend('xcn','xc','nino34','Location','best'); xlabel('year') ylabel('Normalized Variable') figure end % ipolothis xc=xcn; end % ensoremove2 % OK this is where we are going to start with our SAM compositing We will want to cycle over all months starting in 1979, % ending in 2018 40 years % read in SAM index *** These come from samenso.m where the effect of enso has been removed. level='1000'; var='u10m'; % load (fileout,'lat1','lat2','lat3','lat4','xc'); % fileout=['SAM5-',var,'150-290.mat']; %fileout=['SAM5-Nino-6--0.18.mat']; % fileout=['SAM5-',var,'.mat']; %fileout=['SAM5-Nino-6--0.28.mat']; %fileout=['SAM5-60u10m0-359.mat']; %fileout=['SAM5-Nino-3-1-12--0.16.mat']; %fileout=['SAM5-Nino-0-4-9-0.069.mat']; %fileout=['SAM5-Nino-0-1-12--0.1.mat']; %fileout=['SAM5-Nino-3-10-15--0.27.mat']; %im1=4; im2=9; %fileout=['SAM5-Nino-0-4-9--0.069.mat']; %fileout=['SAM5-Nino-0-10-13--0.21.mat']; %im1=1; im2=12; idof=81; %fileout=['SAM5-Nino-0-1-12--0.1.mat']; %fileout=['SAM5-Nino-0-10-15--0.2.mat']; % New Section to regress nino34 out of SST field prior to compositing on SAM %Load SST index % filenino='/home/disk/eos10/dennis/nino/Nino_3-4.mat'; %load (filenino, 'nino34','year'); %ibeg2=29*12+1; iend2=length(nino34)-3; % nino34 starts in 1950 and goes through march of 2022 %nino34=nino34(ibeg2:iend2); %nino34=nino34(ibeg:iend); %year=year(ibeg2:iend2); %year=year(ibeg:iend); if ensoremove==1 xx=(nino3l-mean(nino3l))/std(nino3l); xx=xx'; cof=xx'*ssts'/length(xx); % should be regressions of nino34 onto SST field en34=xx*cof; en34=en34'; ssts=ssts-en34; % here I remove the regressed enso time series from the original data iplothis=0; if iplothis==1 cof=reshape(cof,lonx,latx); % now plot the Nino 3.4 regression coef maxx=squeeze(max(max(abs(cof)))) cint=maxx/5; % maxx=2.0 sstdp=cat(1,cof,cof(1,:)); lon=double(lon); dlon=lon(2)-lon(1); lonp=cat(1,lon,lon(lonx)+dlon); lat=double(lat); cint=maxx/7.; conts=[-maxx:cint:maxx]; maxe=maxx; sstdp=sm121(sstdp,lonp,lat,nsmooth); pslat1=-80.; pslat2=80.0; iplotthis=0 if iplotthis==1 FH = figure; % Try doing a polar stereographic projection ax = axesm('MapProjection','robinson','Grid','on','Frame','on','FLineWidth',1.0,'MapLatLimit',[pslat1 pslat2],'MapLonLimit',[30 390]); [cmap]=buildcmap('bwr'); colormap(cmap); [C, H] = contourfm(lat,lonp,sstdp',conts,'linewidth',0.5); caxis([-maxe maxe]); hold on geoshow(coast.lat,coast.long,'DisplayType','polygon','FaceColor','black') set(gca,'FontSize',12); colorbar('Location','SouthOutside'); ax = axesm('MapProjection','robinson','Grid','on','Frame','on','FLineWidth',1.0,'MapLatLimit',[pslat1 pslat2],'MapLonLimit',[30 390]); title(['SST Nino3.4',num2str(yrbeg),'-',' Season ',num2str(1),'-',num2str(12),', ',num2str(yrend,4),' L',num2str(lag)]); geoshow(coast.lat,coast.long,'DisplayType','polygon','FaceColor','black') set(gca,'Color','none') hold off saveas(gca,['SSTE34_G',var1,num2str(yrbeg),'_',num2str(jfilt),'_',num2str(yrend,4),'_',num2str(im1),'-',num2str(im2),'-',num2str(icor),num2str(gmremove),num2str(lag),'.epsc']); export_fig ['SSTE34_G',var1,num2str(yrbeg),'_',num2str(jfilt),'_',num2str(yrend,4),'_',num2str(im1),'-',num2str(im2),'-',num2str(icor),num2str(gmremove),num2str(lag),'.epsc'] -transparent; end % iplotthis end % iplothis NINO34 Pattern end % ensoremove if idothis=0; if idothis==1 load (fileout,'xc','xcr') if ensoremove==0 xc=xcr(ibeg:iend); end if var0=='nino' xc=-nino34; end mnx=length(xc); yrmx=mnx/12; % Need to normalized xc xc=xc-mean(xc); end % if idothis xc=xc/std(xc,0); xc=detrend(squeeze(xc)); [ac, lags]=autocorr(xc); acorel=ac(2); imx=mnx; sdcrit= 1.0; %make storage arrays n=zeros(2,1); xm=zeros(lonx*latx,2); % Compute autocorrelation array for composited field acxc=autocorr(xc); ac=zeros(lonx*latx,1); yrmx=timx/12; for i=1:lonx*latx xx=autocorr(squeeze(ssts(i,:))','NumLags',1); ac(i)=xx(2); end hits=zeros(yrmx*12,1); imp=-1 imn=-1 % changed this since only the months I want are now present for im=1:mnx if xc(im)>sdcrit hits(im)=1; %if im>imp+1 % test to avoid sequential months xm(:,1)=xm(:,1)+ssts(:,im); n(1)=n(1)+1; nn(n(1),1)=(iy-1)*12+im; imp=im; %end elseif xc(im)<-sdcrit hits(im)=-1; %if im>imn+1 % test to avoid sequential months xm(:,2)=xm(:,2)+ssts(:,im); n(2)=n(2)+1; nn(n(2),2)=(iy-1)*12+im; imn=im; %end end % if on sdcrit end %imx idof=min(n)*(1.-acorel)/(1+acorel); %figure %plot(year,hits) %title('Hits vs year') %figure % OK now take difference and plot xmd=xm(:,1)/n(1)-xm(:,2)/n(2); xmp=reshape(xmd,lonx,latx); % Let's plot the variance of SST, Just for fun sstd=nanstd(ssts,0,2); % standard deviation of SST % standard error % sstv=sstv/sqrt(double(idof)); % I am assuming a composit size of 40 for sdcrit=1 sste=sstd*(sqrt(1/n(1))+sqrt(1./n(2))); % use t-statistic assuming each sample is independent % sstd=sstv*(sqrt(1/n(1))); %+sqrt(1./n(2))); % use t-statistic assuming each sample is independent maxx=squeeze(max(max(abs(sstd)))) disp('maxx std Deviation') cint=maxx/6; conts=[cint:cint:maxx]; % contour interval for standard deviation % maxx=2.0 sstd=reshape(sstd, lonx,latx); sste=reshape(sste, lonx,latx); sstdp=cat(1,sstd,sstd(1,:)); sstep=cat(1,sste,sste(1,:)); maxx=max(max(sstep)) disp('Max standard error') maxx=0.6; cint=0.15; conte=[cint:cint:maxx]; % contour interval for standard error maxe=maxx*1.1; pslat1=-80.; pslat2=80.0; sstdp=sm121(sstdp,lonp,lat,nsmooth); sstep=sm121(sstep,lonp,lat,nsmooth); % stderp= sstdp; % here I save the putative standard error for later division into the difference iplothis=1; % plotting Stanard Error if iplothis==1 figure; % Try doing a polar stereographic projection ax = axesm('MapProjection','robinson','Grid','on','Frame','on','FLineWidth',1.0,'MapLatLimit',[pslat1 pslat2],'MapLonLimit',[30 390]); ax = axesm('MapProjection','robinson','Grid','on','Frame','on','FLineWidth',1.0,'MapLatLimit',[pslat1 pslat2],'MapLonLimit',[30 390]); %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,sstep',conte,'linewidth',0.5); caxis([0 maxe]); %,'LineWidth',0.5,'Fill','on','LineColor','black'); %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','black') set(gca,'FontSize',12); colorbar('Location','SouthOutside'); ax = axesm('MapProjection','robinson','Grid','on','Frame','on','FLineWidth',1.0,'MapLatLimit',[pslat1 pslat2],'MapLonLimit',[30 390]); title(['STD Error SST ',var1,num2str(yrbeg),'-',' Season ',num2str(im1),'-',num2str(im2),', ',num2str(yrend,4),'Filt=',num2str(jfilt),num2str(icor),' Lag=',num2str(lag)]); geoshow(coast.lat,coast.long,'DisplayType','polygon','FaceColor','black') set(gca,'Color','none') hold off saveas(gca,['SSTSER_G',var1,num2str(yrbeg),'_',num2str(jfilt),'_',num2str(yrend,4),'_',num2str(im1),'-',num2str(im2),'-',num2str(icor),num2str(gmremove),num2str(lag),'.epsc']); export_fig ['SSTSER_G',var1,num2str(yrbeg),'_',num2str(jfilt),'_',num2str(yrend,4),'_',num2str(im1),'-',num2str(im2),'-',num2str(icor),num2str(gmremove),num2str(lag),'.epsc'] -transparent; iplotthis=0 if iplotthis==1 FH = figure; % Plot a second time because thefirst one is always FU ax = axesm('MapProjection','robinson','Grid','on','Frame','on','FLineWidth',1.0,'MapLatLimit',[pslat1 pslat2],'MapLonLimit',[30 390]); %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); maxe=conts(length(conts)); caxis([0 maxe]); %,'LineWidth',0.5,'Fill','on','LineColor','black'); %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','black') set(gca,'FontSize',12); colorbar('Location','SouthOutside'); ax = axesm('MapProjection','robinson','Grid','on','Frame','on','FLineWidth',1.0,'MapLatLimit',[pslat1 pslat2],'MapLonLimit',[30 390]); title(['STD SST ',var1,num2str(yrbeg),'-',' Season ',num2str(im1),'-',num2str(im2),', ',num2str(yrend,4),'Filt=',num2str(jfilt),num2str(icor),' Lag=',num2str(lag)]); geoshow(coast.lat,coast.long,'DisplayType','polygon','FaceColor','black') set(gca,'Color','none') hold off saveas(gca,['SSTSIG_G',var1,num2str(yrbeg),'_',num2str(jfilt),'_',num2str(yrend,4),'_',num2str(im1),'-',num2str(im2),'-',num2str(icor),num2str(gmremove),num2str(lag),'.epsc']); export_fig ['SSTSIG_G',var1,num2str(yrbeg),'_',num2str(jfilt),'_',num2str(yrend,4),'_',num2str(im1),'-',num2str(im2),'-',num2str(icor),num2str(gmremove),num2str(lag),'.epsc'] -transparent; end % iplotthis if end end % iplothis iff end % now plot the SAM Composite difference ************************************* maxx=squeeze(max(max(abs(xmp)))) cint=maxx/10; % maxx=2.0 sstdp=cat(1,xmp,xmp(1,:)); lon=double(lon); dlon=lon(2)-lon(1); lonp=cat(1,lon,lon(lonx)+dlon); lat=double(lat); cint=maxx/10.; if var1=='sst' cint=0.1; maxx=1.5; elseif var1=='u10' if lag==0 cint=1.0; maxx=5.; else cint=1.0; maxx=5.; end end conts=[-maxx:cint:-cint,cint:cint:maxx]; % conts=[-maxx:cint:maxx];p maxe=maxx; if var0=='nino' conts=[-4. -3. -2. -1. -0.5 -0.25 0.25 0.5 1. 2. 3. 4.]; maxe=5 end sstdp=sm121(sstdp,lonp,lat,nsmooth); pslat1=-80.; pslat2=80.0; FH = figure; % Try doing a polar stereographic projection ax = axesm('MapProjection','robinson','Grid','on','Frame','on','FLineWidth',1.0,'MapLatLimit',[pslat1 pslat2],'MapLonLimit',[30 390]); [cmap]=buildcmap('bwr'); colormap(cmap); [C, H] = contourfm(lat,lonp,sstdp',conts,'linewidth',0.5); caxis([-maxe maxe]); hold on geoshow(coast.lat,coast.long,'DisplayType','polygon','FaceColor','black') set(gca,'FontSize',12); colorbar('Location','SouthOutside'); ax = axesm('MapProjection','robinson','Grid','on','Frame','on','FLineWidth',1.0,'MapLatLimit',[pslat1 pslat2],'MapLonLimit',[30 390]); title(['SST SAM ',var1,'-',var0,num2str(yrbeg),'-',num2str(yrend),'-',num2str(im1),'-',num2str(im2),'-',num2str(ensoremove),' L',num2str(lag)]); geoshow(coast.lat,coast.long,'DisplayType','polygon','FaceColor','black') set(gca,'Color','none') hold off saveas(gca,['SSTSAM_G',var1,var0,num2str(yrbeg),'_',num2str(jfilt),'_',num2str(yrend,4),'_',num2str(im1),'-',num2str(im2),'-',num2str(ensoremove),num2str(gmremove),num2str(lag),'.epsc']); export_fig ['SSTSAM_G',var1,var0,num2str(yrbeg),'_',num2str(jfilt),'_',num2str(yrend,4),'_',num2str(im1),'-',num2str(im2),'-',num2str(ensoremove),num2str(gmremove),num2str(lag),'.epsc'] -transparent; % ok let's try to do statistic significance right conts=[-maxx:cint:-cint,cint:cint:maxx]; sstdp = sstdp./sstep; % here I take the ratio of the change to the standard error kind of a Z statistic maxx=max(max(abs(sstdp))) disp('Max of T-statistic') if var1=='sst' maxx=7.2; cint=1.96; elseif var1=='u10' if lag==0 maxx=15; cint=1.96; else maxx=15; cint=1.96; end end if var0=='nino' maxx=15; cint=1.96; end %maxx=cint*4.; maxe=maxx; %conts=[-maxx:cint:cint,cint:cint:maxx]; conts=[-3.18,-2.33,-1.96,1.96,2.33,3.18]; % these are 0.05, 0.01 and 0.001 probability levels. pslat1=-80.; pslat2=80.0; FH = figure; % Try doing a polar stereographic projection ax = axesm('MapProjection','robinson','Grid','on','Frame','on','FLineWidth',1.0,'MapLatLimit',[pslat1 pslat2],'MapLonLimit',[30 390]); [cmap]=buildcmap('bwr'); colormap(cmap); [C, H] = contourfm(lat,lonp,sstdp',conts,'linewidth',0.5); caxis([-maxe maxe]); hold on geoshow(coast.lat,coast.long,'DisplayType','polygon','FaceColor','black') set(gca,'FontSize',12); colorbar('Location','SouthOutside'); ax = axesm('MapProjection','robinson','Grid','on','Frame','on','FLineWidth',1.0,'MapLatLimit',[pslat1 pslat2],'MapLonLimit',[30 390]); title(['SST P<0.01 ',var1,'-',var0,num2str(yrbeg),'-',num2str(im1),'-',num2str(im2),'-',num2str(yrend,4),' L',num2str(lag)]); geoshow(coast.lat,coast.long,'DisplayType','polygon','FaceColor','black') set(gca,'Color','none') hold off saveas(gca,['SSTSAM_Pstat',var1,var0,num2str(yrbeg),'_',num2str(jfilt),'_',num2str(yrend,4),'_',num2str(im1),'-',num2str(im2),'-',num2str(ensoremove),num2str(gmremove),num2str(lag),'.epsc']); export_fig ['SSTSAM_Pstat',var1,var0,num2str(yrbeg),'_',num2str(jfilt),'_',num2str(yrend,4),'_',num2str(im1),'-',num2str(im2),'-',num2str(ensoremove),num2str(gmremove),num2str(lag),'.epsc'] -transparent; pause 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',[30 390]); %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','black'); %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','black') set(gca,'FontSize',12); colorbar('Location','SouthOutside'); ax = axesm('MapProjection','robinson','Grid','on','Frame','on','FLineWidth',1.0,'MapLatLimit',[pslat1 pslat2],'MapLonLimit',[30 390]); title ['STD SST ',var1,var0,num2str(yrbeg),'-',num2str(im1),'-',num2str(im2),', ',num2str(yrbeg,4),'Filt=',num2str(jfilt),num2str(icor),' Lag=',num2str(lag)]); geoshow(coast.lat,coast.long,'DisplayType','polygon','FaceColor','black') set(gca,'Color','none') hold off saveas(gca,['SSTSIG_G',var1,var0,num2str(yrbeg),'_',num2str(jfilt),'_',num2str(yrbeg,4),'_',num2str(im1),'-',num2str(im2),'-',num2str(icor),num2str(gmremove),num2str(lag),'.epsc']); export_fig ['SSTSIG_G',var1,var0,num2str(yrbeg),'_',num2str(jfilt),'_',num2str(yrbeg,4),'_',num2str(im1),'-',num2str(im2),'-',num2str(icor),num2str(gmremove),num2str(lag),'.epsc'] -transparent; % ******* ******** ************* ********* ********* 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=im2-im1; if im2 > 12 yrmx=yrmx-1; end im1=1; im2=im2-im1+1; imo1=im1; imo2=im2; 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 if topoz(i,j) < 1.0e-5 ss(ind,:)=sst(i,j,:); mask(i,j)=1.0; ind=ind+1; end end end ind if lag == 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 1900 to be consisten as El Nino if yrbeg == 1900 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_',expname,num2str(xlat1),'_',num2str(jfilt),'_',num2str(yrbeg,4),'_',num2str(im1),'-',num2str(im2),'-',num2str(icor),num2str(gmremove),num2str(lag),'.epsc']); export_fig ['EigenSpec_',expname,num2str(xlat1),'_',num2str(jfilt),'_',num2str(yrbeg,4),'_',num2str(im1),'-',num2str(im2),'-',num2str(icor),num2str(gmremove),num2str(lag),'.epsc'] -transparent % OK, time for regression % loop on eigenvalues %year=([1:timxz])/12. + 1900.; for iv = 1:ivmx % loop on eigenmodes *** ***** **** *** *** *** pc=V(:,iv); pc=(pc-mean(pc))/std(pc); if yrbeg == 1900 if xlat2 == -70. pc=-pc; end end if yrbeg == 1900 if xlat2 == -30. if iv < 3 pc=-pc; end end end if yrbeg == 1900 if xlat2 == -30. if iv == 3 pc=-pc; end end end if yrbeg == 1900 if xlat2 == -30 if iv ==3 pc=-pc; end end end if yrbeg == 1880 if xlat2 == -20 if iv == 1 pc=-pc; end end end if yrbeg == 1900 if xlat2 == -20 if iv < 5 pc=-pc; end end end if yrbeg == 1880 if xlat2 == 20 if iv == 2 pc=-pc; end end end if yrbeg == 1880 if detrend1 == 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 == 1900 pc=pc; end end if iv == 2 if yrbeg == 1900 pc=-pc; end end if yrbeg == 1900 pc=-pc; end if xlat1 == 30. if iv < 3 pc=-pc; end end if xlat2 == -70. if iv == 2 pc=-pc; end end if lag == 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(['G_1900_V(:,1)',' PC#',num2str(iv)]) set(gca,'FontSize',14); ylabel('PC Amplitude'); xlabel('Year') axis([yrbeg yrend -3 4]); pbaspect([2 1 1]); set(gca,'Color','none') saveas(gca,['PC_G_',var1,num2str(yrbeg),'_TimeSeries_',num2str(iv),'_',var1,num2str(yrend),'.epsc']); export_fig ['PC_G_',var1,num2str(yrbeg),'_TimeSeries_',num2str(iv),'_',var1,num2str(yrend),'.epsc'] -transparent figure ykmxm=ykmx-72; plot(years(ykmxm:ykmx5),pc(ykmxm:ykmx5)); title(['G_',var1,num2str(yrbeg),'_V(:,1)',' PC#',num2str(iv)]) set(gca,'FontSize',20); ylabel('PC Amplitude'); xlabel('Year') axis([2010 2016 -3 4]); pbaspect([2 1 1]); end % end of plot if lag if lag == 0 pc=pc'; fileout=['PDO_',expname,'_',var1,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(yrbeg,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,['PDO_',expname,num2str(yrbeg,4),'_Power_',num2str(iv),'_',num2str(m),'_',var1,num2str(yrbeg),'.epsc']) export_fig ['PDO_',expname,num2str(yrbeg,4),'_Power_',num2str(iv),'_',num2str(m),'_',var1,num2str(yrbeg),'.epsc'] -transparent end % of plot lag 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 lag == 0 regr=sstr*pc'/NT; elseif lag > 0 clear sstrs; clear pcs; NTM = NT -lag; sstrs(:,1:NTM)=sstr(:,lag+1:NT); pcs(1:NTM)=pc(1:NTM); regr=sstrs*pcs'/NTM; elseif lag < 0 clear sstrs; clear pcs; jlag = abs(lag) 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,:)); % Marc, % The code that makes the figures is below. I think you can save either the data % or the located data on the plot. Maybe just saving the data is easier. % I guess you just need the lat, long and value % Thanks for taking a look at this. 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',[30 390]); %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','black'); %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','black') set(gca,'FontSize',12); ax = axesm('MapProjection','robinson','Grid','on','Frame','on','FLineWidth',1.0,'MapLatLimit',[pslat1 pslat2],'MapLonLimit',[30 390]); title(['PDO Regression',expname,', Season ',num2str(im1),'-',num2str(im2),' mode=',num2str(iv),', ',num2str(yrbeg,4),'Filt=',num2str(jfilt),num2str(icor),num2str(gmremove),' Lag=',num2str(lag)]); geoshow(coast.lat,coast.long,'DisplayType','polygon','FaceColor','black') set(gca,'Color','none') hold off saveas(gca,['PDO_',expname,'_',num2str(iv),'_',num2str(jfilt),'_',num2str(yrbeg,4),'_',num2str(im1),'-',num2str(im2),'-',num2str(icor),num2str(gmremove),num2str(lag),'.epsc']); export_fig ['PDO_',expname,'_',num2str(iv),'_',num2str(jfilt),'_',num2str(yrbeg,4),'_',num2str(im1),'-',num2str(im2),'-',num2str(icor),num2str(gmremove),num2str(lag),'.epsc'] -transparent R = georasterref('RasterSize',size(regp'),'Latitudelimits',[-90 90],'LongitudeLimits',[0 360]); geotiffwrite(['PDO_Atlantic_65N-65S_1900_',num2str(iv),'_',num2str(jfilt),'_',num2str(yrbeg,4),'_',num2str(im1),'-',num2str(im2),'-',num2str(icor),num2str(gmremove),num2str(lag),'.tif'],regp',R); end % end of eigenvalue loop