- 积分
- 10054
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2014-3-8
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看，没有帐号？立即注册
x
本帖最后由 海soso 于 2017-4-6 18:34 编辑
大家好 用NCL画台风路径和高度场风场叠加的图后发现台风路径不显示了 NCL小白 不知道是为啥 希望看出问题的盆友们能指导下 炒鸡蟹蟹 贴ncl程序（基本上就是把官网上overlay和一个画台风路径的例子合起来了）：
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
;;read u,v,hgt from the data at 500hPa
f=addfile("/nuist/scratch/luzhangcloud/zl/xiong/hgt_6hr/hgt.2000.nc","r")
hgt = short2flt( f->hgt)
hgt1=hgt(0,5,:,:)
hgt2=hgt1/10
copy_VarCoords(hgt1,hgt2)
f1=addfile("/nuist/scratch/luzhangcloud/zl/xiong/uwind_6hr/uwnd.2000.nc","r")
uwnd = short2flt( f1->uwnd)
uwnd1=uwnd(0,5,:,:)
f2=addfile("/nuist/scratch/luzhangcloud/zl/xiong/vwind_6hr/vwnd.2000.nc","r")
vwnd = short2flt( f2->vwnd)
vwnd1=vwnd(0,5,:,:)
;;read pth from the data
fiTY = "/nuist/scratch/luzhangcloud/zl/xiong/pth/20005/20005.txt"
; 获取文本文件的行数，相应的还有numAsciiCol函数用于获取列数
nrow = 29
;YYYYMMDDHH = new(nrow, "string")
MM = new(nrow, "string")
DD = new(nrow, "string")
HH = new(nrow, "string")
lat = new(nrow, "float")
lon = new(nrow, "float")
vmax = new(nrow, "integer")
mslp = new(nrow, "integer")
; print(nrow)
data = asciiread(fiTY, -1, "string")
;YYYYMMDDHH = str_get_field(data, 3, " ")
MM = str_get_field(data, 3, " ")
DD = str_get_field(data, 4, " ")
HH = str_get_field(data, 5, " ")
lat = stringtofloat(str_get_field(data, 6, " "))
lon = stringtofloat(str_get_field(data, 7, " "))
vmax = stringtoint(str_get_field(data, 8, " "))
mslp = stringtoint(str_get_field(data, 9, " "))
;DateChar = stringtochar(YYYYMMDDHH)
;MM = chartostring(DateChar(:,5:6))
;DD = chartostring(DateChar(:,7:8))
;HH = chartostring(DateChar(:,9:10))
;HHi = mod(stringtoint(YYYYMMDDHH), 100)
;DDi = mod(stringtoint(YYYYMMDDHH)/100, 100)
;MMi = mod(stringtoint(YYYYMMDDHH)/10000, 100)
;;create plots;;
wks = gsn_open_wks("png","wnd_hgt_pth") ; send graphics to PNG file
; gsn_define_colormap(wks,"BkBlAqGrYeOrReViWh200")
res = True
res@gsnDraw = False
res@gsnFrame = False
res@gsnMaximize = True
res@tmXTOn = False
res@tmYROn = False
res@gsnLeftString = ""
res@gsnRightString = ""
res@tiMainString = "20005track-000702-000709"
;;set map;;
mpres = res
mpres@mpDataSetName = "Earth..4"
mpres@mpDataBaseVersion = "MediumRes"
mpres@mpOutlineOn = True
mpres@mpOutlineSpecifiers = (/"China:states","Taiwan"/)
mpres@mpOutlineBoundarySets = "National"
mpres@mpGeophysicalLineThicknessF = 2
mpres@mpNationalLineThicknessF = 2
mpres@mpFillDrawOrder = "PostDraw"
mpres@mpFillOn = True
mpres@mpFillAreaSpecifiers = (/"water", "land" /)
mpres@mpSpecifiedFillColors = (/"deepskyblue2","white"/)
; mpres@mpSpecifiedFillColors = (/100,0/)
mpres@mpMaskAreaSpecifiers = (/"China:states","Taiwan"/)
; 绘制洋面经纬网格
mpres@mpGridAndLimbOn = "True"
mpres@mpGridMaskMode = "MaskNotOcean"
mpres@mpGridLineDashPattern = 15
mpres@mpGridSpacingF = 5.0
;;set area;;
mpres@mpMinLatF = 10
mpres@mpMaxLatF = 50
mpres@mpMinLonF = 110
mpres@mpMaxLonF = 160
;;set contour;;
cnres = res
cnres@cnFillDrawOrder = "PreDraw"
cnres@cnFillOn = True
cnres@cnLinesOn = True
cnres@cnLineThicknesses = True
cnres@cnLineThicknessF = 3
cnres@cnLineLabelsOn = True ; turn off line labels
cnres@pmLabelBarWidthF = 0.4
cnres@pmLabelBarHeightF = 0.05
cnres@pmLabelBarOrthogonalPosF = 0.1
cnres@lbLabelFontHeightF = 0.006
cnres@lbLabelAngleF = 45
; Older way to subset a color map
; cnres@cnFillPalette = "BkBlAqGrYeOrReViWh200"
; cnres@gsnSpreadColorStart = 50
; cnres@gsnSpreadColorEnd = 120
; Newer way to subset a color map
cmap = read_colormap_file("BkBlAqGrYeOrReViWh200")
cnres@cnFillPalette = cmap(25:120,:)
cnres@gsnLeftString = "Temp"
;;set vector;;
res_vc = res
res_vc@vcGlyphStyle = "LineArrow"
res_vc@vcLineArrowThicknessF = 5
res_vc@vcMinDistanceF = 0.01
res_vc@vcRefLengthF = 0.03
;;wind barb resources don't apply
;; res_vc@vcGlyphStyle = "WindBarb"
;; res_vc@vcWindBarbLineThicknessF = 5
;; res_vc@vcWindBarbColor = "Gray40"
res_vc@vcRefAnnoOn = True
res_vc@vcRefMagnitudeF = 30
res_vc@vcRefAnnoString1 = "30"
res_vc@vcRefAnnoSide = "Top"
res_vc@vcRefAnnoString2On = False
res_vc@vcRefAnnoPerimOn = False
res_vc@vcRefAnnoOrthogonalPosF = -0.12
res_vc@vcRefAnnoParallelPosF = 0.999
res_vc@vcRefAnnoBackgroundColor = "Purple"
res_vc@vcVectorDrawOrder = "PostDraw"
res_vc@gsnRightString = "Wind"
;;plot;;
;map = gsn_csm_map(wks,mpres)
pth = gsn_csm_map_ce(wks,mpres)
contour = gsn_csm_contour(wks,hgt2,cnres)
vector = gsn_csm_vector(wks,uwnd1,vwnd1,res_vc)
;;pth
;================
; 按照vmax(单位：节 knot，海里/小时)变量提供的风速大小对台风强度进行区分，
; 并在绘制时用不同颜色标出
; 34~63 热带风暴
; 64~82 台风
; 83~95 强台风
; 96~112 超强台风
;113~136
colours = (/3, 20, 60, 120, 180/)
pthresDot = True
pthresLine = True
dumDot= new(nrow, graphic)
dumLine = new(nrow, graphic)
; 绘制线
pthresLine@gsLineThicknessF = 3
do i = 0, nrow-2
xx = (/ lon(i), lon(i+1)/)
yy = (/ lat(i), lat(i+1)/)
if (vmax(i) .ge. 34 .and. vmax(i).le.63) then
pthresLine@gsLineColor = colours(0)
end if
if (vmax(i) .ge. 64 .and. vmax(i).le.82) then
pthresLine@gsLineColor = colours(1)
end if
if (vmax(i).ge.83 .and. vmax(i) .le. 95) then
pthresLine@gsLineColor = colours(2)
end if
if (vmax(i).ge.96 .and. vmax(i) .lt. 112) then
pthresLine@gsLineColor = colours(3)
end if
if (vmax(i).ge.113 .and. vmax(i) .lt. 136) then
pthresLine@gsLineColor = colours(4)
end if
dumLine(i) = gsn_add_polyline(wks, pth, xx, yy, pthresLine)
end do
; 绘制00时的点
pthresDot@gsMarkerIndex = 1
pthresDot@gsMarkerSizeF = 0.02
pthresDot@gsMarkerColor = "black"
do i = 0, nrow-1
if (HH(i) .eq. "00") then
dumDot(i) = gsn_add_polymarker(wks, pth, lon(i), lat(i), pthresDot)
end if
end do
; 绘制图例
pthresLg = True
pthresLg@lgItemType = "MarkLines"
pthresLg@lgMonoMarkerIndex = True
pthresLg@lgMarkerColors = colours
pthresLg@lgMarkerIndex = 1
pthresLg@lgMarkerSizeF = 0.02
pthresLg@lgMonoDashIndex = True
pthresLg@lgDashIndex = 0
pthresLg@lgLineColors = colours
pthresLg@lgLineThicknessF = 3
pthresLg@vpWidthF = 0.14
pthresLg@vpHeightF = 0.13
pthresLg@lgPerimFill = 0
pthresLg@lgPerimFillColor = "Background"
pthresLg@lgLabelFontHeightF = 0.08
pthresLg@lgTitleFontHeightF = 0.015
pthresLg@lgTitleString = "Best Track"
lbid = gsn_create_legend(wks, 5, (/" 34kt to 63kt", \
" 64 to 82kt", " 83 to 95kt", " 96 to 112kt", \
" 113 to 136kt"/), pthresLg)
; 将图例放置在图中
ampthres = True
ampthres@amParallelPosF = -0.3
ampthres@amOrthogonalPosF = -0.3
dumLg = gsn_add_annotation(pth, lbid, ampthres)
; 标注00时的日期
dumDate = new(nrow,graphic)
pthresTx = True
pthresTx@txFontHeightF = 0.012
pthresTx@txFontColor = "black"
pthresTx@txJust = "BottomLeft"
do i = 1, nrow-1
if (HH(i) .ne. "00" ) then
continue
end if
dumDate(i) = gsn_add_text(wks,pth, MM(i)+DD(i)+HH(i), lon(i)+0.5, lat(i), pthresTx)
end do
;;overlay filled contours and vectors on the map;;
overlay(pth,contour)
overlay(pth,vector)
;;add text;;
txres = True
txres@txFontHeightF = 0.02
txres@txFontColor = "Purple"
txres@txBackgroundFillColor = "White"
txres@txFontOpacityF = 0.8
txres@txFontThicknessF = 4.0
dum = gsn_add_text(wks,(/vector/),"500hPa",105,52.7,txres)
;;drawing "map" will draw everything: map, contours, vectors, and text;;
draw(pth)
frame(wks)
end
附我的错的图：
叠加之前的的路径图：（很明显，叠加以后这个路径就没有了T.T）
求救啊！！
|
|
窥视卡
雷达卡
发表于 2017-4-6 18:33:39
提升卡
置顶卡
沉默卡
喧嚣卡
变色卡
抢沙发
千斤顶
显身卡
{:5_213:}
楼主