- 积分
- 15764
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2015-3-8
- 最后登录
- 1970-1-1
|
发表于 2016-3-29 19:58:04
|
显示全部楼层
版主，可以麻烦您帮我看下我的求u-veof脚本哪错了吗？
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
begin
latS = 0.0
latN =30.0
lonL = 110.0
lonR = 160.0
yrStrt = 1871
yrLast = 2000
season = "MAM" ; choose chun seasonal
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
f1=addfile("uwnd.mon.mean.nc","r")
f2=addfile("vwnd.mon.mean.nc","r")
TIME = f1->time
YYYY = cd_calendar(TIME,-1)/100 ; entire file
iYYYY = ind(YYYY.ge.yrStrt .and. YYYY.le.yrLast)
uwnd= short2flt(f1->uwnd)
vwnd= short2flt(f2->vwnd)
u=uwnd(iYYYY,3,:,:)
v=vwnd(iYYYY,3,:,:)
printVarSummary(u)
printVarSummary(v)
lat=f1->lat
lon=f1->lon
u1=month_to_season(u,"MAM" )
v1=month_to_season(v,"MAM" )
;nyrs = dimsizes(lh&time)
printVarSummary(u1)
printVarSummary(v1)
u2=dtrend_leftdim(u1,False)
v2=dtrend_leftdim(v1,False)
copy_VarMeta(u1,u2)
copy_VarMeta(v1,v2)
printVarSummary(u2)
printVarSummary(v2)
; u3=u2({lev|850},lat|:, lon|:, time|:)
CF1 = ezfftf_n(u2,0) ; ezfftf works on right dim
printVarSummary(CF1) ; [2] x [12] x [121] x [120]
cf1 = CF1
cf1(:,:,:,0:12)= 0 ; set waves 2 and higher to 0.0
u4 = ezfftb_n (cf1, 0,1) ; [12] x [121] x [240]
copy_VarMeta(u2, u4)
; uu=u4(time|:,{lev|850},lat|:, lon|:) ; [time | 12] x [lat | 121] x [lon | 240]
printVarSummary(u4)
; v3=v2({lev|850},lat|:, lon|:, time|:)
CF2 = ezfftf_n(v2,0) ; ezfftf works on right dim
printVarSummary(CF2) ; [2] x [12] x [121] x [120]
cf2 = CF2
cf2(:,:,:,0:12)= 0.0 ; set waves 2 and higher to 0.0
v4 = ezfftb_n (cf2, 0.0,1) ; [12] x [121] x [240]
copy_VarMeta(v2, v4)
; vv=v4(time|:,{lev|850},lat|:, lon|:) ; [time | 12] x [lat | 121] x [lon | 240]
printVarSummary(v4)
;YYYY = cd_calendar(time,-1)
;sv= uv2vr_cfd (u4(:,:,:),v4(:,:,:),lat,lon, 3)
; sv=sv*-1000000
;copy_VarMeta(u4,sv)
s=dim_standardize_n_Wrap(u4,1,0)
t=dim_standardize_n_Wrap(v4,1,0)
printVarSummary(s)
printVarSummary(t)
;s=dim_rmvmean_n_Wrap(snew,0)
ss=s(0:28,:,:)
;y=t(0:28,:,:)
; copy_VarMeta(s, ss)
;printVarSummary(ss)
ss(0,:,:)=s(time|7,lat|:,lon|:)
ss(1,:,:)=s(time|10,lat|:,lon|:)
ss(2,:,:)=s(time|15,lat|:,lon|:)
ss(3,:,:)=s(time|18,lat|:,lon|:)
ss(4,:,:)=s(time|26,lat|:,lon|:)
ss(5,:,:)=s(time|29,lat|:,lon|:)
ss(6,:,:)=s(time|32,lat|:,lon|:)
ss(7,:,:)=s(time|35,lat|:,lon|:)
ss(8,:,:)=s(time|41,lat|:,lon|:)
ss(9,:,:)=s(time|44,lat|:,lon|:)
ss(10,:,:)=s(time|48,lat|:,lon|:)
ss(11,:,:)=s(time|55,lat|:,lon|:)
ss(12,:,:)=s(time|60,lat|:,lon|:)
ss(13,:,:)=s(time|71,lat|:,lon|:)
ss(14,:,:)=s(time|81,lat|:,lon|:)
ss(15,:,:)=s(time|83,lat|:,lon|:)
ss(16,:,:)=s(time|87,lat|:,lon|:)
ss(17,:,:)=s(time|93,lat|:,lon|:)
ss(18,:,:)=s(time|95,lat|:,lon|:)
ss(19,:,:)=s(time|99,lat|:,lon|:)
ss(20,:,:)=s(time|102,lat|:,lon|:)
ss(21,:,:)=s(time|106,lat|:,lon|:)
ss(22,:,:)=s(time|107,lat|:,lon|:)
ss(23,:,:)=s(time|109,lat|:,lon|:)
ss(24,:,:)=s(time|112,lat|:,lon|:)
ss(25,:,:)=s(time|117,lat|:,lon|:)
ss(26,:,:)=s(time|121,lat|:,lon|:)
ss(27,:,:)=s(time|124,lat|:,lon|:)
ss(28,:,:)=s(time|127,lat|:,lon|:)
printVarSummary(ss)
y=t(0:28,:,:) ; note the longitude coord
y(0,:,:)=t(time|7,{lat|:},{lon|:})
y(1,:,:)=t(time|10,lat|:,lon|:)
y(2,:,:)=t(time|15,lat|:,lon|:)
y(3,:,:)=t(time|18,lat|:,lon|:)
y(4,:,:)=t(time|26,lat|:,lon|:)
y(5,:,:)=t(time|29,lat|:,lon|:)
y(6,:,:)=t(time|32,lat|:,lon|:)
y(7,:,:)=t(time|35,lat|:,lon|:)
y(8,:,:)=t(time|41,lat|:,lon|:)
y(9,:,:)=t(time|44,lat|:,lon|:)
y(10,:,:)=t(time|48,lat|:,lon|:)
y(11,:,:)=t(time|55,lat|:,lon|:)
y(12,:,:)=t(time|60,lat|:,lon|:)
y(13,:,:)=t(time|71,lat|:,lon|:)
y(14,:,:)=t(time|81,lat|:,lon|:)
y(15,:,:)=t(time|83,lat|:,lon|:)
y(16,:,:)=t(time|87,lat|:,lon|:)
y(17,:,:)=t(time|93,lat|:,lon|:)
y(18,:,:)=t(time|95,lat|:,lon|:)
y(19,:,:)=t(time|99,lat|:,lon|:)
y(20,:,:)=t(time|102,lat|:,lon|:)
y(21,:,:)=t(time|106,lat|:,lon|:)
y(22,:,:)=t(time|107,lat|:,lon|:)
y(23,:,:)=t(time|109,lat|:,lon|:)
y(24,:,:)=t(time|112,lat|:,lon|:)
y(25,:,:)=t(time|117,lat|:,lon|:)
y(26,:,:)=t(time|121,lat|:,lon|:)
y(27,:,:)=t(time|124,lat|:,lon|:)
y(28,:,:)=t(time|127,lat|:,lon|:)
snew=ss(0:28,:,0:179)
snew(0:28,:,0:89)=ss(0:28,:,0:89)
snew(0:28,:,90:179)=y(0:28,:,0:89)
printVarSummary(snew)
; ==============================================================
; compute desired global seasonal mean: month_to_season (contributed.ncl)
; ==============================================================
;OLR = month_to_season (olr, season)
nyrs = dimsizes(snew&time)
; printVarSummary(OLR)
; =================================================================
; create weights: sqrt(cos(lat)) [or sqrt(gw) ]
; =================================================================
rad = 4.*atan(1.)/180.
clat = snew(time|0,lat|:,lon|0)
;clat = f1->lat(50)
clat = sqrt( cos(rad*clat) ) ; gw for gaussian grid
; =================================================================
; weight all observations
; =================================================================
wss = snew ; copy meta data
wss = snew*conform(snew, clat, 1)
wss@long_name = "Wgt: "+wss@long_name
; =================================================================
; Reorder (lat,lon,time) the *weighted* input data
; Access the area of interest via coordinate subscripting
; =================================================================
x = wss({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 = -1.0*eof_ts/sumWgt
fnc = "eofmam.uv.mnmean.nc"
system("/bin/rm -f "+fnc) ; rm any pre-existing file
fout = addfile(fnc, "c") ; new netCDF file
fout@title = "EOF of Vora: MAM"
fout->EOF = eof
fout->EOF_TS = 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
;;yrfrac = yyyymm_to_yyyyfrac(yyyymm, 0.0); not used here
;============================================================
; PLOTS
;============================================================
wks = gsn_open_wks("ps","eofuvanommamsta1")
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@mpCenterLonF = 180.0
res@cnInfoLabelOn = False
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@txString = "ReVorAnom: "+season
;*******************************************
; 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 = "" ; 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@txString = "U-V Anom: "+season
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
出来的结果是这样的
Variable: snew
Type: float
Total Size: 1900080 bytes
475020 values
Number of Dimensions: 3
Dimensions and sizes: [time | 29] x [lat | 91] x [lon | 180]
Coordinates:
time: [685896..1737792]
lat: [90..-90]
lon: [ 0..178]
Number Of Attributes: 20
NMO : 3
statistic_method : Ensemble mean is calculated by averaging over all 56 ensemble members at each time step and then averaging mean over all time steps in a month
actual_range : ( -81.62, 86.64999 )
valid_range : ( -280, 350 )
standard_name : northward_wind
parent_stat : Individual Obs
statistic : Ensemble Mean
level_desc : Pressure Levels
dataset : NOAA-CIRES 20th Century Reanalysis version 2 Monthly Averages
var_desc : v-wind
GRIB_name : VGRD
GRIB_id : 132
precision : 2
units : m/s
long_name : MAM: Monthly V-wind on Pressure Levels
cell_methods : time: mean (monthly from 6-hourly values)
level : 850
missing_value : -9.96921e+36
standardize_op_ncl : dim_standardize_n over dimension(s): time
_FillValue : -9.96921e+36
fatal:NclOneDValGetRangeIndex: Non-monotonic coordinate value being used, can't complete coordinate subscript
fatal:Could not obtain coordinate indexes, unable to perform subscript
fatal:["Execute.c":8575]:Execute: Error occurred at or near line 171 in file eofuvanommamsta1.ncl
|
|
窥视卡
雷达卡
显身卡