Hi all,
I am a newer in using NCL. Recently I used NCL to analyze the cmip5 GFDL precipitation data and I want to get the EOF result of the Indo-Pacific region. But I encountered some errors that confused me.
Below is the scripts:
==============================================================
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 = -25.
latN = 25.
lonL = 40.
lonR = 160.
yrStrt = 1979
yrLast = 2008
season = "DJF" ; choose Dec-Jan-Feb seasonal mean
neof = 1 ; 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
; ==============================================================
fils = systemfunc("ls pr_Amon_GFDL-HIRAM-C180_amip_r1i1p1_*.nc")
f = addfiles (fils,"r")
; printVarSummary(f)
ListSetType (f,"cat") ; concatenate (=default)
TIME = f[:]->time
YYYY = cd_calendar(TIME,-1)/100 ; entire file
iYYYY = ind(YYYY.ge.yrStrt .and. YYYY.le.yrLast)
prate = f[:]->pr(iYYYY,:,:)
printVarSummary(prate) ; variable overview
; ==============================================================
; compute desired global seasonal mean: month_to_season (contributed.ncl)
; ==============================================================
prate=rmMonAnnCycTLL(prate)
PRA= month_to_season (prate, season)
nyrs = dimsizes(PRA&time)
PRA= PRA*24*3600
; =================================================================
; create weights: sqrt(cos(lat)) [or sqrt(gw) ]
; =================================================================
rad = 4.*atan(1.)/180.
clat = f[1]->lat
clat = sqrt( cos(rad*clat) ) ; gw for gaussian grid
printVarSummary(clat)
; =================================================================
; weight all observations
; =================================================================
wPRA = PRA ; copy meta data
wPRA = (/dble2flt(PRA*conform(PRA, clat, 1))/)
printVarSummary(wPRA)
wPRA@long_name = "Wgt: "+wPRA@long_name
; =================================================================
; Reorder (lat,lon,time) the *weighted* input data
; Access the area of interest via coordinate subscripting
; =================================================================
x = wPRA({lat|latS:latN},{lon|lonL:lonR},time|:)
eof = eofunc_Wrap(x, neof, optEOF)
eof_ts = eofunc_ts_Wrap (x, eof, optETS)
printVarSummary( eof ) ; examine EOF variables
printVarSummary( eof_ts )
; =================================================================
; Normalize time series: Sum spatial weights over the area of used
; =================================================================
dimx = dimsizes( x )
mln = dimx(1)
sumWgt = mln*sum( clat({lat|latS:latN}) )
eof_ts = eof_ts/sumWgt
; =================================================================
; Extract the YYYYMM from the time coordinate
; associated with eof_ts [same as PRA&time]
; =================================================================
yyyymm = cd_calendar(eof_ts&time,-2)/100
;;yrfrac = yyyymm_to_yyyyfrac(yyyymm, 0.0); not used here
;============================================================
; PLOTS
;============================================================
wks = gsn_open_wks("ps","eof_test1")
gsn_define_colormap(wks,"BlWhRe") ; choose colormap
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
;---This resource not needed in V6.1.0
res@gsnSpreadColors = True ; spread out color table
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@lbLabelBarOn = False ; turn off individual lb's
; 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
resP@lbLabelAutoStride = True ; auto stride on labels
yStrt = yyyymm(0)/100
yLast = yyyymm(nyrs-1)/100
resP@tiMainString = "PRA: "+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_ce(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 rtsources 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@tiYAxisString = "mm/day" ; y-axis label
rts@gsnYRefLine = 0. ; reference line
rts@gsnXYBarChart = True ; create bar chart
rts@gsnAboveYRefLineColor = "red" ; above ref line fill red
rts@gsnBelowYRefLineColor = "blue" ; below ref line fill blue
; panel plot only resources
rtsP = True ; modify the panel plot
rtsP@gsnMaximize = True ; large format
rtsP@tiMainString = "PRA: "+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
************************************************************
and the result is :
Variable: prate
Type: float
Total Size: 298598400 bytes
74649600 values
Number of Dimensions: 3
Dimensions and sizes: [time | 360] x [lat | 360] x [lon | 576]
Coordinates:
time: [1111.5..12038.5]
lat: [-89.75..89.75]
lon: [0.3125..359.6875]
Number Of Attributes: 11
long_name : Precipitation
units : kg m-2 s-1
cell_methods : time: mean
interp_method : conserve_order1
missing_value : 1e+20
_FillValue : 1e+20
standard_name : precipitation_flux
original_units : kg/m2/s
original_name : precip
cell_measures : area: areacella
associated_files : baseURL: http://cmip-pcmdi.llnl.gov/CMIP5/dataLocation areacella: areacella_fx_GFDL-HIRAM-C180_amip_r0i0p0.nc
Variable: clat
Type: double
Total Size: 2880 bytes
360 values
Number of Dimensions: 1
Dimensions and sizes: [lat | 360]
Coordinates:
lat: [-89.75..89.75]
Number Of Attributes: 5
long_name : latitude
units : degrees_north
bounds : lat_bnds
standard_name : latitude
axis : Y
warning:Dimension (0) has not been defined
warning:Dimension (1) has not been defined
warning:Dimension (2) has not been defined
Variable: wPRA
Type: float
Total Size: 24883200 bytes
6220800 values
Number of Dimensions: 3
Dimensions and sizes: [time | 30] x [lat | 360] x [lon | 576]
Coordinates:
time: [1111.5..11703.5]
lat: [-89.75..89.75]
lon: [0.3125..359.6875]
Number Of Attributes: 13
NMO : 0
_FillValue : 1e+20
anomaly_op_ncl : Annual Cycle Removed: rmMonAnnCycTLL: contributed.ncl
associated_files : baseURL: http://cmip-pcmdi.llnl.gov/CMIP5/dataLocation areacella: areacella_fx_GFDL-HIRAM-C180_amip_r0i0p0.nc
cell_measures : area: areacella
original_name : precip
original_units : kg/m2/s
standard_name : precipitation_flux
missing_value : 1e+20
interp_method : conserve_order1
cell_methods : time: mean
units : kg m-2 s-1
long_name : DJF: Precipitation
Variable: eof
Type: float
Total Size: 76800 bytes
19200 values
Number of Dimensions: 3
Dimensions and sizes: [evn | 1] x [lat | 100] x [lon | 192]
Coordinates:
evn: [1..1]
lat: [-24.75..24.75]
lon: [40.3125..159.6875]
Number Of Attributes: 7
eval_transpose : 19.68391
eval : 13030.75
pcvar : 13.04149
matrix : covariance
method : transpose
_FillValue : 1e+20
long_name : EOF: Wgt: DJF: Precipitation
Variable: eof_ts
Type: float
Total Size: 120 bytes
30 values
Number of Dimensions: 2
Dimensions and sizes: [evn | 1] x [time | 30]
Coordinates:
evn: [1..1]
time: [1111.5..11703.5]
Number Of Attributes: 4
ts_mean : -0.7275095
matrix : covariance
_FillValue : 1e+20
long_name : EOF: Amplitude: Wgt: DJF: Precipitation
fatal:["NclVar.c":1376]:Assignment type mismatch, right hand side can't be coerced to type of left hand side
fatal:Execute: Error occurred at or near line 95 in file eof_pra_gfdl.ncl
****************************************
I would appreciate that if you could help me solve this problem. Thank you very much.
Best regards,
Hai
----------------
Hai Wang
College of Physical and Environmental Oceanography,
Ocean University of China
238, Songling Road, Qingdao, Shandong, China, 266100
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Mon Jul 9 08:43:46 2012
This archive was generated by hypermail 2.1.8 : Mon Jul 09 2012 - 10:45:32 MDT