NCL绘制CALIPSO L2 VFM图像

灭火器

读取数据主要参考了CALIPSO官网、NCL的dim_gbits函数页面，和hdfeos网站的例子：
[1]https://www-calipso.larc.nasa.gov/resources/calipso_users_guide/data_summaries/vfm/
[2]http://www.ncl.ucar.edu/Document/Functions/Built-in/dim_gbits.shtml
[3]hdfeos.org/zoo/index_openLaRC_Examples.php#CALIPSO

需要注意的点有：
（1）读取Latitude数据时，数组大小是(4224,1)。但如果直接用f->Latitude读取，NCL会用hdf文件中大小为63360的fakeDim2作为Latitude的第一维，导致原本的Latitude数组后面多出了许多0值。因此手动用f->Latitude(0:4223, 0)读取。
（2）变量Feature_Classification_Flags的类型是ushort，长度为16 bit。对于一个二进制的数，函数dim_gbits是从左往右读的，而CALIPSO官网所说的1-3 bits表示类型，是从右往左读的。因此使用函数前需要换算一下。
（3）作图时还用上了aerosol subtype中的dust。
（4）图中把-0.5~8.2km，8.2~20.2km的廓线都画了出来。因为直接用数组来画contour时，NCL依据数组的大小而不是coordinate的相对大小来画图，所以把8.2~20.2km内的一个bin分为两个相同的bin来画图。也许用res@sfXArray和res@sfYArray来设置，就不用这么麻烦。
代码在NCL V6.40下运行，具体如下：

1. ;----------------------------------------------------------------
   
2. ;# 2019/04/17
   
3. ;# Draw CALIOP feature type and sub type between -0.5km to 20.2 km.
   
4. ;# Reference:
   
5. ;  https://www-calipso.larc.nasa.gov/resources/calipso_users_guide/data_summaries/vfm/
   
6. ;----------------------------------------------------------------
   
7. begin
   
8.     filename = "../CALIOP/CAL_LID_L2_VFM-Standard-V4-20.2015-04-18T04-57-10ZD.hdf"
   
9.     file_vfm = addfile(filename, "r")
   
10.    
    
11.     ; Dim of Latitude is 0:4223, total 4224.
    
12.     latmin = 30.0
    
13.     latmax = 40.0
    
14.     lat_raw = file_vfm->Latitude(0:4223,0)
    
15.     ii = ind(lat_raw .ge. latmin .and. lat_raw .le. latmax)
    
16.     lat = lat_raw(ii)
    
17.    
    
18.     ; -0.5km to 8.2km, profile 0, subscripts are 1165:5514.
    
19.     fcf1 = file_vfm->Feature_Classification_Flags(ii,1165:1454)
    
20.     dims1 = dimsizes(fcf1)
    
21. 
    
22.     ; 8.2km to 20.2km, profile 0, subscripts are 165:364.
    
23.     fcf2 = file_vfm->Feature_Classification_Flags(ii, 165:364)
    
24.     dims2 = dimsizes(fcf2)
    
25. 
    
26.     ; Select 1-3 bits, Feature Type
    
27.     ; 0 = invalid
    
28.     ; 1 = clear air
    
29.     ; 2 = cloud
    
30.     ; 3 = aerosol
    
31.     ; 4 = stratospheric feature
    
32.     ; 5 = surface
    
33.     ; 6 = subsurface
    
34.     ; 7 = no signal
    
35.     fType1 = toint(dim_gbits(fcf1, 13, 3, 13, dims1(1)))
    
36.     subType1 = toint(dim_gbits(fcf1, 4, 3, 13, dims1(1)))
    
37.     fType2 = toint(dim_gbits(fcf2, 13, 3, 13, dims2(1)))
    
38.     subType2 = toint(dim_gbits(fcf2, 4, 3, 13, dims2(1)))
    
39.     ; When Feature type = 3 and subtype = 2, then type is dust, let it be 8.
    
40.     fType1 = where((fType1 .eq. 3 .and. subType1 .eq. 2), 8, fType1)
    
41.     fType2 = where((fType2 .eq. 3 .and. subType2 .eq. 2), 8, fType2)
    
42.     ; Reverse the height dim, set index 0 = surface.
    
43.     fType1 = fType1(:,::-1)
    
44.     fType2 = fType2(:,::-1)
    
45. 
    
46.     ; Use array data to store fType1 and fType2.
    
47.     ; Note that the vertical resolution of fType2 is twice as long as fType1.
    
48.     data = new((/dims1(0), dims1(1)+2*dims2(1)/), integer)
    
49.     data(:,0:dims1(1)-1) = fType1
    
50.     ; Because NCL draws plot according to the array size rather than coordinate
    
51.     ; size, so divide one bin in fType2 into 2 bins.
    
52.     do i = 0, dims2(1)-1
    
53.         data(:,dims1(1)+2*i) = fType2(:,i)
    
54.         data(:,dims1(1)+2*i+1) = fType2(:,i)
    
55.     end do
    
56. 
    
57.     ; Create altitude.
    
58.     ; Height is from -0.5km to 8.2km.
    
59.     dz1 = 30.0/1000.0    ; Units is km.
    
60.     hgt1 = (ispan(1, 290, 1)-0.5)*dz1 - 0.5
    
61. 
    
62.     ; Height2 is from 8.2km to 20.2km.
    
63.     dz2 = 60.0/1000.0    ; Units is km.
    
64.     hgt2 = (ispan(1, 400, 1)-0.5)*dz1 + 8.2
    
65. 
    
66.     hgt = new(690, float)
    
67.     hgt(0:289) = hgt1
    
68.     hgt(290:689) = hgt2
    
69.     hgt@long_name = "Height (km)"
    
70.     hgt@units = "km"
    
71.    
    
72.     ; Create latitude.
    
73.     lat@long_name = "Latitude"
    
74.     lat@units = "degrees_north"
    
75.    
    
76.     data!0 = "lat"
    
77.     data&lat = lat
    
78.     data!1 = "hgt"
    
79.     data&hgt = hgt
    
80.    
    
81.     ; Draw latitude VS height plot.
    
82.     picname = "./pic/VFM_v_dust"
    
83.     wks = gsn_open_wks("png", picname)
    
84.     res = True
    
85.     res@gsnDraw = False
    
86.     res@gsnFrame = False
    
87.     colormap = (/"white", "aliceblue", "lightskyblue", "gold", "red", \
    
88.                  "seagreen", "palegreen", "black", "orange4"/)
    
89.     levels = (/0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5/)
    
90.     lbstring = (/"invalid", "clear air", "cloud", "aerosol", "stratospheric", \
    
91.                  "surface", "subsurface", "no signal", "dust"/)
    
92. 
    
93.     res@gsnMaximize = True
    
94.     res@cnFillOn = True
    
95.     res@cnLinesOn = False
    
96.     res@cnLineLabelsOn = False
    
97.     res@cnFillMode = "RasterFill"
    
98.     res@cnMissingValFillColor = 0
    
99.     res@cnMissingValFillPattern = 0
    
100.     res@cnFillPalette = colormap
     
101.     res@cnLevelSelectionMode = "ExplicitLevels"
     
102.     res@cnLevels = levels
     
103. 
     
104.     res@tmYROn = False
     
105.     res@tmXTOn = False
     
106.     ; res@tmYLMode = "Explicit"
     
107.     ; res@tmYLValues = ispan(0, 20, 5)*1.0
     
108.     ; res@tmYLLabels = res@tmYLValues
     
109.    
     
110.     res@lbLabelAutoStride = True
     
111.     res@lbBoxLinesOn = False
     
112.     res@lbLabelAlignment = "BoxCenters"
     
113.     res@lbLabelStrings = lbstring
     
114.     res@lbBoxMinorExtentF = 0.2
     
115.     res@lbLabelFontHeightF = 0.009
     
116.    
     
117.     res@tiMainString = ""
     
118.     res@vpWidthF = 0.8
     
119.     res@vpHeightF = 0.2
     
120.    
     
121.     ; Transpose fType so x=lat, y=hgt.
     
122.     plot = gsn_csm_contour(wks, transpose(data), res)  
     
123. 
     
124.     draw(plot)
     
125.     frame(wks)
     
126. end
     

复制代码

Wetter

写的很赞

meng_xlei

请问为什么第19行是1454不是5514呢？

meng_xlei

Wetter 发表于 2019-7-20 00:39
写的很赞

请问为什么第19行是1454不是5514呢？

灭火器

meng_xlei 发表于 2019-7-26 21:11
请问为什么第19行是1454不是5514呢？

如图，1165到1454是最左边的一个剖面，5225到5514是最右边的一个剖面。由于CALIPSO的这几个剖面靠得很近，选取一个剖面即可。程序中就选取的最左边的剖面。

meng_xlei

灭火器 发表于 2019-7-31 09:43
如图，1165到1454是最左边的一个剖面，5225到5514是最右边的一个剖面。由于CALIPSO的这几个剖面靠得很近 ...

哇，学到了，厉害！！！
另外还有一个问题想请教您，用您的程序画了同一轨道、不同纬度范围的两张图，一张范围为35度到42度，另一张范围为42度到48度，但两张图完全相同，这是什么原因呢？

freekiller

谢谢分享，图文并茂

柒HHH

你好，如果想画的是沿卫星轨迹来画的话，用ncl代码应该怎么实现呢？

灭火器

柒HHH 发表于 2019-9-22 21:16
你好，如果想画的是沿卫星轨迹来画的话，用ncl代码应该怎么实现呢？

hdfeos.org/zoo/index_openLaRC_Examples.php#CALIPSO
就照抄这个网址中给出的水平的例子，选取某个高度的bin中的一条ray来画

no.2

请教一下，５５１４表示什么意思？是高度吗？