MeteoInfoLab脚本示例：CloudSAT HDF Swath数据

MeteoInfo

读取CloudSAT HDF Swath数据，绘图分上下两部分，上面是时间和高度维的Radar Reflectivity Factor二维图，下面是卫星轨迹图。

示例程序：

1. # Add file
   
2. f = addfile('D:/Temp/hdf/2010128055614_21420_CS_2B-GEOPROF_GRANULE_P_R04_E03.hdf')
   
3. 
   
4. # Read data
   
5. vname = 'Radar_Reflectivity'
   
6. v_data = f[vname]
   
7. data = v_data[:,:]
   
8. v_height = f['Height']
   
9. height = v_height[0,:]
   
10. time = f['Profile_time'][:]
    
11. lon = f['Longitude'][:]
    
12. lat = f['Latitude'][:]
    
13. 
    
14. # Read attributes
    
15. long_name = v_data.attrvalue('long_name')[0]
    
16. scale_factor = v_data.attrvalue('factor')[0]
    
17. valid_min = v_data.attrvalue('valid_range')[0]
    
18. valid_max = v_data.attrvalue('valid_range')[1]
    
19. units = v_data.attrvalue('units')[0]
    
20. units_h = v_height.attrvalue('units')[0]
    
21. 
    
22. # Apply scale factor
    
23. valid_max = valid_max / scale_factor
    
24. valid_min = valid_min / scale_factor
    
25. data = data / scale_factor
    
26. data[data>valid_max] = nan
    
27. data[data<valid_min] = nan
    
28. data = transpose(data)
    
29. data = data[::-1,:]
    
30. 
    
31. # Make a split window plot
    
32. subplot(2, 1, 1)
    
33. 
    
34. # Contour the data
    
35. levs = arange(-38, 50, 2)
    
36. layer = imshow(time, height[::-1], data, levs)
    
37. colorbar(layer)
    
38. title('Radar Reflectivity Factor')
    
39. xlabel('Seconds since the start of the granule. (seconds)')
    
40. ylabel('Height (m)')
    
41. 
    
42. # The 2nd plot is the trajectory
    
43. subplot(2, 1, 2)
    
44. axesm(newaxes=False)
    
45. lworld = shaperead('D:/Temp/map/country1.shp')
    
46. geoshow(lworld, edgecolor='k')
    
47. plotm(lon, lat, '-b', linewidth=4)
    
48. #scatterm(lon, lat, lon, size=4, edge=False, facecolor='b')
    
49. scatterm(lon[0], lat[0], size=6, facecolor='r')
    
50. xlim(-180, 180)
    
51. ylim(-90, 90)
    
52. title('Trajectory of Flight Path (starting point in red)')

# Add file<br /> f = addfile('D:/Temp/hdf/2010128055614_21420_CS_2B-GEOPROF_GRANULE_P_R04_E03.hdf')<br /> <br /> # Read data<br /> vname = 'Radar_Reflectivity'<br /> v_data = f[vname]<br /> data = v_data[:,:]<br /> v_height = f['Height']<br /> height = v_height[0,:]<br /> time = f['Profile_time'][:]<br /> lon = f['Longitude'][:]<br /> lat = f['Latitude'][:]<br /> <br /> # Read attributes<br /> long_name = v_data.attrvalue('long_name')[0]<br /> scale_factor = v_data.attrvalue('factor')[0]<br /> valid_min = v_data.attrvalue('valid_range')[0]<br /> valid_max = v_data.attrvalue('valid_range')[1]<br /> units = v_data.attrvalue('units')[0]<br /> units_h = v_height.attrvalue('units')[0]<br /> <br /> # Apply scale factor<br /> valid_max = valid_max / scale_factor<br /> valid_min = valid_min / scale_factor<br /> data = data / scale_factor<br /> data[data>valid_max] = nan<br /> data[data<valid_min] = nan<br /> data = transpose(data)<br /> data = data[::-1,:]<br /> <br /> # Make a split window plot<br /> subplot(2, 1, 1)<br /> <br /> # Contour the data<br /> levs = arange(-38, 50, 2)<br /> layer = imshow(time, height[::-1], data, levs)<br /> colorbar(layer)<br /> title('Radar Reflectivity Factor')<br /> xlabel('Seconds since the start of the granule. (seconds)')<br /> ylabel('Height (m)')<br /> <br /> # The 2nd plot is the trajectory<br /> subplot(2, 1, 2)<br /> axesm(newaxes=False)<br /> lworld = shaperead('D:/Temp/map/country1.shp')<br /> geoshow(lworld, edgecolor='k')<br /> plotm(lon, lat, '-b', linewidth=4)<br /> #scatterm(lon, lat, lon, size=4, edge=False, facecolor='b')<br /> scatterm(lon[0], lat[0], size=6, facecolor='r')<br /> xlim(-180, 180)<br /> ylim(-90, 90)<br /> title('Trajectory of Flight Path (starting point in red)')

freekiller

赞王老师，和之前matlab的几乎一样，功能越来越强大了。把第一幅图的横坐标设置为经纬就更好了。

大花园

{:eb502:}

5xx

file:///C:/Users/cma5x/Desktop/捕获.PNG 请问下用这个脚本老提示这个是什么原因

MeteoInfo

5xx 发表于 2016-1-4 09:41
请问下用这个脚本老提示这个是什么原因

你修改过这个脚本吗？怎么修改的？

5xx

MeteoInfo 发表于 2016-1-4 09:58
你修改过这个脚本吗？怎么修改的？

我下载的也是这个反射率数据，他里面数据存放的格式不都应该一样的么？

MeteoInfo

5xx 发表于 2016-1-4 10:28
我下载的也是这个反射率数据，他里面数据存放的格式不都应该一样的么？

你运行脚本后在Console里输入：
>>> f

然后回车可以看到文件的详细信息。看看变量名还有属性之类的是不是和脚本不同？

597442388

。。。。。。。。楼主

5xx

MeteoInfo 发表于 2016-1-4 10:34
你运行脚本后在Console里输入：
>>> f

从读取属性那里开始报错，但是我属性名都是相同的

MeteoInfo

5xx 发表于 2016-1-4 11:26
从读取属性那里开始报错，但是我属性名都是相同的

要不你把hdf文件作为附件附上，我下载试试。