MeteoInfo脚本示例：读取FY3A AOD HDF文件

FY3A卫星有AOD产品数据，HDF格式，这里示例用MeteoInfo脚本程序读取和显示该类数据。

 

脚本程序如下：

#-----------------------------------------------------
# Author: Yaqiang Wang
# Date: 2015-3-18
# Purpose: Read FY3A AOD hdf5 data
# Note: Sample
#-----------------------------------------------------
from org.meteoinfo.data.meteodata import MeteoDataInfo
from org.meteoinfo.data.meteodata import Dimension
from org.meteoinfo.data.meteodata import DimensionType
from org.meteoinfo.data.meteodata import DrawMeteoData
from org.meteoinfo.legend import LegendScheme
from org.meteoinfo.shape import ShapeTypes
from ucar.nc2 import NetcdfFile
import os.path
import jarray

#Set data directory
dataDir = 'D:/Temp/Hdf'

#Create MeteoDataInfo object
mdi = MeteoDataInfo()

#Read hdf data file
fn = os.path.join(dataDir, 'FY3A_MERSI_GBAL_L3_ASL_MLT_GLL_20140331_AOAM_5000M_MS.HDF')
if os.path.isfile(fn):
    print fn
    mdi.openNetCDFData(fn)
    dataInfo = mdi.getDataInfo()
    #print mdi.getInfoText()
    xmin = -30.0
    ymin = -90.0
    xnum = 7200
    ynum = 3600
    xdelt = 0.05
    ydelt = 0.05
    xlist = []
    ylist = []
    for i in range(0,xnum):
        xlist.append(xmin + xdelt * i)
    for i in range(0,ynum):
        ylist.append(ymin + ydelt * i)

    X = jarray.array(xlist, 'd')
    Y = jarray.array(ylist, 'd')
    xDim = Dimension(DimensionType.X)
    xDim.setValues(X)
    dataInfo.setXDimension(xDim)
    yDim = Dimension(DimensionType.Y)
    yDim.setValues(Y)
    dataInfo.setYDimension(yDim)
    var = dataInfo.getVariable('Aerosol_Optical_Thickness_of_MERSI_550nm')
    print var.getName()
    dimList = [yDim, xDim]
    var.setDimensions(dimList)
    gData = mdi.getGridData(var.getName())        
    gData.yReverse()
    gData = gData.div(10000.0)
    gData.missingValue = -3.2767
    #aLS = LegendManage.createLegendSchemeFromGridData(gData, LegendType.GraduatedColor, ShapeTypes.Polygon)
    aLS = LegendScheme(ShapeTypes.Polygon)
    aLS.importFromXMLFile(os.path.join(dataDir, 'AOD.lgs'))
    aLayer = DrawMeteoData.createRasterLayer(gData, "Test_HDF", aLS)
    mf = miapp.getMapDocument().getActiveMapFrame()
    mf.addLayer(aLayer)
    mf.moveLayer(aLayer, 0)

print 'Finished!'

 

转载于:https://www.cnblogs.com/yaqiang/p/4347538.html