; ============================================================== ; eof_SST.ncl ; ; Concepts illustrated: ; - Calculating EOFs ; - Drawing a time series plot ; - Using coordinate subscripting to read a specified geographical region ; - Rearranging longitude data to span -180 to 180 ; - Calculating symmetric contour intervals ; - Drawing filled bars above and below a given reference line ; - Drawing subtitles at the top of a plot ; - Reordering an array ; ============================================================== ; NCL V6.4.0 has new functions eofunc_n_Wrap and ; eofunc_ts_n_Wrap that allow you to calculate the EOFs without ; first having to first reorder the data. See eof_1_640.ncl. ; ============================================================== ; Calculate EOFs of the Sea Surface Temperature over the Tropical Indian Ocean ; ============================================================== ; ============================================================== ; These files are loaded by default in NCL V6.2.0 and newer ; 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" ;============================================================================== ;A function for averaging june july august september undef ("season_from_month") function season_from_month (xMon:numeric, seastart:numeric, seaend:numeric, season:string) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; This function written by Gibies George is adapted from ; D Shea's month_to_seasonN() and month_to_season12() functions, ; with some slight difference in script, but enterly different usage ; which enables objective definision of season. ; ; Compute the seasonal average for user specified seasons. ; ; xMon(time) or xMon(time,lat,lon) or xMon(time,lev,lat,lon) ; ; The input "x" are assumed to contain monthly mean data ; The size of "time" MUST be divisible by 12. ; Also, it is assumed the "Jan" is the 1st month. ; ; first DJF season is a 2-month average (DJF=JF) ; ; USAGE : xSea = season_from_month (xMon, 6, 9, JJAS) ; gives "JJAS : Asian Summer Monsoon Season" ; : xSea = season_from_month (xMon, 10, 12, OND) ; gives "OND : Asian Post-monsoon Season" ; : xSea = season_from_month (xMon, 4, 5, pre-monsoon) ; gives "April-May : Asian Pre-monsoon Season" ; ; RESULT xSea(time/12,lat,lon) or xSea(time/12,lev,lat,lon) ; ; The above would return: ; xSea(time/12,lat,lon) or xSea(time/12,lev,lat,lon) ; ; NOTE: the "time" dimension may have to be altered to the user's desires. ; it may correspond to those associated with the 1st month. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; local dimx, rank, nmos, ntim, nyrs, xSea12, nlat, mlon \ , xSeaN, dName, cv begin isdef=isdefined((/"dimsize","dims","rank"/)) if(isdef(0)) delete(dimsize) end if if(isdef(1)) delete(dims) end if if(isdef(2)) delete(rank) end if delete (isdef) dimsize=dimsizes(xMon) dims=getvardims(xMon) rank = dimsizes(dimsize) if (rank.eq.2 .or. rank.ge.5) then print ("contributed: season_from_month: rank="+rank) print ("----- rank currently not handled -----") end if nmos = 12 ntim = dimsize(0) ;print ("contributed: season_from_month: ntim="+ntim+" nmos="+nmos) modCheck ("season_from_month", ntim, nmos) nyrs = ntim/nmos if (seaend.lt.seastart) then seaend = 12 + seaend end if sealen = seaend - seastart ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; ; xSea12 = month_to_season12 (xMon) ; compute the 12 seasons in month_to_seasonN() is here replaced ; by the following script adapted from month_to_season12(), ; with an additional the variable "sealen", ; which stands for the length of the season ; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; if (rank.eq.1) then ; (time) xSea = xMon ; transfer meta and reorder xSea = runave (xSea ,sealen, 0 ) ; overwrite with seasonal means xSea(0) = (xMon(0) + xMon(1) )*0.5 xSea(ntim-1) = (xMon(ntim-2) + xMon(ntim-1) )*0.5 xSea@long_name = season+" means: "+getLongName(xMon) xSea@season = season xSea12 = xSea($dims(0)$|:) end if if (rank.eq.3) then ; (time,lat,lon) ;xSea = xMon(lat|:,lon|:,time|:) ; transfer meta and reorder xSea = xMon($dims(1)$|:,$dims(2)$|:,$dims(0)$|:) ; transfer meta and reorder xSea = runave (xSea ,sealen, 0 ) ; overwrite with seasonal means xSea(:,:,0) = (xMon(0,:,:) + xMon(1,:,:) )*0.5 xSea(:,:,ntim-1) = (xMon(ntim-2,:,:) + xMon(ntim-1,:,:) )*0.5 xSea@long_name = season+" means: "+getLongName(xMon) xSea@season = season xSea12 = xSea($dims(0)$|:,$dims(1)$|:,$dims(2)$|:) end if if (rank.eq.4) then ; (time,lev,lat,lon) ;xSea = xMon(lev|:,lat|:,lon|:,time|:) ; transfer meta and reorder xSea = xMon($dims(1)$|:,$dims(2)$|:,$dims(3)$|:,$dims(0)$|:) xSea = runave (xMon ,sealen, 0 ) ; overwrite with seasonal means xSea(:,:,:,0) = (xMon(0,:,:,:) + xMon(1,:,:,:) )*0.5 xSea(:,:,:,ntim-1) = (xMon(ntim-2,:,:,:) + xMon(ntim-1,:,:,:) )*0.5 xSea@long_name = season+" means: "+getLongName(xMon) xSea@season = season xSea12 = xSea($dims(0)$|:,$dims(1)$|:,$dims(2)$|:,$dims(3)$|:) end if ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; if (rank.ge.3) then nlat = dimsize(rank-2) mlon = dimsize(rank-1) end if if (rank.eq.1) then ; (time,lat,lon) xSeaN = new ( nyrs, typeof(xSea12), getFillValue(xMon)) ;printVarSummary(xSea12) xSeaN(:) = (/ xSea12(seastart:ntim-1:nmos) /) end if if (rank.eq.3) then ; (time,lat,lon) xSeaN = new ( (/nyrs,dimsize(1),dimsize(2)/), typeof(xSea), \ getFillValue(xMon)) ;printVarSummary(xSea12) xSeaN(:,:,:) = (/ xSea12(seastart:ntim-1:nmos,:,:) /) end if if (rank.eq.4) then ; (time,lat,lon) xSeaN = new ( (/nyrs,dimsize(1),dimsize(2),dimsize(3)/), typeof(xSea), \ getFillValue(xMon)) ;printVarSummary(xSea12) xSeaN(:,:,:,:) = (/ xSea12(seastart:ntim-1:nmos,:,:,:) /) end if if (rank.eq.5) then ; (time,lat,lon) xSeaN = new ( (/nyrs,dimsize(1),dimsize(2),dimsize(3),dimsize(4)/), typeof(xSea), \ getFillValue(xMon)) ;printVarSummary(xSea12) xSeaN(:,:,:,:,:) = (/ xSea12(seastart:ntim-1:nmos,:,:,:,:) /) end if copy_VarAtts (xMon, xSeaN) xSeaN@season = season if (isatt(xMon,"long_name") .or. isatt(xMon,"description") .or. \ isatt(xMon,"standard_name") ) then xSeaN@long_name = season+" Means: "+getLongName(xMon) end if ; copy dimension stuff dimName = xSea12!0 xSeaN!0 = dimName if(iscoord(xSea,dimName)) then cv = xSea12&$dimName$(seastart:ntim-1:nmos) xSeaN&$dimName$ = cv ; possibly override if (isatt(cv,"units") .and. \ (cv@units.eq."YYYYMM" .or. cv@units.eq."YYMM")) then cv = cv/100 cv@units = "YYYY" xSeaN&$dimName$ = cv end if if (isatt(cv,"units") .and. cv@units.eq."YYYYMMDD") then cv = cv/10000 cv@units = "YYYY" xSeaN&$dimName$ = cv end if end if if (rank.gt.1) then do i=1,rank-1 ; copy spatial coords dimName = xSea12!i xSeaN!(i) = dimName if(iscoord(xSea,dimName)) then xSeaN&$dimName$ = xSea12&$dimName$ end if end do end if ;printVarSummary(xSeaN) return (xSeaN) end ; ------------------------------------------------------- ;================================================================= begin ; ============================================================== ; User defined parameters that specify region of globe and ; ============================================================== latS = -30. latN = 30. lonL = 40 lonR = 120. yrStrt = 1981 yrLast = 2019 season = "OND" ; choose oct-nov-dec seasonal mean neof = 3 ; number of EOFs optEOF = True optEOF@jopt = 0 ; This is the default; most commonly used; no need to specify. ;;optEOF@jopt = 1 ; **only** if the correlation EOF is desired optETS = False ; ============================================================== ; Open the file: Read only the user specified period ; ============================================================== f = addfile ("hadsstano_1981-2019.nc", "r") TIME = f->time YYYY = cd_calendar(TIME,-1)/100 ; entire file iYYYY = ind(YYYY.ge.yrStrt .and. YYYY.le.yrLast) sst = f->sst(iYYYY,:,:) printVarSummary(sst) ; variable overview lat = f->latitude lon = f->longitude ; ============================================================== ; compute desired global seasonal mean: month_to_season (contributed.ncl) ; ============================================================== SST = season_from_month (sst, 10, 12, "JJAS") nyrs = dimsizes(SST&time) printVarSummary(SST) ; ================================================================= ; create weights: sqrt(cos(lat)) [or sqrt(gw) ] ; ================================================================= rad = 4.*atan(1.)/180. clat = f->latitude clat = sqrt( cos(rad*clat) ) ; gw for gaussian grid ; ================================================================= ; weight all observations ; ================================================================= wSST = SST ; copy meta data wSST = SST*conform(SST, clat, 1) wSST@long_name = "Wgt: "+wSST@long_name ; ================================================================= ; Reorder (lat,lon,time) the *weighted* input data ; Access the area of interest via coordinate subscripting ; ================================================================= xw = wSST({latitude|latS:latN},{longitude|lonL:lonR},time|:) x = wSST(time|:,{latitude|latS:latN},{longitude|lonL:lonR}) eof = eofunc_Wrap(xw, neof, optEOF) eof_ts = eofunc_ts_Wrap (xw, eof, optETS) printVarSummary( eof ) ; examine EOF variables printVarSummary( eof_ts ) ; ================================================================= ; Normalize time series: Sum spatial weights over the area of used ; ================================================================= dimxw = dimsizes( xw ) mln = dimxw(1) sumWgt = mln*sum( clat({latitude|latS:latN}) ) eof_ts = eof_ts/sumWgt ; ================================================================= ; Extract the YYYYMM from the time coordinate ; associated with eof_ts [same as SST&time] ; ================================================================= yyyymm = cd_calendar(eof_ts&time,-2)/100 ;;yrfrac = yyyymm_to_yyyyfrac(yyyymm, 0.0); not used here ;=================================================================== ;Calculate DMIV ;================================================================== weio = wgt_areaave_Wrap(sst(:, -10:10, 50:70), 1.0, 1.0, 0) ;============================================================ ; PLOTS ;============================================================ wks = gsn_open_wks("png","eof_OND") ; send graphics to PNG file plot = new(neof,graphic) ; create graphic array ; only needed if paneling ; EOF patterns res = True res@gsnDraw = False ; don't draw yet res@gsnFrame = False ; don't advance frame yet res@gsnAddCyclic = False ; plotted dataa are not cyclic res@mpFillOn = False ; turn off map fill res@mpMinLatF = latS ; zoom in on map res@mpMaxLatF = latN res@mpMinLonF = lonL res@mpMaxLonF = lonR res@cnFillOn = True ; turn on color fill res@cnLinesOn = False ; True is default ;res@cnLineLabelsOn = False ; True is default res@cnFillPalette = "MPL_RdylGn" ; set color map res@lbLabelBarOn = False ; turn off individual lb's res@cnMinLevelValF = -0.06 res@cnMaxLevelValF = 0.06 ; set symmetric plot min/max symMinMaxPlt(eof, 16, False, res) ; contributed.ncl ; panel plot only resources resP = True ; modify the panel plot resP@gsnMaximize = True ; large format resP@gsnPanelLabelBar = True ; add common colorbar yStrt = yyyymm(0)/100 yLast = yyyymm(nyrs-1)/100 resP@gsnPanelMainString = "SST: "+season+": "+yStrt+"-"+yLast ;******************************************* ; first plot ;******************************************* do n=0,neof-1 res@gsnLeftString = "EOF "+(n+1) res@gsnRightString = sprintf("%5.1f", eof@pcvar(n)) +"%" plot(n)=gsn_csm_contour_map(wks,eof(n,:,:),res) end do gsn_panel(wks,plot,(/neof,1/),resP) ; now draw as one plot ;******************************************* ; second plot ;******************************************* ; EOF time series [bar form] rts = True rts@gsnDraw = False ; don't draw yet rts@gsnFrame = False ; don't advance frame yet rts@gsnScale = True ; force text scaling ; these four resources allow the user to stretch the plot size, and ; decide exactly where on the page to draw it. rts@vpHeightF = 0.40 ; Changes the aspect ratio rts@vpWidthF = 0.85 rts@vpXF = 0.10 ; change start locations rts@vpYF = 0.75 ; the plot rts@gsnXYLineChart = True rts@tiYAxisString = "SST" ; y-axis label rts@gsnYRefLine = 0. ; reference line rts@gsnXYBarChart = False ; create bar chart rts@gsnAboveYRefLineColor = "tomato" ; above ref line fill red rts@gsnBelowYRefLineColor = "skyblue" ; below ref line fill blue ; panel plot only resources rtsP = True ; modify the panel plot rtsP@gsnMaximize = True ; large format rtsP@gsnPanelMainString = "SST: "+season+": "+yStrt+"-"+yLast year = yyyymm/100 ; create individual plots do n=0,neof-1 rts@gsnLeftString = "EOF "+(n+1) rts@gsnRightString = sprintf("%5.1f", eof@pcvar(n)) +"%" plot(n) = gsn_csm_xy (wks,year,eof_ts(n,:),rts) end do gsn_panel(wks,plot,(/neof,1/),rtsP) ; now draw as one plot end