- 积分
- 26454
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2017-9-4
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
本帖最后由 灭火器 于 2019-4-18 08:22 编辑
读取数据主要参考了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下运行,具体如下:
- ;----------------------------------------------------------------
- ;# 2019/04/17
- ;# Draw CALIOP feature type and sub type between -0.5km to 20.2 km.
- ;# Reference:
- ; https://www-calipso.larc.nasa.gov/resources/calipso_users_guide/data_summaries/vfm/
- ;----------------------------------------------------------------
- begin
- filename = "../CALIOP/CAL_LID_L2_VFM-Standard-V4-20.2015-04-18T04-57-10ZD.hdf"
- file_vfm = addfile(filename, "r")
-
- ; Dim of Latitude is 0:4223, total 4224.
- latmin = 30.0
- latmax = 40.0
- lat_raw = file_vfm->Latitude(0:4223,0)
- ii = ind(lat_raw .ge. latmin .and. lat_raw .le. latmax)
- lat = lat_raw(ii)
-
- ; -0.5km to 8.2km, profile 0, subscripts are 1165:5514.
- fcf1 = file_vfm->Feature_Classification_Flags(ii,1165:1454)
- dims1 = dimsizes(fcf1)
- ; 8.2km to 20.2km, profile 0, subscripts are 165:364.
- fcf2 = file_vfm->Feature_Classification_Flags(ii, 165:364)
- dims2 = dimsizes(fcf2)
- ; Select 1-3 bits, Feature Type
- ; 0 = invalid
- ; 1 = clear air
- ; 2 = cloud
- ; 3 = aerosol
- ; 4 = stratospheric feature
- ; 5 = surface
- ; 6 = subsurface
- ; 7 = no signal
- fType1 = toint(dim_gbits(fcf1, 13, 3, 13, dims1(1)))
- subType1 = toint(dim_gbits(fcf1, 4, 3, 13, dims1(1)))
- fType2 = toint(dim_gbits(fcf2, 13, 3, 13, dims2(1)))
- subType2 = toint(dim_gbits(fcf2, 4, 3, 13, dims2(1)))
- ; When Feature type = 3 and subtype = 2, then type is dust, let it be 8.
- fType1 = where((fType1 .eq. 3 .and. subType1 .eq. 2), 8, fType1)
- fType2 = where((fType2 .eq. 3 .and. subType2 .eq. 2), 8, fType2)
- ; Reverse the height dim, set index 0 = surface.
- fType1 = fType1(:,::-1)
- fType2 = fType2(:,::-1)
- ; Use array data to store fType1 and fType2.
- ; Note that the vertical resolution of fType2 is twice as long as fType1.
- data = new((/dims1(0), dims1(1)+2*dims2(1)/), integer)
- data(:,0:dims1(1)-1) = fType1
- ; Because NCL draws plot according to the array size rather than coordinate
- ; size, so divide one bin in fType2 into 2 bins.
- do i = 0, dims2(1)-1
- data(:,dims1(1)+2*i) = fType2(:,i)
- data(:,dims1(1)+2*i+1) = fType2(:,i)
- end do
- ; Create altitude.
- ; Height is from -0.5km to 8.2km.
- dz1 = 30.0/1000.0 ; Units is km.
- hgt1 = (ispan(1, 290, 1)-0.5)*dz1 - 0.5
- ; Height2 is from 8.2km to 20.2km.
- dz2 = 60.0/1000.0 ; Units is km.
- hgt2 = (ispan(1, 400, 1)-0.5)*dz1 + 8.2
- hgt = new(690, float)
- hgt(0:289) = hgt1
- hgt(290:689) = hgt2
- hgt@long_name = "Height (km)"
- hgt@units = "km"
-
- ; Create latitude.
- lat@long_name = "Latitude"
- lat@units = "degrees_north"
-
- data!0 = "lat"
- data&lat = lat
- data!1 = "hgt"
- data&hgt = hgt
-
- ; Draw latitude VS height plot.
- picname = "./pic/VFM_v_dust"
- wks = gsn_open_wks("png", picname)
- res = True
- res@gsnDraw = False
- res@gsnFrame = False
- colormap = (/"white", "aliceblue", "lightskyblue", "gold", "red", \
- "seagreen", "palegreen", "black", "orange4"/)
- levels = (/0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5/)
- lbstring = (/"invalid", "clear air", "cloud", "aerosol", "stratospheric", \
- "surface", "subsurface", "no signal", "dust"/)
- res@gsnMaximize = True
- res@cnFillOn = True
- res@cnLinesOn = False
- res@cnLineLabelsOn = False
- res@cnFillMode = "RasterFill"
- res@cnMissingValFillColor = 0
- res@cnMissingValFillPattern = 0
- res@cnFillPalette = colormap
- res@cnLevelSelectionMode = "ExplicitLevels"
- res@cnLevels = levels
- res@tmYROn = False
- res@tmXTOn = False
- ; res@tmYLMode = "Explicit"
- ; res@tmYLValues = ispan(0, 20, 5)*1.0
- ; res@tmYLLabels = res@tmYLValues
-
- res@lbLabelAutoStride = True
- res@lbBoxLinesOn = False
- res@lbLabelAlignment = "BoxCenters"
- res@lbLabelStrings = lbstring
- res@lbBoxMinorExtentF = 0.2
- res@lbLabelFontHeightF = 0.009
-
- res@tiMainString = ""
- res@vpWidthF = 0.8
- res@vpHeightF = 0.2
-
- ; Transpose fType so x=lat, y=hgt.
- plot = gsn_csm_contour(wks, transpose(data), res)
- draw(plot)
- frame(wks)
- end
复制代码
|
-
-
评分
-
查看全部评分
|