绘制MODIS AOD数据

MeteoInfo

MODIS AOD（气溶胶光学厚度）数据应用比较广泛，下载的数据格式有HDF, NetCDF等，NetCDF格式数据可以直接在MeteoInfo中打开绘制出正确的图形，HDF数据由于没有遵循EOS协议，需要写脚本来确定坐标信息以便正确绘图，这里给出一个例子：


脚本程序代码：

1. #-----------------------------------------------------
   
2. # Author: Yaqiang Wang
   
3. # Date: 2013-12-5
   
4. # Purpose: Read MODIS AOD hdf4 data
   
5. # Note: Sample
   
6. #-----------------------------------------------------
   
7. from org.meteoinfo.data.meteodata import MeteoDataInfo
   
8. from org.meteoinfo.data.meteodata import Dimension
   
9. from org.meteoinfo.data.meteodata import DimensionType
   
10. from org.meteoinfo.data.meteodata import DrawMeteoData
    
11. from org.meteoinfo.legend import LegendScheme
    
12. from org.meteoinfo.shape import ShapeTypes
    
13. from ucar.nc2 import NetcdfFile
    
14. import os.path
    
15. import jarray
    
16. 
    
17. #Set data directory
    
18. dataDir = 'D:/Temp/toMingJing/AOD/'
    
19. 
    
20. #Create MeteoDataInfo object
    
21. mdi = MeteoDataInfo()
    
22. 
    
23. #Read hdf data file
    
24. fn = os.path.join(dataDir, 'MOD08_M3.A2000061.051.2010273213159.pscs_000500501908.Optical_Depth_Land_And_Ocean_Mean_Mean.G3.HDF')
    
25. if os.path.isfile(fn):
    
26.         print fn
    
27.         mdi.openNetCDFData(fn)
    
28.         dataInfo = mdi.getDataInfo()
    
29.         #print mdi.getInfoText()
    
30.         xmin = -180.0
    
31.         ymin = -90.0
    
32.         xnum = 360
    
33.         ynum = 180
    
34.         xdelt = 1.0
    
35.         ydelt = 1.0
    
36.         xlist = []
    
37.         ylist = []
    
38.         for i in range(0,xnum):
    
39.                 xlist.append(xmin + xdelt * i)
    
40.         for i in range(0,ynum):
    
41.                 ylist.append(ymin + ydelt * i)
    
42. 
    
43.         X = jarray.array(xlist, 'd')
    
44.         Y = jarray.array(ylist, 'd')
    
45.         xDim = Dimension(DimensionType.X)
    
46.         xDim.setValues(X)
    
47.         dataInfo.setXDimension(xDim)
    
48.         yDim = Dimension(DimensionType.Y)
    
49.         yDim.setValues(Y)
    
50.         dataInfo.setYDimension(yDim)
    
51.         var = dataInfo.getVariable('Optical_Depth_Land_And_Ocean_Mean_Mean')
    
52.         print var.getName()
    
53.         dimList = [yDim, xDim]
    
54.         var.setDimensions(dimList)
    
55.         gData = mdi.getGridData(var.getName())       
    
56.         #aLS = LegendManage.createLegendSchemeFromGridData(gData, LegendType.GraduatedColor, ShapeTypes.Polygon)
    
57.         aLS = LegendScheme(ShapeTypes.Polygon)
    
58.         aLS.importFromXMLFile(os.path.join(dataDir, 'AOD.lgs'))
    
59.         aLayer = DrawMeteoData.createRasterLayer(gData, "Test_HDF", aLS)
    
60.         mf = miapp.getMapDocument().getActiveMapFrame()
    
61.         mf.addLayer(aLayer)
    
62.         mf.moveLayer(aLayer, 0)
    
63. 
    
64. print 'Finished!'

#-----------------------------------------------------<br /> # Author: Yaqiang Wang<br /> # Date: 2013-12-5<br /> # Purpose: Read MODIS AOD hdf4 data<br /> # Note: Sample<br /> #-----------------------------------------------------<br /> from org.meteoinfo.data.meteodata import MeteoDataInfo<br /> from org.meteoinfo.data.meteodata import Dimension<br /> from org.meteoinfo.data.meteodata import DimensionType<br /> from org.meteoinfo.data.meteodata import DrawMeteoData<br /> from org.meteoinfo.legend import LegendScheme<br /> from org.meteoinfo.shape import ShapeTypes<br /> from ucar.nc2 import NetcdfFile<br /> import os.path<br /> import jarray<br /> <br /> #Set data directory<br /> dataDir = 'D:/Temp/toMingJing/AOD/'<br /> <br /> #Create MeteoDataInfo object<br /> mdi = MeteoDataInfo()<br /> <br /> #Read hdf data file<br /> fn = os.path.join(dataDir, 'MOD08_M3.A2000061.051.2010273213159.pscs_000500501908.Optical_Depth_Land_And_Ocean_Mean_Mean.G3.HDF')<br /> if os.path.isfile(fn):<br />         print fn<br />         mdi.openNetCDFData(fn)<br />         dataInfo = mdi.getDataInfo()<br />         #print mdi.getInfoText()<br />         xmin = -180.0<br />         ymin = -90.0<br />         xnum = 360<br />         ynum = 180<br />         xdelt = 1.0<br />         ydelt = 1.0<br />         xlist = []<br />         ylist = []<br />         for i in range(0,xnum):<br />                 xlist.append(xmin + xdelt * i)<br />         for i in range(0,ynum):<br />                 ylist.append(ymin + ydelt * i)<br /> <br />         X = jarray.array(xlist, 'd')<br />         Y = jarray.array(ylist, 'd')<br />         xDim = Dimension(DimensionType.X)<br />         xDim.setValues(X)<br />         dataInfo.setXDimension(xDim)<br />         yDim = Dimension(DimensionType.Y)<br />         yDim.setValues(Y)<br />         dataInfo.setYDimension(yDim)<br />         var = dataInfo.getVariable('Optical_Depth_Land_And_Ocean_Mean_Mean')<br />         print var.getName()<br />         dimList = [yDim, xDim]<br />         var.setDimensions(dimList)<br />         gData = mdi.getGridData(var.getName())        <br />         #aLS = LegendManage.createLegendSchemeFromGridData(gData, LegendType.GraduatedColor, ShapeTypes.Polygon)<br />         aLS = LegendScheme(ShapeTypes.Polygon)<br />         aLS.importFromXMLFile(os.path.join(dataDir, 'AOD.lgs'))<br />         aLayer = DrawMeteoData.createRasterLayer(gData, "Test_HDF", aLS)<br />         mf = miapp.getMapDocument().getActiveMapFrame()<br />         mf.addLayer(aLayer)<br />         mf.moveLayer(aLayer, 0)<br /> <br /> print 'Finished!'

rceclx

王老师的新作，强大啊！持续关注和学习！
周末愉快！

zihan6688

继续支持

rceclx

王老师，周末好！按照您的方法试了一下，可以读取HDF格式的AOD数据，但也遇到了几个问题：1、java版里的显示-屏蔽外部图形打开后，没有可供选择的shp文件，图层中已经添加了shp；2、java版选中hdf图层后，显示-图层属性数据不能显示，shp图层可以；3、AOD.lgs文件方便能提供一下么？4、C#版能不能读取HDF格式的AOD数据？谢谢您。

MeteoInfo

rceclx 发表于 2013-12-7 10:50

                               
登录/注册后可看大图

王老师，周末好！按照您的方法试了一下，可以读取HDF格式的AOD数据，但也遇到了几个问题：1、java版里的显示 ...

1、屏蔽用的图层必须是多边形图层，如果你要屏蔽中国以外的区域可以用china.shp
2、上面的例子生成的是RasterLayer（栅格图层），是没有属性数据表的。
3、 AOD.lgs (1.6 KB, 下载次数: 46)

2013-12-7 11:09 上传

点击文件名下载附件 4、C#版只能读取HDF5格式数据，这个例子的数据是HDF4格式，C#版读不了。

rceclx

MeteoInfo 发表于 2013-12-7 11:10

                               
登录/注册后可看大图

1、屏蔽用的图层必须是多边形图层，如果你要屏蔽中国以外的区域可以用china.shp
2、上面的例子生成的是R ...

王老师，非常感谢！
我就是准备屏蔽中国以外的区域。还有点问题，1、例子中生成的RasterLayer能输出么，不只是画图显示出来，想按照“也说格点数据插值为站点数据
http://bbs.06climate.com/forum.p ... 91&fromuid=5721”的方法，将格点数据插值为站点数据？2、在添加了图层后，meteo data打开后不能添加数据了，但重新打开java版，没有任何图层的时候，就没问题。

MeteoInfo

rceclx 发表于 2013-12-7 11:29

                               
登录/注册后可看大图

王老师，非常感谢！
我就是准备屏蔽中国以外的区域。还有点问题，1、例子中生成的RasterLayer能 ...

1、可以用脚本程序来进行插值
或者在上面的脚本程序里将获取的GridData (gData)存为Surfer ASCII 格点数据，然后在GUI中操作。

2、可能是Java Swing的bug，我也偶尔发现类似的现象，不知道问题在哪里？

rceclx

MeteoInfo 发表于 2013-12-7 12:43

                               
登录/注册后可看大图

1、可以用脚本程序来进行插值
或者在上面的脚本程序里将获取的GridData (gData)存为Surfer ASCII 格点数 ...

好的，我按您的提示先试验一下，有问题再请教王老师！

rceclx

MeteoInfo 发表于 2013-12-7 12:43

                               
登录/注册后可看大图

1、可以用脚本程序来进行插值
或者在上面的脚本程序里将获取的GridData (gData)存为Surfer ASCII 格点数 ...

王老师，试着修改了一下脚本，将获取的GridData (gData)存为Surfer ASCII 格点数据。
把C#中版中用过的outfn = dataDir + "test.grd"
gData.SaveAsSurferASCIIFile(outfn)，放到 java版的脚本里，提示属性错误，不能SaveAsSurferASCIIFile(outfn)，请问一下java中这个命令应该怎么写？有没有可供查找此类命令的指南？谢谢。

MeteoInfo

rceclx 发表于 2013-12-7 16:52

                               
登录/注册后可看大图

王老师，试着修改了一下脚本，将获取的GridData (gData)存为Surfer ASCII 格点数据。
把C#中版中用过的o ...

SaveAsSurferASCIIFile第一个字母小写saveAsSurferASCIIFile

可以下载一个jd-gui来查看MeteoInfoLib.jar文件中的类和方法等，源代码也可以看到呦。