load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl" begin y1 = 1951 ; year start y2 = 2016 ; year end nsvd = 2 ; number of modes season = "OND" infile1 = "cru_ts4.03.1901.2018.pre.dat.nc" ; data for variable 1(Z500) infile2 = "ersst.mnmean.nc" ; data for variable 2(SST) out = "eps" ; output form ;;; Read the data of variable 1(Z500) in1 = addfile(infile1,"r") ; get file to read time = in1->time ; get time array YYYY = cd_calendar(time,-1)/100 ; Time in YYYYMM format (option -1). ; The division by 100 shows only the years in YYYY format. iYYYY = ind(YYYY.ge.y1 .and. YYYY.le.y2) ; Select indices between 1979 and 2003. pre = in1->pre(iYYYY,:,:) ; This corresponds to cropping only the values ​​from January 1951 to 2016 ;********************************************************************************************************************************* ;create low pass filter ;********************************************************************************************************************************* ihp =0 ;low pass filter sigma =1.0 ;Lanczos sigma nWgt =11 ;loose 120 months at each end fca =1./10. ;cut-off frequency, indicating decadal signal wgtp =filwgts_lanczos(nWgt,ihp,fca,-999.,sigma) ; ============================================================== ; detrend pre ; ============================================================== pre = dtrend_msg_n(ispan(0,dimsizes(pre&time)-1,1),pre,True,False,0) ; average and trend removal clima=clmMonTLL(pre) preanom=calcMonAnomTLL(pre,clima) ; remove the time average to the deviation printVarSummary(pre) ;exit ; ============================================================== ; compute desired global seasonal mean: month_to_season (contributed.ncl) ; ============================================================== pre_sm = month_to_season (preanom, season) ; Extracts only the average of OND. ; ================================================================= ; Apply Lanczos filter ; ============================================================== pre_filter = wgt_runave_n_Wrap (pre_sm, wgtp, 0,0) lat1 = in1->lat ; get an array of latitudes corresponding to Z500 lon1 = in1->lon ; get an array of longitude array corresponding to Z500 pre_filter!0 = "time" ; add coordinates to pre pre_filter!1 = "lat" pre_filter!2 = "lon" delete([/time,YYYY/]) ; delete time and YYYYMM array ; (because the same variable name is used when reading SST data) printVarSummary(pre_filter) ;exit ;;; Reading variable 2(SST) data as above in2 = addfile(infile2,"r") time = in2->time YYYY = cd_calendar(time,-1)/100 iYYYY = ind(YYYY.ge.y1 .and. YYYY.le.y2) sst = in2->sst(iYYYY,:,:) ; ============================================================== ; detrend sst ; ============================================================== sst = dtrend_msg_n(ispan(0,dimsizes(sst&time)-1,1),sst,True,False,0) climo=clmMonTLL(sst) sstanom=calcMonAnomTLL(sst,climo) printVarSummary(sst) ;exit ; ============================================================== ; compute desired global seasonal mean: month_to_season (contributed.ncl) ; ============================================================== sst_sm = month_to_season (sstanom, season) ; ================================================================= ; Apply Lanczos filter ; ============================================================== sst_filter = wgt_runave_n_Wrap (sst_sm, wgtp, 0,0) lat2 = in2->lat lon2 = in2->lon sst_filter!0 = "time" sst_filter!1 = "lat" sst_filter!2 = "lon" delete([/time,YYYY/]) printVarSummary(sst_filter) ;exit ; ============================================================== ;;; Prepare SVD analysis ; ============================================================== ;;; First, for variable 1(Z500) LAT1 = ind(lat1.ge.-12.and.lat1.lt.1) ; specify the SVD area (latitude range) of variable 1 LON1 = ind(lon1.ge.28.and.lon1.le.42) ; specify the SVD area (longitude range) of variable 1 var1 = pre_filter(lat|LAT1,lon|LON1,time|:) ; cut out necessary area of Z500 ; however, the order of dimension of time and space is changed dim1 = dimsizes(var1) ; the size of each dimension is NY1,NX1,NT1 from the left NY1 = dim1(0) NX1 = dim1(1) NT1 = dim1(2) ;;; weight each grid area ;;; in SVD analysis, the area weight should be applied to each component of the covariance matrix. ;;; Note that the square root of the area is applied to the deviation wgt1 = sqrt(cos(lat1(LAT1)*atan(1.)/45.)) ; vcos(latitude) weight(latitude width is uniform) NW1 = dimsizes(wgt1) do j = 0, NW1-1 var1(j,:,:) = var1(j,:,:) * wgt1(j) end do ;;; change the shape of var1. Combine the dimensions of latitude and longitude into a two-dimensional array of space and time var1 := reshape(var1,(/NY1*NX1,NT1/)) ;;; SVD functions do not work well if they contain undefined values. Therefore, it is necessary to remove undefined values. ;;; Now, assume that undefined values are included in the same spatial position regardless of time (like the land of SST data) NMS1 = ind(.not.ismissing(var1(:,0))) ; get non-defined location in spatial dimension printVarSummary(NMS1) var1 := var1(NMS1,:) ; cut out undefined values and make the array smaller ; NMS1 defined here is necessary to restore the original size later printVarSummary(var1) ;exit ;;; Similarly,prepare for variable 2(SST) LAT2 = ind(lat2.ge.-20.and.lat2.le.20) LON2 = ind(lon2.ge.40.and.lon2.le.120) var2 = sst_filter(lat|LAT2,lon|LON2,time|:) dim2 = dimsizes(var2) NY2 = dim2(0) NX2 = dim2(1) NT2 = dim2(2) ; must be the same as NT1 wgt2 = sqrt(cos(lat2(LAT2)*atan(1.)/45.)) NW2 = dimsizes(wgt2) do j = 0, NW2-1 var2(j,:,:) = var2(j,:,:) * wgt2(j) end do var2 := reshape(var2,(/NY2*NX2,NT2/)) NMS2 = ind(.not.ismissing(var2(:,0))) var2 := var2(NMS2,:) ;;; Run SVD.This time, both singular vectors and homogeneous/heterogeneous correlation maps are obtained ;;; First, calculate singular vectors using svdcov_sv var1_svc = new((/nsvd,dimsizes(NMS1)/),float) ; singular vector of variable 1(Z500) var2_svc = new((/nsvd,dimsizes(NMS2)/),float) ; singular vectors of variable 2(SST) pcvar_a = svdcov_sv(var1,var2,nsvd,var1_svc,var2_svc) ; SVD var1_hom = new((/nsvd,dimsizes(NMS1)/),float) ; homogenous correlation map of variable 1(Z500) var1_het = new((/nsvd,dimsizes(NMS1)/),float) ; homogeneous correlation map of variable 2(SST) var2_hom = new((/nsvd,dimsizes(NMS2)/),float) ; heterogrneous correlation map of variable 1(Z500) var2_het = new((/nsvd,dimsizes(NMS2)/),float) ; heterogeneous correlation map of varuable 2(SST) pcvar_b = svdcov(var1,var2,nsvd,var1_hom,var1_het,var2_hom,var2_het) ; SVD ;;; pcvar_a and pcvar_b should have the same values because they are doing the same SVD (attributes are different) ;;; pcvar_b has a time function of variable 1 as @ak and a time funstion of variable 2 as @bk. ;;; Note that the time function is a one-dimensional array of nsvd×time length. pre_pc = reshape(pcvar_b@ak,(/nsvd,NT1/)) ; time function of variable 1(Z500) sst_pc = reshape(pcvar_b@bk,(/nsvd,NT2/)) ; time function of variable 2(SST) ;;; By the way, since the time function is the projection of the original data onto the onto the singular vector, it can be calculated as follows. ;;; pre_pc = var1_svc#var1 ;;; sst_pc = var2_svc#var2 ;;; However, since the time functions of different modes are not orthogonal in SVD, ;;; Note that the singular vector cannot be obtained by projecting the original data onto the time function. pre_pc@pcvar = pcvar_b ; Co-dispersion rate sst_pc@pcvar = pcvar_b pre_svc = new((/nsvd,NY1*NX1/),float) ; Z500 singular vectors(space dimension is defined by the size including undefined values) ; By default, all elements have undefined values pre_svc(:,NMS1) = var1_svc ; Substitute the singular vector obtained by SVD to the place that was not an undefined value pre_svc := reshape(pre_svc,(/nsvd,NY1,NX1/)) ; Return spatial dimension to two dimensions of latitude×longitude sst_svc = new((/nsvd,NY2*NX2/),float) ; SST singular vectors sst_svc(:,NMS2) = var2_svc sst_svc := reshape(sst_svc,(/nsvd,NY2,NX2/)) pre_hom = new((/nsvd,NY1*NX1/),float) ; Z500 homogeneous correlation map pre_hom(:,NMS1) = var1_hom pre_hom := reshape(pre_hom,(/nsvd,NY1,NX1/)) sst_hom = new((/nsvd,NY2*NX2/),float) ; SST homogeneous correlation map sst_hom(:,NMS2) = var2_hom sst_hom := reshape(sst_hom,(/nsvd,NY2,NX2/)) pre_het = new((/nsvd,NY1*NX1/),float) ; Z500 heterogeneous correlation map pre_het(:,NMS1) = var1_het pre_het := reshape(pre_het,(/nsvd,NY1,NX1/)) sst_het = new((/nsvd,NY2*NX2/),float) ; SST heterogeneous correlation map sst_het(:,NMS2) = var2_het sst_het := reshape(sst_het,(/nsvd,NY2,NX2/)) delete([/var1,var2,var1_svc,var2_svc,var1_hom,var1_het,var2_hom,var2_het/]) ;;; I will return it by the amount of the weight do j = 0, NW1-1 pre_svc(:,j,:) = pre_svc(:,j,:) / wgt1(j) end do do j = 0, NW2-1 sst_svc(:,j,:) = sst_svc(:,j,:) / wgt2(j) end do ;;; Set coordinates for each variable in preparation for drawing pre_svc!0 = "svd" pre_svc!1 = "lat_pre" pre_svc!2 = "lon_pre" pre_svc&svd = ispan(1,nsvd,1) pre_svc&lat_pre = lat1(LAT1) pre_svc&lon_pre = lon1(LON1) copy_VarCoords(pre_svc,pre_hom) copy_VarCoords(pre_svc,pre_het) pre_pc!0 = "svd" pre_pc!1 = "year" pre_pc&svd = ispan(1,nsvd,1) pre_pc&year = ispan(y1,y2,1) sst_svc!0 = "svd" sst_svc!1 = "lat_sst" sst_svc!2 = "lon_sst" sst_svc&svd = ispan(1,nsvd,1) sst_svc&lat_sst = lat2(LAT2) sst_svc&lon_sst = lon2(LON2) copy_VarCoords(sst_svc,sst_hom) copy_VarCoords(sst_svc,sst_het) sst_pc!0 = "svd" sst_pc!1 = "year" sst_pc&svd = ispan(1,nsvd,1) sst_pc&year = ispan(y1,y2,1) ;;; Although not always necessary, normalize the time function and convert the singular vectors into quantities with dimensions. do i = 0, nsvd-1 pre_svc(i,:,:) = pre_svc(i,:,:)*stddev(pre_pc(i,:)) sst_svc(i,:,:) = sst_svc(i,:,:)*stddev(sst_pc(i,:)) pre_pc(i,:) = pre_pc(i,:)/stddev(pre_pc(i,:)) sst_pc(i,:) = sst_pc(i,:)/stddev(sst_pc(i,:)) end do printMinMax(pre_svc,0) printMinMax(sst_svc,0) printMinMax(pre_pc,0) printMinMax(sst_pc,0) ;exit ;;; Plot setting res1 = True ; resourcer used to draw Z500 res1@gsnDraw = False ; do not draw when executing gsn_csm_contour_map_ce res1@gsnFrame = False ; do not update Frame when gsn_csm_contour_map_ce is xecuted res1@cnFillOn = True ; use shade res1@cnLinesOn = False ; do not use contours res1@cnLineLabelsOn = False ; no contour label res1@cnInfoLabelOn = False ; no contour information res1@lbLabelBarOn = False ; do not draw a color bar for each panel (for now) ; later temporarily becomes True res1@lbTitleOn = True ; add a title to the color bar res1@lbTitlePosition = "Right" ; color bar title position res1@lbTitleDirection = "Across" ; color bar title direction res1@pmLabelBarOrthogonalPosF = 0.3 ; move the color bar slightly down res1@mpGeophysicalLineThicknessF = 1. ; coastline thickness res1@gsnAddCyclic = False ; not periodic in the longitude direction res1@tmXTOn = False ; do not use the upper x-axis (turn off the scale) res1@cnLevelSelectionMode = "ManualLevels" ; set contour level to "ManualLevels" mode ;;; Actually, this is the same as the resource used foe SST drawing. ;;; So copy the resource used for SST drawing as res3 res3 = res1 ;;; After that, we gonna make individual settings ;;; Against Z500 res1@mpMinLatF = -12. ; southern edge of map res1@mpMaxLatF = 1. ; northern edge of map res1@mpMinLonF = 28. ; western edge of map res1@mpMaxLonF = 42. ; eastern edge of map ;res1@mpCenterLonF = 180. ; has longitude coordinates from 0° to 360° res1@cnMinLevelValF = -50. ; contour min value res1@cnMaxLevelValF = 50. ; contour max value res1@cnLevelSpacingF = 10. ; contour interval res1@lbTitleString = "mm" ; color bar title res1@lbTitleFontHeightF = 0.021 ; font size of color bar res1@lbLabelFontHeightF = 0.021 ; font size of color label res1@gsnLeftString = "PRE" ; LeftString title res1@gsnLeftStringFontColor = "red" ; LeftString font color res1@gsnLeftStringFontHeightF = 0.038 ; LeftString font height res1@tmYLOn = True ; use the left y-axis res1@tmYLLabelsOn = True ; use the left y-axis for vertical axis label res1@tmYROn = False ; dont use right y-axis res1@tmYRLabelsOn = False ; the vertical axis label does not use the right y-axis res1@tmXBLabelFontHeightF = 0.025 ; horizontal axis label text size res1@tmYLLabelFontHeightF = 0.025 ; size of vertical axis label text res1@mpDataBaseVersion = "MediumRes" ;res1@mpDataSetName = "Earth..4" res1@mpOutlineOn = True ; Use outlines res1@cnFillDrawOrder = "preDraw" res1@cnLineDrawOrder = "preDraw" res1@mpAreaMaskingOn = True res1@mpOutlineOn = True res1@mpMaskAreaSpecifiers = (/"tanzania"/) res1@mpLandFillColor = "gray" res1@mpInlandWaterFillColor= "white" res1@mpOceanFillColor = "white" res1@mpOutlineBoundarySets = "National" ;;; For SST(basically the same, but put different values) res3@mpMinLatF = -20. res3@mpMaxLatF = 20. res3@mpMinLonF = 40. res3@mpMaxLonF = 120. ;res3@mpCenterLonF = 180. res3@cnMinLevelValF = -0.3 res3@cnMaxLevelValF = 0.3 res3@cnLevelSpacingF = 0.1 res3@lbTitleString = "~S~o~N~C" res3@lbTitleFontHeightF = 0.021 res3@lbLabelFontHeightF = 0.021 res3@gsnLeftString = "SST" res3@gsnLeftStringFontColor = "blue" res3@gsnLeftStringFontHeightF = 0.038 res3@tmYLOn = False ; dont use left y-axis res3@tmYLLabelsOn = False res3@tmYROn = True ; use the right y-axis res3@tmYRLabelsOn = True res3@tmXBTickSpacingF = 20 res3@tmXBLabelFontHeightF = 0.025 res3@tmYLLabelFontHeightF = 0.025 res3@mpShapeMode = "FreeAspect" ; variable map plot aspect ratio res3@vpWidthF = 0.6 ; figure width (horizontal length changes when attached) ; the map becomes taller along with this res3@tmXBTickSpacingF = 20 ; the horizontal scale interval is 45° res2 = True ; resource used for drawing time funstions res2@gsnDraw = False ; do not draw when executing gsn_csm_xy res2@gsnFrame = False ; do not update Frame when executing gsn_csm_xy res2@trXMinF = y1 ; horizontal axis range res2@trXMaxF = y2 ; res2@trYMinF = -4 ; vertical axis range res2@trYMaxF = 4 ; res2@tmXTOn = False ; do not use the upper x-axis(turn off the scale) res2@vpHeightF = 0.3 res2@tmYRLabelsOn = False ; do not draw labels on the right y-axis res2@xyLineColors = (/"red","blue"/) ; line color res2@xyDashPattern = 0 ; line type of solid line (solid line) res2@gsnYRefLine = 0 ; draw a line with y=0 res2@tmXBLabelFontHeightF = 0.045 ; horizontal axis label text size res2@gsnLeftStringFontHeightF = 0.045 ; LeftString font height pres = True ; panel plot resource pres@gsnPanelTop = 0.9 ; make a margin in the top 10% of the panel plot pres@gsnPanelBottom = 0.1 ; make a margin at the bottom 10% of the panel plot pres@gsnPanelYWhiteSpacePercent = 10 ; insert 10% margin between each panel pres@gsnPanelLabelBar = False ; do not draw color bars (for now) ; become True later pres@pmLabelBarWidthF = 0.6 ; color bar width pres@pmLabelBarHeightF = 0.05 ; color bar height pres@lbLabelFontHeightF = 0.012 ; color bar font height ;;; attach resources attach_res1 = True attach_res2 = True attach_res3 = True ;;; First:time and singular function wks = gsn_open_wks(out,"svd_sv") plot1 = new(nsvd,graphic) plot2 = new(nsvd,graphic) plot3 = new(nsvd,graphic) dum = new(nsvd,graphic) ;;; drawing for each mode do i = 0, nsvd-1 ;; draw color bars for,Z500 and SST only in the last mode if (i.eq.nsvd-1) then res1@lbLabelBarOn = True res3@lbLabelBarOn = True end if ;; by attaching to the title and,res2 of each panel, it will be located near the center ;; " ~C~ " is a line feed res2@gsnLeftString = "SVD "+sprinti("%1.1i",pre_svc&svd(i)) + \ " : "+sprintf("%3.1f",pre_pc@pcvar(i)) + "% ~C~ " ;; draw a singular vector of Z500 plot1(i) = gsn_csm_contour_map_ce(wks,pre_svc(i,:,:),res1) ;; draw time function plot2(i) = gsn_csm_xy(wks,pre_pc&year,(/pre_pc(i,:),sst_pc(i,:)/),res2) ;; draw a singular vector of SST plot3(i) = gsn_csm_contour_map_ce(wks,sst_svc(i,:,:),res3) ;; attach SST singular vector graph to time function graph dum(i) = gsn_attach_plots(plot2(i),plot3(i),attach_res2,attach_res3) ;; then attach it to the Z500 singular vector graph dum(i) = gsn_attach_plots(plot1(i),plot2(i),attach_res1,attach_res2) end do gsn_panel(wks,plot1,(/nsvd,1/),pres) ; panel plot delete([/wks,plot1,plot2,plot3,dum/]) ; delete for next drawing ;;; Particular changes in preparation for drawing homogeneous/heterogeneous correlation map res1@lbLabelBarOn = False ; do not draw color bars foe individual panels res3@lbLabelBarOn = False ; pres@gsnPanelLabelBar = True ; instead, attach a color bar to the whole res1@cnMinLevelValF = -0.9 ; change the maximum/minimum/interval of contour res1@cnMaxLevelValF = 0.9 ; correlation coefficient is -1 to +1 res1@cnLevelSpacingF = 0.1 res3@cnMinLevelValF = -0.9 res3@cnMaxLevelValF = 0.9 res3@cnLevelSpacingF = 0.1 ;;; Second: time function and homogeneous correlation map wks = gsn_open_wks(out,"svd_hom") plot1 = new(nsvd,graphic) plot2 = new(nsvd,graphic) plot3 = new(nsvd,graphic) dum = new(nsvd,graphic) do i = 0, nsvd-1 res2@gsnLeftString = "SVD "+sprinti("%1.1i",pre_hom&svd(i)) + \ " : "+sprintf("%3.1f",pre_pc@pcvar(i)) + "% ~C~ " plot1(i) = gsn_csm_contour_map_ce(wks,pre_hom(i,:,:),res1) plot2(i) = gsn_csm_xy(wks,pre_pc&year,(/pre_pc(i,:),sst_pc(i,:)/),res2) plot3(i) = gsn_csm_contour_map_ce(wks,sst_hom(i,:,:),res3) dum(i) = gsn_attach_plots(plot2(i),plot3(i),attach_res2,attach_res3) dum(i) = gsn_attach_plots(plot1(i),plot2(i),attach_res1,attach_res2) end do gsn_panel(wks,plot1,(/nsvd,1/),pres) delete([/wks,plot1,plot2,plot3,dum/]) ;;; Third:time function and heterogeneous correlation map wks = gsn_open_wks(out,"svd_het") plot1 = new(nsvd,graphic) plot2 = new(nsvd,graphic) plot3 = new(nsvd,graphic) dum = new(nsvd,graphic) do i = 0, nsvd-1 res2@gsnLeftString = "SVD "+sprinti("%1.1i",pre_het&svd(i)) + \ " : "+sprintf("%3.1f",pre_pc@pcvar(i)) + "% ~C~ " plot1(i) = gsn_csm_contour_map_ce(wks,pre_het(i,:,:),res1) plot2(i) = gsn_csm_xy(wks,pre_pc&year,(/pre_pc(i,:),sst_pc(i,:)/),res2) plot3(i) = gsn_csm_contour_map_ce(wks,sst_het(i,:,:),res3) dum(i) = gsn_attach_plots(plot2(i),plot3(i),attach_res2,attach_res3) dum(i) = gsn_attach_plots(plot1(i),plot2(i),attach_res1,attach_res2) end do gsn_panel(wks,plot1,(/nsvd,1/),pres) delete([/wks,plot1,plot2,plot3,dum/]) end