NCL读取TRMM2A12（TMI）资料画其中水成物的垂直分布

denglin19901030

按照2A12资料的介绍，理解了其中水成物应该是两个变量相乘的结果，主要提到的是cluster(nspecies*nfindex*nlayer*ncluster)与clusterScale(nspecies*npixel*nscan)相乘，但是写完脚本运行，总是提示超出的设定的范围，但是反复查过维数和大小没有问题，请教有大神做过这类数据处理吗？

脚本如下：load "$NCARG_ROOT/lib/ncarg/nclex/gsun/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
begin

; Read file.
file_name = "0720/2A12.20140718.94953.7.HDF"
hdf4_file=addfile(file_name, "r")

; Print metadata information.
print(hdf4_file)

; Read data to plot.
cluster=hdf4_file->cluster
clusterN=hdf4_file->clusterNumber
clusterScale=hdf4_file->clusterScale
freezingHI=hdf4_file->freezingHeightIndex
clusterNumber=byte2flt(clusterN)
freezingHeightIndex=byte2flt(freezingHI)

nlayer=28
nscan=2918
npixel=208
printVarSummary(cluster)
printVarSummary(clusterNumber)
printVarSummary(clusterScale)
printVarSummary(freezingHeightIndex)
; Read lat/lon data.
longitude=hdf4_file->Longitude
latitude=hdf4_file->Latitude
rainwater=new((/nscan,npixel,nlayer/),"float")
printVarSummary(rainwater)

do j=0,nscan-1
do i=0,npixel-1
do z=0,nlayer-1
a=freezingHeightIndex(j,i)
aa=floattoint(a)
b=floattoint(clusterNumber(j,i,2))
bb=floattoint(b)
if((aa.eq.157))
rainwater(z,j,i)=-9999
continue
end if
rainwater(j,i,z)=cluster(bb,z,aa,2)*clusterScale(j,i,2)
delete(a)
delete(b)
end do
end do
end do

;print(max(kk))
; Correct units to follow the CF conventions.
; In the HDF4 file, the attribute values are "degrees".
longitude@units = "degrees_east"
latitude@units = "degrees_north"
xwks=gsn_open_wks("png","rain")
gsn_define_colormap(xwks,"prcp_2")
;plts = new(3,"graphic")
res=True ; plot mods desired
res@gsnDraw = True
res@gsnFrame = True
res@cnFillOn=True ; enable contour fill
res@gsnMaximize=True; make plot large
res@gsnPaperOrientation =        "portrait" ; force portrait orientation
res@cnLinesOn=True ; turn off contour lines
res@cnLineLabelsOn       =  False; turn off contour line labels
res@gsnSpreadColors=True ; use the entire color spectrum
plot=gsn_contour(xwks,rainwater(3,:,:),res)
end

longlivehj

62行是哪一行啊？是rainwater(j,i,z)=cluster(bb,z,aa,2)*clusterScale(j,i,2)这一行吗？aa和bb都是从1开始的，作为数组下标要减1，估计是这个原因。

denglin19901030

longlivehj 发表于 2015-3-16 17:04
62行是哪一行啊？是rainwater(j,i,z)=cluster(bb,z,aa,2)*clusterScale(j,i,2)这一行吗？aa和bb都是从1开始 ...

大神~~~就是这一行，我试试看~

puyiliang

      朋友，最近碰巧要用到2A12的轨道卫星资料，但是由于平时都用NCL处理格点资料，轨道资料的信息实在不懂啊，nlayer=28是层数，（ nscan=2918，npixel=208这两个是什么意思实在不懂？）而且我需要用到某个区域的数据，比如说南海区域，那我应该怎样来选取位置呐？因为实在找不到经纬度信息？
抱歉，看到了一个相关帖子就赶紧问问，身边的人都没有用过这个轨道资料。先谢过了

denglin19901030

puyiliang 发表于 2015-3-16 21:32
朋友，最近碰巧要用到2A12的轨道卫星资料，但是由于平时都用NCL处理格点资料，轨道资料的信息实在不 ...

我也是才开始接触的，可以和你交流一下啊，但是不确定我的想法是不是正确，是这样的，2A12里面有很多数据，比如你是想要surfacerain/surfaceprecipitation这种，它是二维的，但是你可以直接叠加地理信息，就比如经度纬度：
longitude=hdf4_file->Longitude
latitude=hdf4_file->Latitude
longitude@units = "degrees_east"
latitude@units = "degrees_north"
data@lon2d=longitude
data@lat2d=latitude

这样你画的时候，用plts=gsn_csm_contour_map_ce(xwks,data,res)，就可以画出有全球范围的带有经纬度的图了。
还有一种数据就是水成物的，里面的变量是cluster，但是这个数据处理起来有点麻烦，我也还在试验中~~~
不知道能不能帮到你~

puyiliang

denglin19901030 发表于 2015-3-17 08:32
我也是才开始接触的，可以和你交流一下啊，但是不确定我的想法是不是正确，是这样的，2A12里面有很多数据 ...

很感谢指点，我也慢慢揣摩

粲粲倾城

denglin19901030 发表于 2015-3-17 08:32
我也是才开始接触的，可以和你交流一下啊，但是不确定我的想法是不是正确，是这样的，2A12里面有很多数据 ...

                               
登录/注册后可看大图

如果是画某一点的垂直廓线，如上图，怎么在程序中确定这一点的经纬度？

粲粲倾城

denglin19901030 发表于 2015-3-17 08:32
我也是才开始接触的，可以和你交流一下啊，但是不确定我的想法是不是正确，是这样的，2A12里面有很多数据 ...

还有，V6版本水成物是包括云水云冰可降水可降冰四个物理量，V7版本的要用哪几个物理量来计算？

denglin19901030

puyiliang 发表于 2015-3-17 09:55
很感谢指点，我也慢慢揣摩

恩，我还没有尝试成功，要是你成功的话可以教教我哈

denglin19901030

粲粲倾城 发表于 2015-3-17 10:13
如果是画某一点的垂直廓线，如上图，怎么在程序中确定这一点的经纬度？

我是下载了一个叫orbitviewer的hdf文件阅读器，打开之后如下图，可以任意点击某个点，会出来对应的经纬度和nscan和npixel的数值，但是latentheat这个在2A12V7版本里是不能直接读的，要像我上面读水成物一样，其中的nspecis=5就是latentheat.