[ncl-talk] Average Spatial Correlations

[Melissa Lazenby]

Dear NCL Users

I am wanting to create 1 plot of averaged spatial correlations.
I currently have 40 global spatial correlation plots but would like to average those all out to 1 global plot with the overall average spatial correlations at each gridpoint.

I am not sure whether to use the average function or to sum them up first and divide by 40.

Your help would be very much appreciated.

Below is the code I am trying to create this average plot from.

Kindest Regards
Melissa

; ==============================================================
; eof_1.ncl
; ==============================================================

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"

begin
; ==============================================================
; User defined parameters that specify region of globe and
; ==============================================================
  latS   =  -40.
  latN   =  0.
  lonL   =  10.
  lonR   =  60.

  yrStrt = 1975
  yrLast = 2004

  season = "DJF"    ; choose Dec-Jan-Feb seasonal mean

  neof   = 1        ; number of EOFs
  optEOF = True
  optEOF at jopt = 0   ; This is the default; most commonly used; no need to specify.
  optETS = False

model = (/"ACCESS1-0","ACCESS1-3","bcc-csm1-1","bcc-csm1-1-m","BNU-ESM","CanESM2","CCSM4","CESM1-BGC","CESM1-CAM5","CESM1-FASTCHEM","CESM1-WACCM","CMCC-CESM","CMCC-CM","CMCC-CMS","CNRM-CM5","CSIRO-Mk3-6-0","EC-EARTH","FGOALS-g2","GFDL-CM3","GFDL-ESM2G","GFDL-ESM2M","GISS-E2-H","GISS-E2-H-CC","GISS-E2-R","GISS-E2-R-CC","HadGEM2-AO","HadGEM2-CC","HadGEM2-ES","inmcm4","IPSL-CM5A-LR","IPSL-CM5A-MR","IPSL-CM5B-LR","MIROC5","MIROC-ESM","MIROC-ESM-CHEM","MPI-ESM-LR","MPI-ESM-MR","MPI-ESM-P","MRI-CGCM3","NorESM1-M","NorESM1-ME"/)


wks  = gsn_open_wks("X11","ALL_40_MMM_Eof_cor_SST")              ; specifies a plot
gsn_define_colormap(wks,"ncl_default")  ; choose color map
plot = new (dimsizes(model),"graphic")
do gg = 0,dimsizes(model)-1

; read in model data

  in=addfile("/mnt/geog/ml382/melphd/eof_sicz/eof90models/pr_"+model(gg)+"_eof_DJF_Anom.nc","r")

  lat    = in->lat
  lon    = in->lon
  time   = in->time
  YYYY   = cd_calendar(time,-1)/100                 ; entire file
  iYYYY  = ind(YYYY.ge.yrStrt .and. YYYY.le.yrLast)

  if((gg.eq.19) .or. (gg.eq.20)) then   ;(gg.eq.17) .or.
        PR = in->pr(iYYYY,:,:)
    else
  PR    = in->pr(iYYYY,0,:,:)
  end if
;
  nyrs   = dimsizes(PR&time)
; =================================================================
; normalize data at each gridpoint by local standard deviation at each grid pt
; =================================================================
  PR = dim_standardize_n(PR,1,0)

; =================================================================
; Reorder (lat,lon,time) the *weighted* input data
; Access the area of interest via coordinate subscripting
; =================================================================
  x      = PR({lat|latS:latN},{lon|lonL:lonR},time|:)

  eof    = eofunc_Wrap(x, neof, optEOF)
  eof_ts = eofunc_ts_Wrap (x, eof, optEOF)

  printVarSummary( eof )                         ; examine EOF variables
  printVarSummary( eof_ts )

; =================================================================
; Extract the YYYYMM from the time coordinate
; associated with eof_ts [same as SLP&time]
; =================================================================

  yyyymm = cd_calendar(eof_ts&time,-2)/100

; ==============================================================
; Open the file: Read only the user specified period
; ==============================================================
  s      = addfile ("/mnt/nfs2/geog/ml382/melphd/global/sstmodel90/tos_Omon_"+model(gg)+"_historical_safrica_climDJF.nc", "r")


  TIME1   = s->time
  YYYYZ   = cd_calendar(TIME1,-1)/100                 ; entire file
  iYYYYZ  = ind(YYYYZ.ge.yrStrt .and. YYYYZ.le.yrLast)

  sst    = s->tos(iYYYYZ,:,:)
  printVarSummary(sst)

  sst1 =sst(:,:,:)


  printVarSummary(eof_ts)                              ; variable overview
  eof1 = eof_ts(0,:)

 sst1 = lonFlip(sst1)
; =================================================================
; Reorder (latitude,longitude,time) the *weighted* input data
; Access the area of interest via coordinate subscripting
; =================================================================

  q      = eof_ts(0,:)

  y      = sst1(lat|:,lon|:,time|:)

     y&lat at units = "degrees_north"
     y&lon at units = "degrees_east"

     x&lat at units = "degrees_north"
     x&lon at units = "degrees_east"

 printVarSummary(q)
 printVarSummary(y)


 ccr = escorc(q, y)

 printVarSummary(ccr)

  ccr!0 = "lat" ; name dimensions
  ccr!1 = "lon"
  ccr&lat = y&lat ; assign coordinate values and
  ccr&lon = y&lon ; units attributes

 printVarSummary(ccr)


 ;ensccrdimavg = dim_avg_Wrap( ccr(lat|:, lon|:)
 ensccr = new((/dimsizes(lat),dimsizes(lon)/), typeof(ccr))
    ;ensccr!0="lat"
    ;ensccr&lat=lat
    ;ensccr!1="lon"
    ;ensccr&lon=lon
;ensccr(gg) = ensccrsum

;ensccrsum(gg) = sum( ccr(lat|:, lon|: ))
 ;ensccr = ensccrsum / 40

ensccr = dim_avg_n_Wrap(ccr,0)


  ensccr!0 = "lat" ; name dimensions
  ensccr!1 = "lon"
  ensccr&lat = lat ; assign coordinate values and
  ensccr&lon = lon ; units attributes

 printVarSummary(ensccr)
;============================================================
; PLOT
;============================================================


;*******************************************
; first plot
;*******************************************

  rescn                       = True
  rescn at cnFillOn              = True
  rescn at cnLinesOn            = False        ; True is default
  rescn at gsnDraw              = False              ; don't draw
  rescn at gsnFrame             = False              ; don't advance frame
  rescn at cnLineLabelsOn       = False        ; True is default
  rescn at lbLabelBarOn         = False        ; turn off individual lb's

  rescn at cnLevelSelectionMode = "ManualLevels"     ; set manual contour levels
  rescn at cnMinLevelValF       = -0.5                ; set min contour level
  rescn at cnMaxLevelValF       =  0.5                ; set max contour level
  rescn at cnLevelSpacingF      =  0.1               ; set contour spacing

  ;rescn at lbOrientation        = "Vertical"         ; vertical label bar
  rescn at tiMainString = ""+model(gg)+""
;---This resource defaults to True in NCL V6.1.0
  rescn at lbLabelAutoStride    = True               ; optimal label stride

  rescn at gsnSpreadColors      = True               ; use full range of colors
  rescn at mpCenterLonF         = 180.               ; center plot at 180
  rescn at lbLabelFontHeightF = 0.015

  rescn at gsnAddCyclic         = True
 ; reverse the first two colors
  setvalues wks
    "wkColorMap"        : "ncl_default"
    "wkForegroundColor" : (/0.,0.,0./)
    "wkBackgroundColor" : (/1.,1.,1./)
  end setvalues


end do

plot = gsn_csm_contour_map(wks,ensccr,rescn)


draw(plot)
frame(wks)

end


-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20160613/18a949e3/attachment.html