ncl画垂直剖面，显示剖面线的问题

chongzika

1. load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
   
2. load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
   
3. 
   
4. begin
   
5. ;
   
6. ; The WRF ARW input file.
   
7. ; This needs to have a ".nc" appended, so just do it.
   
8.   files     = systemfunc("ls -1 wrfout_d02*")
   
9.   a = addfiles(files,"r")
   
10. xlon = wrf_user_getvar(a, "XLONG",0)
    
11. xlat = wrf_user_getvar(a, "XLAT",0)
    
12.     ; gsn_define_colormap(wks,"MPL_gist_yarg")
    
13. 
    
14. ; Set some basic resources
    
15. 
    
16. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    
17. ll_res= True
    
18. loc0= wrf_user_ll_to_ij(a[0],70,37,ll_res)
    
19. print(loc0)
    
20. loc1= wrf_user_ll_to_ij(a[0],80,37,ll_res)
    
21. print(loc1)
    
22. 
    
23. ;---------------------------------------------------------------
    
24. 
    
25.   plane = new(4,float)
    
26.   plane = (/loc0(0),loc1(0),loc0(1), loc1(1)/)  ;(/58,65/)       ;
    
27.   print(plane)
    
28.   opts = True
    
29.   angle = 0.
    
30. 
    
31.   lat_plane = wrf_user_intrp2d(xlat,plane,angle,opts)
    
32.   lon_plane = wrf_user_intrp2d(xlon,plane,angle,opts)
    
33. ;gsn_polyline(wks,plot,(/lon_plane(0),lon_plane(-1)/),(/lat_plane(0),lat_plane(-1)/),lnres)
    
34. end
    

复制代码

输出的起始值经纬度
70,37—— 80,37
读取后插值画图用的
70.21637，36.09188 —— 81.05466，37.04618结果画出来的就只是红色的线，而不是想要的蓝色线条官网脚本 http://www2.mmm.ucar.edu/wrf/OnL ... Section_smooth4.htm

chongzika

1. <div class="blockcode"><blockquote>
   

复制代码

1. load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
   
2. load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
   
3. begin
   
4. files = systemfunc("ls -1 wrfout_d02*")
   
5. a = addfiles(files,"r")
   
6. xlon = wrf_user_getvar(a, "XLONG",0)
   
7. xlat = wrf_user_getvar(a, "XLAT",0)
   
8. 
   
9. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   
10. ll_res= True
    
11. lon0=79.5
    
12. lat0=37
    
13. lon1=80
    
14. lat1=37
    
15. loc0= wrf_user_ll_to_ij(a[0],lon0,lat0,ll_res)
    
16. ;print(loc0)
    
17. loc1= wrf_user_ll_to_ij(a[0],lon1,lat1,ll_res)
    
18. 
    
19. ;---------------------------------------------------------------
    
20. 
    
21. plane = new(4,float)
    
22. plane = (/loc0(0),loc1(0),loc0(1), loc1(1)/) ;(/58,65/) ;
    
23. opts = True
    
24. angle=0.
    
25. lat_plane = wrf_user_intrp2d(xlat,plane,angle,opts)
    
26. lon_plane = wrf_user_intrp2d(xlon,plane,angle,opts)
    
27. dimsX = dimsizes(lat_plane)
    
28. loc2=wrf_user_ij_to_ll(a[0],lon_plane(0),lat_plane(0),ll_res)
    
29. loc3=wrf_user_ij_to_ll(a[0],lon_plane(dimsX-1),lat_plane(dimsX-1),ll_res)
    
30. loc4= wrf_user_ll_to_ij(a[0],loc2(0),loc2(1),ll_res)
    
31. loc5= wrf_user_ll_to_ij(a[0],loc3(0),loc3(1),ll_res)
    
32. print("lon0 lat0 lon1 lat1")
    
33. print (sprintf("%6.2f", lon0) +" " \
    
34. +sprintf("%6.2f", lat0) +" " \
    
35. +sprintf("%6.2f", lon1) +" " \
    
36. +sprintf("%6.2f", lat1) )
    
37. print("loc0_x loc0_y loc1_x loc1_y")
    
38. print (sprintf("%6.2f", loc0(0)) +" " \
    
39. +sprintf("%6.2f", loc0(1)) +" " \
    
40. +sprintf("%6.2f", loc1(0)) +" " \
    
41. +sprintf("%6.2f", loc1(1)) )
    
42. print("loc0_x1 loc0_y1 loc1_x1 loc1_y1")
    
43. print (sprintf("%6.2f", lon_plane(0)) +" " \
    
44. +sprintf("%6.2f", lat_plane(0)) +" " \
    
45. +sprintf("%6.2f", lon_plane(1)) +" " \
    
46. +sprintf("%6.2f", lat_plane(1)) )
    
47. print("lon0_x lat0_y lon1_x lat1_y")
    
48. print (sprintf("%6.2f", loc2(0)) +" " \
    
49. +sprintf("%6.2f", loc2(1)) +" " \
    
50. +sprintf("%6.2f", loc3(0)) +" " \
    
51. +sprintf("%6.2f", loc3(1)) )
    
52. print("loc0_x2 loc0_y2 loc1_x2 loc1_y2")
    
53. print (sprintf("%6.2f", loc4(0)) +" " \
    
54. +sprintf("%6.2f", loc4(1)) +" " \
    
55. +sprintf("%6.2f", loc5(0)) +" " \
    
56. +sprintf("%6.2f", loc5(1)) )
    
57. end

复制代码

输出对比一看，经过坐标转化再插值再转换后完全对不上了
  lon0     lat0      lon1   lat1
   79.50  37.00   80.00   37.00
  loc0_x     loc0_y      loc1_x   loc1_y
  132.00 133.00  136.00  132.00
  loc0_x1     loc0_y1      loc1_x1   loc1_y1
   79.52  37.30   79.56   37.21
  lon0_x     lat0_y      lon1_x   lat1_y
   75.27  27.80   75.29   27.78
  loc0_x2     loc0_y2      loc1_x2   loc1_y2
   80.00  37.00   80.00   37.00

freekiller

你一定用的是ncl6.2后的版本

chongzika

freekiller 发表于 2018-3-19 23:38
你一定用的是ncl6.2后的版本

是的6.2，难道之前的都没问题么？