MeteoInfoLab脚本示例：SeaWiFS HDF Grid数据

MeteoInfo

SeaWiFS HDF Grid数据读取，特别是涉及到了文件的众多属性数据的读取，数据取对数后绘图。

脚本程序：

1. #Add data file
   
2. f = addfile('D:/Temp/hdf/S1999001.L3m_DAY_CHL_chlor_a_9km.hdf')
   
3. #Get data variable
   
4. vname = 'l3m_data'
   
5. v = f[vname]
   
6. #Set x/y
   
7. ny = f.attrvalue('Number_of_Lines')[0]
   
8. nx = f.attrvalue('Number_of_Columns')[0]
   
9. sx = f.attrvalue('Westernmost_Longitude')[0]
   
10. ex = f.attrvalue('Easternmost_Longitude')[0]
    
11. sy = f.attrvalue('Southernmost_Latitude')[0]
    
12. ey = f.attrvalue('Northernmost_Latitude')[0]
    
13. x = linspace(sx, ex, nx)
    
14. y = linspace(sy, ey, ny)
    
15. #Set x/y dimensions
    
16. v.setdim('Y', y, 0)
    
17. v.setdim('X', x, 1)
    
18. #Get data array
    
19. fillv = v.attrvalue('Fill')[0]
    
20. scale = v.attrvalue('Slope')[0]
    
21. offset = v.attrvalue('Intercept')[0]
    
22. data = v[::-1,:] * scale + offset
    
23. data.fill_value = fillv
    
24. data = log(data)
    
25. #Plot
    
26. axesm()
    
27. world = shaperead('D:/Temp/map/country1.shp')
    
28. geoshow(world)
    
29. levs = arange(-6, 4, 0.2)
    
30. layer = imshowm(data, levs)
    
31. colorbar(layer)
    
32. title('SeaWiFS - ' + vname + ' (log scale)')

#Add data file<br /> f = addfile('D:/Temp/hdf/S1999001.L3m_DAY_CHL_chlor_a_9km.hdf')<br /> #Get data variable<br /> vname = 'l3m_data'<br /> v = f[vname]<br /> #Set x/y<br /> ny = f.attrvalue('Number_of_Lines')[0]<br /> nx = f.attrvalue('Number_of_Columns')[0]<br /> sx = f.attrvalue('Westernmost_Longitude')[0]<br /> ex = f.attrvalue('Easternmost_Longitude')[0]<br /> sy = f.attrvalue('Southernmost_Latitude')[0]<br /> ey = f.attrvalue('Northernmost_Latitude')[0]<br /> x = linspace(sx, ex, nx)<br /> y = linspace(sy, ey, ny)<br /> #Set x/y dimensions<br /> v.setdim('Y', y, 0)<br /> v.setdim('X', x, 1)<br /> #Get data array<br /> fillv = v.attrvalue('Fill')[0]<br /> scale = v.attrvalue('Slope')[0]<br /> offset = v.attrvalue('Intercept')[0]<br /> data = v[::-1,:] * scale + offset<br /> data.fill_value = fillv<br /> data = log(data)<br /> #Plot<br /> axesm()<br /> world = shaperead('D:/Temp/map/country1.shp')<br /> geoshow(world)<br /> levs = arange(-6, 4, 0.2)<br /> layer = imshowm(data, levs)<br /> colorbar(layer)<br /> title('SeaWiFS - ' + vname + ' (log scale)')