- 积分
- 1216
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2018-9-28
- 最后登录
- 1970-1-1
|
发表于 2020-6-8 10:35:09
|
显示全部楼层
本帖最后由 liweier123 于 2020-6-8 10:39 编辑
- 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/contrib/time_axis_labels.ncl"
- ;----------------------------------------------------------------------
- ; Converts YYYYMM time to "days since" time.
- ;----------------------------------------------------------------------
- function yyyymm_to_time(time)
- local ntime, stime, year, month, day, hour, minute, second, units
- begin
- ntime = dimsizes(time)
- stime = tostring(time)
- year = toint(str_get_cols(stime,0,3))
- month = toint(str_get_cols(stime,4,5))
- day = new(ntime,integer)
- hour = new(ntime,integer)
- minute = new(ntime,integer)
- second = new(ntime,double)
-
- ;----------------------------------------------------------------------
- ; Note: you may need to change the day, hour, minute values depending
- ; on your definition of "YYYYMM". Is this the first day and the
- ; first hour of the month? Is it the middle of the month?
- ;----------------------------------------------------------------------
- day = 1
- hour = 0
- minute = 0
- second = 0.d
- units = "days since " + min(year) + "-01-01 00:00"
- newtime = cd_inv_calendar(year,month,day,hour,minute,second,units,0)
- return(newtime)
- end
- ;----------------------------------------------------------------------
- ; Main code
- ;----------------------------------------------------------------------
- begin
- shpname = "/cygdrive/g/Tibet/add_TP/DBATP_Line.shp"
- ;afiles = systemfunc("ls /cygdrive/g/Tibet/cloudsat/cloud-height/height*.nc")
- files = systemfunc("ls /cygdrive/g/Tibet/cloudsat/data/height_monthly*.nc")
- f = addfiles(files, "r")
- ListSetType(f, "cat")
- ;printVarSummary(f)
- data = f[:]->valavg_monthly
- ;printVarSummary(data)
- dims = dimsizes(data)
- ntim = dims(0)
- time = data&time
- ;time!0 = "time"
- ;time@units = "YYYYMM"
- ;print(time)
- ; convert integer YYYYMM to float
- ;timeF = yyyymm_to_yyyyfrac(time,0)
- ;data&time := timeF
- latS = 25.0
- latN = 40.0
- lonL = 70.0
- lonR = 105.0
- ;======================================== 余弦加权
- lat = data&lat
- clat = sqrt(cos(lat*0.01745329))
- CLAT = conform(data, clat, 1)
- dataw = data*CLAT
- copy_VarCoords(data, dataw)
- ;printVarSummary(dataw)
- ;======================================== 对加权数据进行EOFs分解并显著性检验
- neof = 3
- optEOF = True
- optEOF@jopt = 0 ;默认值,计算使用协方差矩阵
- optETS = False
- eof = eofunc_n_Wrap(dataw, neof, optEOF, 0)
- eof = smth9_Wrap(eof, 0.5, 0.5, True)
- eof_ts = eofunc_ts_n_Wrap(dataw, eof, optETS, 0)
- printVarSummary(eof) ; examine EOF variables
- printVarSummary(eof_ts)
- eof_ts = dim_standardize_n_Wrap(eof_ts, 0, 1);标准化时间序列 0为样本标准差即N-1 1为进行标准化的维度
- pcvar = eof@pcvar ;获取方差贡献属性项
- prinfo = True
- sig_ev = eofunc_north(eof@eval, ntim, prinfo) ;差异显著性
- ;======================================== 定义时间序列eof_ts的属性
- eof_ts!0 = "eofs"
- ;eof_ts!1 = "time"
- eof_ts&eofs = (/1,2,3/) ;三个模态
- ;eof_ts&time = year
- ;print(eof_ts(0,:))
- ;============================================================================
- wks = gsn_open_wks("png", "EOF_MaxdBZ_height_month_pinghua-07-10")
- stringw = (/"(a)","(b)","(c)","(d)","(e)","(f)"/)
- ploteof = new ( 6, "graphic")
- ;=================================================;
- ;空间场相关设置
- ;=================================================;
- res = True ; contour plots
- res@gsnDraw = False
- res@gsnFrame = False
- res@gsnLeftString = ""
- res@gsnRightString = ""
- res@cnFillOn = True ; turn on color
- res@cnFillMode = "RasterFill" ; turn on raster mode
- res@cnLinesOn = False ; turn off contour lines
- res@cnLineLabelsOn = False
- res@cnInfoLabelOn = False
- res@cnFillPalette = "BlWhRe" ; specify colormap
- res@mpShapeMode = "FreeAspect"
- res@vpHeightF = 0.3
- res@vpWidthF = 0.6
- res@gsnAddCyclic = False ; force cyclic value
- res@lbLabelBarOn = False ; turn off individual lb's
- res@mpFillOn = False
- res@mpMinLatF = latS ; range to zoom in on
- res@mpMaxLatF = latN
- res@mpMinLonF = lonL
- res@mpMaxLonF = lonR
- res@tmXBLabelFontHeightF = 0.02
- res@tmYLLabelFontHeightF = 0.02
- res@gsnLeftStringFontHeightF = 0.02
- res@gsnRightStringFontHeightF = 0.02
- symMinMaxPlt(eof, 40, False, res) ; 设置几种颜色
- do n=0,2
- res@gsnLeftString = stringw(n)+"EOF"+(n+1)
- res@gsnRightString = sprintf("%5.1f", pcvar(n)) +"%"
- ploteof(n) = gsn_csm_contour_map(wks,eof(n,:,:),res)
- end do
- tpres =True
- tpres@gsLineColor = "black" ;设置廓线颜色
- tpres@gsLineThicknessF = 3 ;设置廓线宽度
- poly0=gsn_add_shapefile_polylines(wks,ploteof(0),shpname,tpres) ;plot为底图
- poly1=gsn_add_shapefile_polylines(wks,ploteof(1),shpname,tpres)
- poly2=gsn_add_shapefile_polylines(wks,ploteof(2),shpname,tpres)
-
- ;===========================================
- ; 时间场相关设置
- ;===========================================
- ;---Convert the time values and reassign as x's time coordinate array
- ; newtime = yyyymm_to_time(time)
- ; eof_ts&time := newtime
- ; printVarSummary(eof_ts)
- rts = True
- rts@gsnDraw = False ; don't draw yet
- rts@gsnFrame = False ; don't advance frame yet
- ; rts@gsnScale = True
- ;rts@tmXBFormat= "f" ; no unnessary decimal pts
- rts@vpHeightF = 0.40 ; Changes the aspect ratio
- rts@vpWidthF = 0.6
- ;rts@vpXF = 0.20 ; change start locations
- ;rts@vpYF = 0.80 ; the plot
- 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
- rts@trXMinF = min(time)
- rts@trXMaxF = max(time) ; nicer end bound
- rts@tiXAxisString = "Time"
- rts@tiYAxisString = ""
- ; rts@tmXBMode = "Explicit"
- ; rts@tmXBValues = time(::12)
- ; rts@tmXBLabels = (/"Jan ~C~2007","Jan ~C~2008","Jan ~C~2009","Jan ~C~2010"/)
- rts@tmXBLabelFontHeightF = 0.02
- rts@tmYLLabelFontHeightF = 0.02
- res@tmLabelAutoStride = True ; nice stride on labels
- rts@gsnLeftStringFontHeightF = 0.02
- ;---Set resources necessary to nicely format X axis
- ; restick = True
- ; restick@ttmFormat = "%c%Y" ; Jan 2010, Feb 2010, etc.
- ;time_axis_labels(newtime,rts,restick)
- do n=0,2
- rts@gsnLeftString = stringw(n+3)+"PC"+(n+1)
- ploteof(n+3) = gsn_csm_xy(wks,time,eof_ts(n,:),rts)
- end do
- resP = True ; panel plots
- resP@gsnMaximize = True ; large format
- resP@gsnPanelLabelBar = True ; add common colorbar
- resP@gsnPaperOrientation = "portrait" ; force portrait [landscape]
- resP@amJust = "TopLeft"
- gsn_panel(wks,ploteof,(/2,3/),resP) ; draw all 'neof' as one plot
- end
复制代码
图只是eof时间序列,子程序试过也画不出来 |
-
|