NCL做ERA5水汽场的EOF分解

清木77

;********************************************************************************
;**********水汽场 数据ERA5**************1979-2018*********************************
;****************0.25*0.25********************************************************
;***********EOF分析***************************************************************
;********************************************************************************
begin
    yearstart               = 1979
    yearend                 = 2018  

;设置范围

    latmin_p                = 34
    latmax_p                = 56
    lonmin_p                = 45
    lonmax_p                = 88

; 变量
    year                    = ispan(yearstart, yearend, 1)
    month                   = ispan(1, 12, 1)

;路径
    file_path               = "/Users/Desktop/zhongya/data/"
    pic_out                 = "/Users/Desktop/zhongya/eof_ncl/res/eof_ca"
    pic_out2                = "/Users/Desktop/zhongya/eof_ncl/res/eof_xy"

;************data*********************
    ;读取数据
    f                       = addfile(file_path + "ERA5_month_Total_column_water_vapour_1979_2020.nc", "r")   
    ; print(f)
    vapor_all               = short2flt(f->tcwv(0:479, {latmin_p:latmax_p}, {lonmin_p:lonmax_p}) )   ;提取时间范围数据
    vapor_all@units         = "mm"

    utc                     = cd_calendar(vapor_all&time, 0)    ; 转换时间，求出年份

    ; 创建变量存储年总数据
    vap_year                = month_to_annual(vapor_all(:, :, :), 1)    ;年累计0，年平均1
    vap_year!0              = "time"
    vap_year&time           = year      ; 增加年坐标
    ; printVarSummary(vap_year)

    vap_rmv                 = dim_rmvmean_n_Wrap(vap_year, 0) ;距平
;*************计算EOF************
    ;计算权重
    vap_lat                 = doubletofloat(vap_rmv&latitude)   
    wgt                     = sqrt(cos(vap_lat*0.01745329))         ;权重系数
    wp                      = vap_rmv * conform(vap_rmv, wgt, 1)
    copy_VarCoords(vap_rmv, wp)

    ;EOF分解
    x                       = wp(latitude|:,longitude|:,time|:)
    neof                    = 6
    eof                     = eofunc_Wrap(x, neof, False)
    ;north检验
    dims_p                  = dimsizes(vap_rmv)
    ntim                    = dims_p(0)
    sig_pcv                 = eofunc_north(eof@pcvar, ntim, False)  
printVarSummary(eof)
print(sig_pcv)
    ;计算时间序列，标准化
    ;eof_ts                  = eofunc_ts_Wrap(x, eof, False)               ;EOF时间序列
    eof_ts                  = eofunc_ts_n_Wrap(vap_rmv, eof, False, 0)
    eof_t                   = dim_standardize_n_Wrap(eof_ts, 1, 1)
;*******************************************
;  plots about EOF
;*******************************************
    wks                                 = gsn_open_wks("png", pic_out)        
    plot                                = new(neof,graphic)

    res                                 = True         
    res@gsnDraw                         = False         ; don't draw yet
    res@gsnFrame                        = False         ; don't advance frame yet
    res@gsnAddCyclic                    = False         ; data not cyclic

    ; 地图范围
    res@mpMaxLatF                       = latmax_p
    res@mpMinLatF                       = latmin_p
    res@mpMaxLonF                       = lonmax_p
    res@mpMinLonF                       = lonmin_p

    ; 国界设置
    res@mpDataSetName                   = "/Users/Desktop/zhongya/data/map/map_earth/Earth..4"  ; 修正地图数据包
    res@mpDataBaseVersion               = "MediumRes"
    res@mpOutlineBoundarySets           = "National"    ;国界边界
    ; 国界线条设置
    res@mpNationalLineThicknessF        = 2
    res@mpOutlineOn                     = True      
    res@mpLandFillColor                 = "white"
    res@mpInlandWaterFillColor          = "white"
    res@mpFillOn                        = False

    res@cnFillPalette                   = "BlRe"      ; choose colormap
    res@cnFillOn                        = True          ; turn on color fill
    res@cnLinesOn                       = False         ; True is default
    res@lbLabelBarOn                    = False         ;;关闭每个图形公共的色标，以公用一个色标

    res@cnLineLabelsOn                  = False
    res@cnInfoLabelOn                   = False
    ;为用一个公共色标，为每幅图设定相同等值线
    res@cnLevelSelectionMode            = "ExplicitLevels"   ;"ManualLevels"    ;
    res@cnLevels                        = ispan(-16, 16, 2)*0.001
    ;res@cnFillColors                    = (/3, 25, 34, 40, 46, 53, 58, 63, 68, 75, 81, 88, 95/)    ;
    ;res@cnFillColors                    = (/3, 10, 14, 19, 26, 32, 37, 42, 46, 53, 75, 95/)    ;

    resm                                = res

    ;X Y 坐标设置
    Lon_space                           = 8      ;经度间隔
    Lat_space                           = 4      ;纬度间隔

    resm@tmXBMode                       = "Explicit"
    resm@tmXBValues                     = ispan(lonmin_p, lonmax_p, Lon_space)
    resm@tmXBLabels                     = tostring(ispan(lonmin_p, lonmax_p, Lon_space)) + "~S~o~N~E"
    resm@tmXBMinorOn                    = False        ;小刻度显隐开关.默认True
    resm@tmXTOn                         = False   ;X没有顶轴，
    resm@tmXBLabelFont                  = 25
    ;resm@tmXBLabelFontThicknessF        = 0.03
    resm@tmXBLabelFontHeightF           = 0.02       ;设置X轴下标签高度0.02-0.06

    resm@tmYLMode                       = "Explicit"
    resm@tmYLValues                     = ispan(latmin_p, latmax_p, Lat_space)
    resm@tmYLLabels                     = tostring(ispan(latmin_p, latmax_p, Lat_space)) + "~S~o~N~N"
    resm@tmYLMinorOn                    = False
    resm@tmYROn                         = False   ;Y没有右轴
    resm@tmYLLabelFont                  = 25
    ;resm@tmYLLabelFontThicknessF        = 0.03  
    resm@tmYLLabelFontHeightF           = 0.02   ;设置Y轴左标签高度  0.02-0.06

    do n=0, neof-1
        resm@gsnLeftString              = "~F25~EOF "+(n+1)
        resm@gsnLeftStringOrthogonalPosF= 0.02
        resm@gsnRightString             = sprintf("~F25~%5.1f", eof@pcvar(n)) +"%"
        resm@gsnRightStringOrthogonalPosF= 0.02
        plot(n)                         = gsn_csm_contour_map(wks, eof(n, :, :), resm)

    end do

    ; panel plot onlyresources
    resP                                = True         ; modify the panel plot
    resP@gsnMaximize                    = True         ; large format
    resP@gsnPanelLabelBar               = True         ; add common colorbar

    ; resP@gsnPanelMainString             = "EOF"
    ; resP@gsnPanelMainFont               = 25
    resP@lbLabelFontHeightF             = 0.005         ;色标标注字体大小
    resP@lbLabelFont                    = 25
    gsn_panel(wks,plot,(/neof,1/),resP)     ; draw all 'neof' as one plot
;EOF时间序列的res设定
;*******************************************
; time series (principalcomponent) plot
;*******************************************
    wks2                                = gsn_open_wks("png", pic_out2)
    ;eof_t@long_name                     = "Amplitude"

    rts                                 = True
    rts@gsnDraw                         = False       ; don't draw yet
    rts@gsnFrame                        = False       ; don't advanceframe yet           

    rts@vpHeightF                       = 0.40        ; Changes the aspect ratio
    rts@vpWidthF                        = 0.85
    rts@vpXF                            = 0.10        ; change startlocations
    rts@vpYF                            = 0.75        ; the plot

    rts@gsnYRefLine                     = 0.              ; 参考线  
    rts@gsnXYBarChart                   = False          ;True
    ; rts@gsnAboveYRefLineColor           = "red"           ; above refline fill red
    ; rts@gsnBelowYRefLineColor           = "blue"          ; below refline fill blue

    rts@xyMarkLineMode                  = "Marklines"
    rts@xyMarker                        = 16
    rts@xyMarkerColor                   = "black"
    rts@xyMarkerSizeF                   = 0.01

    rts@tmXBLabelFont                   = 25
    rts@tmYLLabelFont                   = 25
    rts@tmXBMinorOn                     = False        ;小刻度显隐开关.默认True
    rts@tmXTOn                          = False   ;X没有顶轴
    rts@tmYLMinorOn                     = False
    rts@tmYROn                          = False

    rts@tiYAxisString                   = ""

    rts@tmXBMode                        = "Explicit"
    rts@tmXBValues                      = (/1,11,21,31/)
    rts@tmXBLabels                      = (/1980,1990,2000,2010/)

    ; rts@trYMinF                         = -3.5
    ; rts@trYMaxF                         = 3.5

    do n=0,neof-1
        rts@gsnLeftString               = "~F25~EOF "+(n+1)
        rts@gsnLeftStringOrthogonalPosF = 0.02
        rts@gsnRightString              = sprintf("~F25~%5.1f", eof@pcvar(n)) +"%"
        rts@gsnRightStringOrthogonalPosF= 0.02
        ;plot(n)                         = gsn_csm_xy(wks2,year,eof_t(n,:),rts)
        plot(n)                         = gsn_csm_y(wks2, eof_t(n,:), rts)
    end do                                    
    ; panel plot onlyresources
    rtsP                                = True             ; modify the panel plot
    rtsP@gsnMaximize                    = True             ; large format
    gsn_panel(wks2,plot,(/neof,1/),rtsP)        ; draw all 'neof' as one plot



end

felin

感谢分享，很好的脚本，收藏啦！