- 积分
- 23225
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2013-10-31
- 最后登录
- 1970-1-1
|
发表于 2021-7-20 10:30:27
|
显示全部楼层
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/esmf/ESMF_regridding.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
begin
a = addfile("$NCARG_ROOT/lib/ncarg/data/cdf/landsea.nc","r")
lsdata = a->LSMASK
data_dir1 = "/disk_raid5/pub/lzb-data/others/lanzhoushuoshi/data/"
f2 = addfile(data_dir1+"hgt-2.nc","r")
pre = f2->z(:,4,:,:)
;pre = (pre_pri*pre_pri@scale_factor+pre_pri@add_offset) * 1000
dim_pre = dimsizes(pre)
pre_4dim = onedtond(ndtooned(pre), (/dim_pre(0)/12,12,dim_pre(1),dim_pre(2)/))
pre_sum_ave = dim_avg_n(pre_4dim(:,5:7,:,:), 1)
copy_VarCoords(pre(0,:,:), pre_sum_ave(0,:,:))
pre_sum_lonf = lonFlip(pre_sum_ave)
lsm = landsea_mask(lsdata,pre_sum_lonf&lat,pre_sum_lonf&lon)
pre_sum_lonf_mask = mask(pre_sum_lonf,conform(pre_sum_lonf,lsm,(/1,2/)).eq.0,False)
copy_VarCoords(pre_sum_lonf(0,:,:),pre_sum_lonf_mask(0,:,:))
;************************select-region*******************************************
pre_field = pre_sum_lonf_mask(:,{0:90},:)
pre_field!0 = "time"
dim_pre_field = dimsizes(pre_field)
pre_field_re = new((/dim_pre_field(0),(dim_pre_field(1)*dim_pre_field(2))/),typeof(pre_field),pre_field@_FillValue)
lat_p = new((/dim_pre_field(1)*dim_pre_field(2)/),"integer")
lon_p = new((/dim_pre_field(1)*dim_pre_field(2)/),"integer")
np = 0
do i = 0,dim_pre_field(1)-1
do j = 0,dim_pre_field(2)-1
ind_miss = ind(ismissing(pre_field(:,i,j)))
if(ismissing(ind_miss(0))) then
pre_field_re(:,np) = pre_field(:,i,j)
lat_p(np) = i
lon_p(np) = j
np = np + 1
end if
delete(ind_miss)
end do
end do
ind_pre_miss = ind(ismissing(pre_field_re(0,:)))
ind_pre_miss_first = ind_pre_miss(0)-2
pre_field_re!0 = "time"
pre_field_re!1 = "grids"
pre_field_re_reorder = pre_field_re(grids|0:ind_pre_miss_first,time|:)
;*************************************************************************************
data_dir2 = "/disk_raid5/pub/lzb-data/others/lanzhoushuoshi/data/"
f5 = addfile(data_dir2+"mrso1-2.nc","r")
f6 = addfile(data_dir2+"mrso2-2.nc","r")
f7 = addfile(data_dir2+"mrso3-2.nc","r")
swvl1_pri = f5->swvl1(24:,:,:)
swvl2_pri = f6->swvl2(24:,:,:)
swvl3_pri = f7->swvl3(24:,:,:)
swvl1 = swvl1_pri*swvl1_pri@scale_factor+swvl1_pri@add_offset
swvl2 = swvl2_pri*swvl2_pri@scale_factor+swvl2_pri@add_offset
swvl3 = swvl3_pri*swvl3_pri@scale_factor+swvl3_pri@add_offset
dim_swvl = dimsizes(swvl1_pri)
swvl_all = new((/3,dim_swvl(0),dim_swvl(1),dim_swvl(2)/), typeof(swvl1), swvl1@_FillValue)
swvl_all(0,:,:,:) = swvl1
swvl_all(1,:,:,:) = swvl2
swvl_all(2,:,:,:) = swvl3
swvl_all_ave = dim_avg_n(swvl_all, 0)
swvl_4dim = onedtond(ndtooned(swvl_all_ave), (/dim_swvl(0)/12,12,dim_swvl(1),dim_swvl(2)/))
swvl_spr_ave = dim_avg_n(swvl_4dim(:,2:4,:,:), 1)
copy_VarCoords(swvl1_pri(0,:,:), swvl_spr_ave(0,:,:))
swvl_lonf = lonFlip(swvl_spr_ave)
;************************************************************swvl1
ksm = landsea_mask(lsdata,swvl_lonf&lat,swvl_lonf&lon)
swvl_lonf_mask = mask(swvl_lonf,conform(swvl_lonf,lsm,(/1,2/)).eq.0,False)
copy_VarCoords(swvl_lonf(0,:,:),swvl_lonf_mask(0,:,:))
swvl_field = swvl_lonf_mask(:,{0:90},:)
dim_swvl_field = dimsizes(swvl_field)
swvl_field_re = new((/dim_swvl_field(0),(dim_swvl_field(1)*dim_swvl_field(2))/),typeof(swvl_field),swvl_field@_FillValue)
lat_s = new((/dim_swvl_field(1)*dim_swvl_field(2)/),"integer")
lon_s = new((/dim_swvl_field(1)*dim_swvl_field(2)/),"integer")
ns = 0
do i = 0,dim_swvl_field(1)-1
do j = 0,dim_swvl_field(2)-1
ind_miss = ind(ismissing(swvl_field(:,i,j)))
if( ismissing(ind_miss(0))) then
swvl_field_re(:,ns) = swvl_field(:,i,j)
lat_s(ns) = i
lon_s(ns) = j
ns = ns + 1
end if
delete(ind_miss)
end do
end do
ind_swvl_miss = ind(ismissing(swvl_field_re(0,:)))
ind_swvl_miss_first = ind_swvl_miss(0)-2
swvl_field_re!0 = "time"
swvl_field_re!1 = "grids"
swvl_field_re_reorder = swvl_field_re(grids|0:ind_swvl_miss_first,time|:)
dim_lft = dimsizes(pre_field_re_reorder)
dim_rgt = dimsizes(swvl_field_re_reorder)
;SVD--------------------------------------------------------------------------------------------
ntime = dim_lft(1) ; # time steps
ncolr = dim_lft(0) ; # columns (stations or grid pts) for S
ncols = dim_rgt(0) ; # columns (stations or grid pts) for Z
nsvd = 3 ; # svd patterns to calculate
xmsg = -1e+30 ; missing value
homlft = new((/nsvd,ncolr/),float)
hetlft = new((/nsvd,ncolr/),float)
homrgt = new((/nsvd,ncols/),float)
hetrgt = new((/nsvd,ncols/),float)
x = svdstd(pre_field_re_reorder,swvl_field_re_reorder,nsvd,homlft,hetlft,homrgt,hetrgt)
print("svdstd: percent variance= " + x)
ak = onedtond(x@ak,(/nsvd,ntime/))
bk = onedtond(x@bk,(/nsvd,ntime/))
ak!0 = "sv"
ak!1 = "time"
ak&time = ispan(1981,2010,1)
bk!0 = "sv"
bk!1 = "time"
bk&time = ispan(1981,2010,1)
dim_ak = dimsizes(ak)
dim_bk = dimsizes(bk)
ak_time = ispan(1, dim_ak(1), 1)
bk_time = ispan(1, dim_bk(1), 1)
;------------------------------------------------------------------------------------------------
homlft_re = new((/nsvd,dim_pre_field(1),dim_pre_field(2)/),float,homlft@_FillValue)
do p = 0,nsvd-1
do i = 0,ncolr-1
homlft_re(p,lat_p(i),lon_p(i)) = homlft(p,i)
end do
end do
copy_VarCoords(pre_field(0,:,:),homlft_re(0,:,:))
hetlft_re = new((/nsvd,dim_pre_field(1),dim_pre_field(2)/),float,hetlft@_FillValue)
do p = 0,nsvd-1
do i = 0,ncolr-1
hetlft_re(p,lat_p(i),lon_p(i))= hetlft(p,i)
end do
end do
copy_VarCoords(pre_field(0,:,:),hetlft_re(0,:,:))
homrgt_re = new((/nsvd,dim_swvl_field(1),dim_swvl_field(2)/),float,homrgt@_FillValue)
do p=0,nsvd-1
do i=0,ncols-1
homrgt_re(p,lat_s(i),lon_s(i))= homrgt(p,i)
end do
end do
copy_VarCoords(swvl_field(0,:,:),homrgt_re(0,:,:))
hetrgt_re = new((/nsvd,dim_swvl_field(1),dim_swvl_field(2)/),float,hetrgt@_FillValue)
do p=0,nsvd-1
do i=0,ncols-1
hetrgt_re(p,lat_s(i),lon_s(i))= hetrgt(p,i)
end do
end do
copy_VarCoords(swvl_field(0,:,:),hetrgt_re(0,:,:))
;------------------------------------------------------------------------------------------------
wks = gsn_open_wks("pdf","fig13")
gsn_define_colormap(wks,"BlWhRe")
plot1 = new(3,graphic)
plot2 = new(3,graphic)
plotaxy = new(3,graphic)
plotbxy = new(3,graphic)
;plott1 = new(6,graphic)
;plott2 = new(6,graphic)
res = True
res@gsnMaximize = True
res@gsnAddCyclic = False
res@gsnMaximize = True ; maximize plot in frame
res@gsnFrame = False
res@gsnDraw = False
res@mpDataSetName = "Earth..4" ; This new database contains
res@mpDataBaseVersion = "LowRes";"MediumRes" ; Medium resolution database
res@mpOutlineOn = True ; Turn on map outlines
res@cnFillDrawOrder = "PreDraw" ; Make sure map fill happens
res@mpAreaMaskingOn = True
res@mpOceanFillColor = 1
res@mpInlandWaterFillColor = 0
res@mpFillOn = False
res@mpOutlineOn = True
res@mpMaxLatF = 80
res@mpMaxLonF = 180
res@mpMinLatF = 0
res@mpMinLonF = -180
res@tmXTOn = False ; Turn off top and right
res@tmYROn = False ; tickmarks.
res@mpGeophysicalLineThicknessF = 1
res@mpNationalLineThicknessF = 1
res@tmBorderThicknessF = 2
res@tmXBLabelFontHeightF = 0.016
res@tmYLLabelFontHeightF = 0.016
res@gsnLeftStringFontHeightF = 0.016
res@gsnRightStringFontHeightF = 0.016
res@tmYLMode = "Explicit"
res@tmYLValues = (/0,20,40,60,80/) ;
res@tmYLLabels = (/"EQ","20~S~o~N~N","40~S~o~N~N","60~S~o~N~N","80~S~o~N~N"/)
res@tmXBMode = "Explicit"
res@tmXBValues = (/-180,-120,-60,0,60,120,180/)
res@tmXBLabels = (/"180~S~o~N","120~S~o~N~W","60~S~o~N~W","0~S~o~N","60~S~o~N~E","120~S~o~N~E","180~S~o~N"/)
res@tmXBMinorOn = False
res@tmYLMinorOn = False
res@mpAreaMaskingOn = True
res@mpMaskAreaSpecifiers = (/"Land"/)
res@mpFillOn = True
res@mpLandFillColor = 0
res@mpInlandWaterFillColor = 0
res@mpOceanFillColor = 0
res@lbLabelBarOn = False ; no individual label bar
cnres = res
cnres@cnFillOn = True
cnres@cnLinesOn = True
cnres@cnLevelSelectionMode = "ExplicitLevels"
;cnres@cnFillColors = (/27,32,37,42,47,57,62,67,72,77/)
;cnres@cnLevels = (/-0.8,-0.6,-0.4,-0.2,0,0.2,0.4,0.6,0.8/)
cnres@pmLabelBarOrthogonalPosF = 0.25
cnres@gsnContourNegLineDashPattern = 2
wrf_smooth_2d(hetlft_re,3)
wrf_smooth_2d(hetrgt_re,3)
cnres@gsnLeftString = "(a)SVD1_MAM SM"
cnres@gsnRightString = sprintf("%5.1f", x(0)) +"%"
plot1(0) = gsn_csm_contour_map(wks,hetlft_re(0,:,:),cnres)
cnres@gsnLeftString = "(a)SVD1_JJA 500hPa hgt"
cnres@gsnRightString = sprintf("%5.1f", x(0)) +"%"
plot2(0) = gsn_csm_contour_map(wks,hetrgt_re(0,:,:),cnres)
cnres@gsnLeftString = "(b)SVD2_MAM SM"
cnres@gsnRightString = sprintf("%5.1f", x(1)) +"%"
plot1(1) = gsn_csm_contour_map(wks,hetlft_re(1,:,:),cnres)
cnres@gsnLeftString = "(b)SVD2_JJA 500hPa hgt"
cnres@gsnRightString = sprintf("%5.1f", x(1)) +"%"
plot2(1) = gsn_csm_contour_map(wks,hetrgt_re(1,:,:),cnres)
cnres@gsnLeftString = "(c)SVD3_MAM SM"
cnres@gsnRightString = sprintf("%5.1f", x(2)) +"%"
plot1(2) = gsn_csm_contour_map(wks,hetlft_re(2,:,:),cnres)
cnres@gsnLeftString = "(c)SVD3_JJA 500hPa hgt"
cnres@gsnRightString = sprintf("%5.1f", x(2)) +"%"
plot2(2) = gsn_csm_contour_map(wks,hetrgt_re(2,:,:),cnres)
res2 = True
res2@gsnMaximize = True ; maximize plot in frame
res2@gsnFrame = False ;False
res2@gsnDraw = False ;False
res2@vpWidthF = 0.8
res2@vpHeightF = 0.2
res2@vpXF = 0.15
res2@xyLineThicknesses = 3
res2@trYMaxF = 40
res2@trYMinF = -40
res2@trXMaxF = 1981
res2@trXMinF = 2010
res2@tmXTOn = False
res2@tmYROn = False
res2@tmXBLabelFontHeightF = 0.016
res2@tmYLLabelFontHeightF = 0.016
res2@gsnLeftStringFontHeightF = 0.016
res2@gsnRightStringFontHeightF = 0.016
res2@pmLegendDisplayMode = "Always" ; turn on legend
res2@pmLegendWidthF = 0.05 ; Change width and
res2@pmLegendHeightF = 0.03 ; height of legend.
res2@lgPerimOn = False
res2@lgLabelFontHeightF = .011
res2@gsnLeftString = "(a)SVD1"
res2@gsnRightString = sprintf("%5.1f", x(0)) +"%"
res2@xyDashPatterns = 1
res2@xyExplicitLegendLabels = (/"500hPa hgt(JJA)"/)
res2@pmLegendOrthogonalPosF = -1.18 ; more neg = down
res2@pmLegendParallelPosF = 0.88 ; move units right
plotaxy(0) = gsn_csm_xy(wks,ispan(1981,2010,1),ak(0,:),res2)
delete(res2@gsnRightString)
res2@pmLegendOrthogonalPosF = -1.25 ; more neg = down
res2@pmLegendParallelPosF = 0.88 ; move units right
res2@xyExplicitLegendLabels = (/"Soil Moisture(MAM)"/) ; labels for the legend
res2@xyDashPatterns = 0
plotbxy(0) = gsn_csm_xy(wks,ispan(1981,2010,1),bk(0,:),res2)
overlay(plotaxy(0),plotbxy(0))
res2@gsnLeftString = "(b)SVD2"
res2@gsnRightString = sprintf("%5.1f", x(1)) +"%"
res2@xyDashPatterns = 1
res2@xyExplicitLegendLabels = (/"500hPa hgt(JJA)"/)
res2@pmLegendOrthogonalPosF = -1.18 ; more neg = down
res2@pmLegendParallelPosF = 0.88 ; move units right
plotaxy(1) = gsn_csm_xy(wks,ispan(1981,2010,1),ak(1,:),res2)
delete(res2@gsnRightString)
res2@pmLegendOrthogonalPosF = -1.25 ; more neg = down
res2@pmLegendParallelPosF = 0.88 ; move units right
res2@xyExplicitLegendLabels = (/"Soil Moisture(MAM)"/) ; labels for the legend
res2@xyDashPatterns = 0
plotbxy(1) = gsn_csm_xy(wks,ispan(1981,2010,1),bk(1,:),res2)
overlay(plotaxy(1),plotbxy(1))
res2@gsnLeftString = "(c)SVD3"
res2@gsnRightString = sprintf("%5.1f", x(2)) +"%"
res2@xyDashPatterns = 1
res2@xyExplicitLegendLabels = (/"500hPa hgt(JJA)"/)
res2@pmLegendOrthogonalPosF = -1.18 ; more neg = down
res2@pmLegendParallelPosF = 0.88 ; move units right
plotaxy(2) = gsn_csm_xy(wks,ispan(1981,2010,1),ak(2,:),res2)
delete(res2@gsnRightString)
res2@pmLegendOrthogonalPosF = -1.25 ; more neg = down
res2@pmLegendParallelPosF = 0.88 ; move units right
res2@xyExplicitLegendLabels = (/"Soil Moisture(MAM)"/) ; labels for the legend
res2@xyDashPatterns = 0
plotbxy(2) = gsn_csm_xy(wks,ispan(1981,2010,1),bk(2,:),res2)
overlay(plotaxy(2),plotbxy(2))
res1_lines = True ; polyline mods desired
res1_lines@gsLineDashPattern = 2.
res1_lines@gsLineThicknessF = 2. ; line thicker
res1_lines@gsLineColor = "black" ; line color
zz = (/-9999,9999/)
uu = (/0,0/)
dumA1 = gsn_add_polyline(wks,plotaxy(0),zz,uu,res1_lines)
dumA2 = gsn_add_polyline(wks,plotaxy(1),zz,uu,res1_lines)
dumA3 = gsn_add_polyline(wks,plotaxy(2),zz,uu,res1_lines)
res1Panel = True
res1Panel@gsnMaximize = True
res1Panel@gsnFrame = False
res1Panel@gsnPanelLabelBar = True
res1Panel@lbLabelFontHeightF = 0.01
res1Panel@gsnPanelBottom = 0.1
res1Panel@gsnPanelXWhiteSpacePercent = 0
res1Panel@gsnPanelYWhiteSpacePercent = 9
gsn_panel(wks,plot2,(/3,1/),res1Panel)
frame(wks)
gsn_panel(wks,plot1,(/3,1/),res1Panel)
frame(wks)
res2Panel = True
res2Panel@gsnMaximize = True
res2Panel@gsnFrame = False
res2Panel@gsnPanelXWhiteSpacePercent = 0
res2Panel@gsnPanelYWhiteSpacePercent = 10
gsn_panel(wks,plotaxy,(/3,1/),res2Panel)
frame(wks)
end
|
|
窥视卡
雷达卡
提升卡
置顶卡
沉默卡
喧嚣卡
变色卡
抢沙发
千斤顶
显身卡
楼主