关于格点数据中指定两点间的剖面图

youxia1873

1. load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"  
   
2. load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"  
   
3. 
   
4. begin
   
5. f = addfile("个人路径/20210725.nc","r")
   
6. pv=short2flt(f->pv)           ;era5产品读取需要short2flt函数
   
7. w=short2flt(f->w)
   
8. pv=pv(:,:,::-1,:)
   
9. w=w(:,:,::-1,:)
   
10. 
    
11. wks = gsn_open_wks("png","profile")
    
12. gsn_define_colormap(wks,"MPL_RdBu")
    
13. time = f->time
    
14. time2 = time+8
    
15. copy_VarMeta(time, time2)
    
16. chn_data = cd_calendar(time2,0)              ;era5产品的时间以秒计算，所以要换算
    
17. year   = tointeger(chn_data(:,0))
    
18. month  = tointeger(chn_data(:,1))
    
19. day    = tointeger(chn_data(:,2))
    
20. hour   = tointeger(chn_data(:,3))
    
21. minute = tointeger(chn_data(:,4))
    
22. ;chn_date = sprinti("%0.4i-", year)+sprinti("%0.2i-", month)+sprinti("%0.2i_", day)+sprinti("%0.2i", hour)
    
23. chn_date = sprinti("%0.2i", day)+sprinti("%0.2i", hour)
    
24. chn = tointeger(chn_date)
    
25. ;print(chn_date)
    
26. lon = f->longitude
    
27. lat = f->latitude
    
28. lat = lat(::-1)
    
29. lev = f->level
    
30. downlat = 29.996        ;设置起始点经纬度
    
31. downlon = 122.223
    
32. toplat  = 39.098
    
33. toplon  = 119.210
    
34. npts = 40
    
35. 
    
36. dist = gc_latlon(downlat,downlon,toplat,toplon,npts,2)   ;插值起止点间的40个点位的经纬度
    
37. lat_new = dist@gclat
    
38. copy_VarAtts(lat,lat_new)
    
39. lon_new = dist@gclon
    
40. w_new  = linint2_points(lon, lat, w, False, lon_new, lat_new, 0)    ;将要绘制的变量进行插值
    
41. pv_new = linint2_points(lon, lat, pv, False, lon_new, lat_new, 0)
    
42. pv_new = smth9(pv_new,0.5,-0.25,False)     ;平滑
    
43. w_profil =w_new(1,:,:)
    
44. w_profil!0 ="level"
    
45. w_profil!1 = "grid"
    
46. w_profil&level = lev
    
47. w_profil&grid  = lat_new
    
48. pv_profil =pv_new(1,:,:)*1000000
    
49. pv_profil!0 ="level"
    
50. pv_profil!1 = "grid"
    
51. pv_profil&level = lev
    
52. pv_profil&grid  = lat_new
    
53. printVarSummary(pv)
    
54. ;print(pv)
    
55. 
    
56. res=True
    
57. res@gsnFrame=False
    
58. res@gsnDraw=False
    
59. res@gsnLeftString  = " "
    
60. res@gsnRightString = " "
    
61. res@tmXBMode       = "Automatic"
    
62. res@tmXBMaxTicks   = 8
    
63. res@pmTickMarkDisplayMode   = "Always"
    
64. res@gsnAddCyclic   = False
    
65. res@tmXTOn         = False
    
66. resw=res
    
67. resw@cnLinesOn = False
    
68. resw@cnFillOn  = True
    
69. resw@cnLevelSelectionMode = "ManualLevels"
    
70. resw@cnMaxLevelValF = 3
    
71. resw@cnMinLevelValF = -3
    
72. resw@gsnAddCyclic   = False
    
73. resw@cnInfoLabelOn  = False
    
74. plota = gsn_csm_pres_hgt(wks,w_profil(10:,:),resw)
    
75. resp =res
    
76. resp@cnFillOn  = False
    
77. resp@cnLinesOn = True
    
78. resp@cnInfoLabelOn  = False
    
79. resp@cnLineLabelsOn = False
    
80. resp@cnMonoLineThickness  = False
    
81. resp@cnLineThicknesses    = (/1,3/)
    
82. resp@cnLevelSelectionMode = "ExplicitLevels"
    
83. resp@cnLevels = (/2, 0.5/)
    
84. plotb = gsn_csm_pres_hgt(wks,pv_profil(10:,:),resp)
    
85. overlay(plota,plotb)             ;叠加
    
86. draw(wks)
    
87. frame(wks)
    
88. end

复制代码