MeteoInfoLab脚本示例：CAPLIPSO数据

MeteoInfo

CALIPSO卫星载激光雷达数据被广泛用于气候和大气物理大气环境的研究中，但其数据格式和处理过程都比较复杂，这里只是给出两个简单的数据处理绘图例子。

1. # Add file
   
2. fn = 'CAL_LID_L2_VFM-ValStage1-V3-02.2011-12-31T23-18-11ZD.hdf'
   
3. f = addfile('D:/Temp/hdf/' + fn)
   
4. 
   
5. # Read data
   
6. vname = 'Feature_Classification_Flags'
   
7. var = f[vname]
   
8. data = var[:,1256]
   
9. lon = f['Longitude'][:,0]
   
10. lat = f['Latitude'][:,0]
    
11. lon = lon[::10]
    
12. lat = lat[::10]
    
13. data = data[::10]
    
14. 
    
15. # Extract Feature Type only through bitmask.
    
16. data = data & 7
    
17. 
    
18. # Plot
    
19. axesm()
    
20. lworld = shaperead('D:/Temp/map/country1.shp')
    
21. geoshow(lworld, edgecolor='k')
    
22. levs = arange(8)
    
23. cols = [(0,0,0),(0,0,255),(255,255,0),(0,255,0),(255,0,0), \
    
24.     (200,100,255),(100,50,255),(127,127,127)]
    
25. ls = makesymbolspec('point', levels=levs, colors=cols)
    
26. layer = scatterm(lon, lat, data, size=5, edge=False, symbolspec=ls)
    
27. colorbar(layer, ticklabels=['invalid', 'clear', 'cloud', 'aerosol', \
    
28.     'strato', 'surface', 'subsurf', 'no signal'])
    
29. xlim(-180, 180)
    
30. ylim(-90, 90)
    
31. title([fn, 'Feature Type at Altitude = 2500m'])

# Add file<br /> fn = 'CAL_LID_L2_VFM-ValStage1-V3-02.2011-12-31T23-18-11ZD.hdf'<br /> f = addfile('D:/Temp/hdf/' + fn)<br /> <br /> # Read data<br /> vname = 'Feature_Classification_Flags'<br /> var = f[vname]<br /> data = var[:,1256]<br /> lon = f['Longitude'][:,0]<br /> lat = f['Latitude'][:,0]<br /> lon = lon[::10]<br /> lat = lat[::10]<br /> data = data[::10]<br /> <br /> # Extract Feature Type only through bitmask.<br /> data = data & 7<br /> <br /> # Plot<br /> axesm()<br /> lworld = shaperead('D:/Temp/map/country1.shp')<br /> geoshow(lworld, edgecolor='k')<br /> levs = arange(8)<br /> cols = [(0,0,0),(0,0,255),(255,255,0),(0,255,0),(255,0,0), \<br />     (200,100,255),(100,50,255),(127,127,127)]<br /> ls = makesymbolspec('point', levels=levs, colors=cols)<br /> layer = scatterm(lon, lat, data, size=5, edge=False, symbolspec=ls)<br /> colorbar(layer, ticklabels=['invalid', 'clear', 'cloud', 'aerosol', \<br />     'strato', 'surface', 'subsurf', 'no signal'])<br /> xlim(-180, 180)<br /> ylim(-90, 90)<br /> title([fn, 'Feature Type at Altitude = 2500m'])

1. # Add file
   
2. fn = 'D:/Temp/hdf/CAL_LID_L2_VFM-ValStage1-V3-02.2011-12-31T23-18-11ZD.hdf'
   
3. f = addfile(fn)
   
4. 
   
5. # Read data
   
6. vname = 'Feature_Classification_Flags'
   
7. var = f[vname]
   
8. data = var[:,:]
   
9. lat = f['Latitude'][:,0]
   
10. 
    
11. # Extract Feature Type only through bitmask.
    
12. data = data & 7
    
13. 
    
14. # Subset latitude values for the region of interest (40N to 62N).
    
15. lat = lat[3500:4000]
    
16. size = lat.shape[0]
    
17. 
    
18. data2d = data[3500:4000, 1165:]  # -0.5km to  8.2km
    
19. data3d = reshape(data2d, (size, 15, 290))
    
20. data = data3d[:,0,:]
    
21. 
    
22. # Focus on cloud (=2) data only.
    
23. data[data > 2] = 0
    
24. data[data < 2] = 0
    
25. data[data == 2] = 1
    
26. 
    
27. # Generate altitude data according to file specification [1].
    
28. alt = zeros(290)
    
29. # -0.5km to 8.2km
    
30. for j in range (0, 290):
    
31.     alt[j] = -0.5 + j*0.03
    
32. 
    
33. # Plot
    
34. levs = arange(2)
    
35. cols = ['w','b']
    
36. ls = makesymbolspec('image', levels=levs, colors=cols)
    
37. layer = imshow(lat, alt, rot90(data, 3), symbolspec=ls)
    
38. colorbar(layer, ticklabels=['Others','Cloud'])
    
39. basename = os.path.basename(fn)
    
40. title([basename, 'Feature Type (Bits 1-3) in Feature Classification Flag'])
    
41. xlabel('Latitude (degrees north)')
    
42. ylabel('Altitude (km)')

# Add file<br /> fn = 'D:/Temp/hdf/CAL_LID_L2_VFM-ValStage1-V3-02.2011-12-31T23-18-11ZD.hdf'<br /> f = addfile(fn)<br /> <br /> # Read data<br /> vname = 'Feature_Classification_Flags'<br /> var = f[vname]<br /> data = var[:,:]<br /> lat = f['Latitude'][:,0]<br /> <br /> # Extract Feature Type only through bitmask.<br /> data = data & 7<br /> <br /> # Subset latitude values for the region of interest (40N to 62N).<br /> lat = lat[3500:4000]<br /> size = lat.shape[0]<br /> <br /> data2d = data[3500:4000, 1165:]  # -0.5km to  8.2km<br /> data3d = reshape(data2d, (size, 15, 290))<br /> data = data3d[:,0,:]<br /> <br /> # Focus on cloud (=2) data only.<br /> data[data > 2] = 0<br /> data[data < 2] = 0<br /> data[data == 2] = 1<br /> <br /> # Generate altitude data according to file specification [1].<br /> alt = zeros(290)<br /> # -0.5km to 8.2km<br /> for j in range (0, 290):<br />     alt[j] = -0.5 + j*0.03<br /> <br /> # Plot<br /> levs = arange(2)<br /> cols = ['w','b']<br /> ls = makesymbolspec('image', levels=levs, colors=cols)<br /> layer = imshow(lat, alt, rot90(data, 3), symbolspec=ls)<br /> colorbar(layer, ticklabels=['Others','Cloud'])<br /> basename = os.path.basename(fn)<br /> title([basename, 'Feature Type (Bits 1-3) in Feature Classification Flag'])<br /> xlabel('Latitude (degrees north)')<br /> ylabel('Altitude (km)')

dhy

[img][/img]

freekiller

{:eb502:}{:eb502:}

3870383

王老师meteinfo能算大气视热源Q1嘛

MeteoInfo

3870383 发表于 2016-1-20 17:27
王老师meteinfo能算大气视热源Q1嘛

没做过这类计算，很难回答这个问题

sunflower$

官网有直接给出matlab和IDL 的程序 稍微改改就能用

MeteoInfo

sunflower$ 发表于 2016-1-21 10:47
官网有直接给出matlab和IDL 的程序 稍微改改就能用

你说的官网指的是哪个网址？

zdx_1982

赞一个，不错哦

citrushe

王老师，想请问一下，按照这个代码运行的时候，为什么会出现以下这种问题
>>> run script...
Traceback (most recent call last):
  File "<iostream>", line 14, in <module>
  File "E:\meteoinfo\MeteoInfo\pylib\mipylib\dimarray.py", line 143, in __getitem__
    raise IndexError()
IndexError

MeteoInfo

citrushe 发表于 2016-9-2 16:14
王老师，想请问一下，按照这个代码运行的时候，为什么会出现以下这种问题
>>> run script...
Traceback ( ...

你把你的脚本程序贴出来看看