用MeteoInfo绘制全球水色分布图

MeteoInfo

看到Aires发的用hdfview画全球水色分布图的帖子（http://bbs.06climate.com/forum.php?mod=viewthread&tid=10426&extra=&page=1），也尝试用MeteoInfo Java版做了一下。数据：S19972441997273.L3m_MO_CHL_chlor_a_9km，下载地址：http://oceandata.sci.gsfc.nasa.gov/SeaWiFS/Mapped/Monthly/9km/chlor/ 。

数据是HDF4格式，没有遵循EOS规范，因此MeteoInfo Java版不能自动读取格点位置信息，不过用软件打开文件后可以查看文件信息：
File Name: D:\Temp\Hdf\S19972441997273.L3m_MO_CHL_chlor_a_9km
File type: Hierarchical Data Format, version 4 (HDF4)
Dimensions: 2
        fakeDim0 = 2160;
        fakeDim1 = 4320;
Global Attributes: [Product_Name = "S19972441997273.L3m_MO_CHL_chlor_a_9km", Sensor_Name = "SeaWiFS", Sensor = "Sea-viewing Wide Field-of-view Sensor (SeaWiFS)", Title = "SeaWiFS Level-3 Standard Mapped Image", Data_Center = "NASA/GSFC SeaWiFS Data Processing Center", Station_Name = "Wallops Flight Facility", Station_Latitude = 37.9272f, Station_Longitude = -75.4753f, Mission = "SeaStar SeaWiFS", Mission_Characteristics = "Nominal orbit: inclination = 98.2 (Sun-synchronous); node = 12 noon local (descending); eccentricity = <0.002; altitude = 705 km; ground speed = 6.75 km/sec", Sensor_Characteristics = "Number of bands = 8; number of active bands = 8; wavelengths per band (nm) = 412, 443, 490, 510, 555, 670, 765, 865; bits per pixel = 10; instantaneous field-of-view = 1.5835 mrad; pixels per scan = 1285; scan rate = 6/sec; sample rate = 7710/sec", Product_Type = "other", Processing_Version = "2010.0", Software_Name = "smigen", Software_Version = "4.17", Processing_Time = "2010272172320000", Input_Files = "S19972441997273.L3b_MO_CHL.main", Processing_Control = "smigen par=S19972441997273.L3m_MO_CHL_chlor_a_9km.param", Input_Parameters = "IFILE = /data4/sdpsoper/vdc/vpu3/workbuf/S19972441997273.L3b_MO_CHL.main|OFILE = S19972441997273.L3m_MO_CHL_chlor_a_9km|PFILE = |PROD = chlor_a|PALFILE = DEFAULT|PROCESSING VERSION = 2010.0|MEAS = 1|STYPE = 0|DATAMIN = 0.000000|DATAMAX = 0.000000|LONWEST = -180.000000|LONEAST = 180.000000|LATSOUTH = -90.000000|LATNORTH = 90.000000|RESOLUTION = 9km|PROJECTION = RECT|GAP_FILL = 0|SEAM_LON = -180.000000|PRECISION=F", L2_Flag_Names = "ATMFAIL,LAND,HILT,HISATZEN,STRAYLIGHT,CLDICE,COCCOLITH,LOWLW,CHLWARN,CHLFAIL,NAVWARN,MAXAERITER,ATMWARN,HISOLZEN,NAVFAIL,FILTER,HIGLINT", Period_Start_Year = 1997S, Period_Start_Day = 247S, Period_End_Year = 1997S, Period_End_Day = 273S, Start_Time = "1997247162633871", End_Time = "1997274055157036", Start_Year = 1997S, Start_Day = 247S, Start_Millisec = 59193871, End_Year = 1997S, End_Day = 274S, End_Millisec = 21117036, Start_Orbit = 519, End_Orbit = 905, Orbit = 712, Map_Projection = "Equidistant Cylindrical", Latitude_Units = "degrees North", Longitude_Units = "degrees East", Northernmost_Latitude = 90.0f, Southernmost_Latitude = -90.0f, Westernmost_Longitude = -180.0f, Easternmost_Longitude = 180.0f, Latitude_Step = 0.083333336f, Longitude_Step = 0.083333336f, SW_Point_Latitude = -89.958336f, SW_Point_Longitude = -179.95833f, Data_Bins = 2750226, Number_of_Lines = 2160, Number_of_Columns = 4320, Parameter = "Chlorophyll a concentration", Measure = "Mean", Units = "mg m^-3", Scaling = "linear", Scaling_Equation = "(Slope*l3m_data) + Intercept = Parameter value", Slope = 1.0f, Intercept = 0.0f, Data_Minimum = 0.00149f, Data_Maximum = 97.93204f, Suggested_Image_Scaling_Minimum = 0.01f, Suggested_Image_Scaling_Maximum = 20.0f, Suggested_Image_Scaling_Type = "LOG", Suggested_Image_Scaling_Applied = "No", _History = "Direct read of HDF4 file through CDM library", HDF4_Version = "4.2.5 (HDF Version 4.2 Release 5, February 17, 2010)"]
        : Product_Name = "S19972441997273.L3m_MO_CHL_chlor_a_9km"
        : Sensor_Name = "SeaWiFS"
        : Sensor = "Sea-viewing Wide Field-of-view Sensor (SeaWiFS)"
        : Title = "SeaWiFS Level-3 Standard Mapped Image"
        : Data_Center = "NASA/GSFC SeaWiFS Data Processing Center"
        : Station_Name = "Wallops Flight Facility"
        : Station_Latitude = 37.9272f
        : Station_Longitude = -75.4753f
        : Mission = "SeaStar SeaWiFS"
        : Mission_Characteristics = "Nominal orbit: inclination = 98.2 (Sun-synchronous); node = 12 noon local (descending); eccentricity = <0.002; altitude = 705 km; ground speed = 6.75 km/sec"
        : Sensor_Characteristics = "Number of bands = 8; number of active bands = 8; wavelengths per band (nm) = 412, 443, 490, 510, 555, 670, 765, 865; bits per pixel = 10; instantaneous field-of-view = 1.5835 mrad; pixels per scan = 1285; scan rate = 6/sec; sample rate = 7710/sec"
        : Product_Type = "other"
        : Processing_Version = "2010.0"
        : Software_Name = "smigen"
        : Software_Version = "4.17"
        : Processing_Time = "2010272172320000"
        : Input_Files = "S19972441997273.L3b_MO_CHL.main"
        : Processing_Control = "smigen par=S19972441997273.L3m_MO_CHL_chlor_a_9km.param"
        : Input_Parameters = "IFILE = /data4/sdpsoper/vdc/vpu3/workbuf/S19972441997273.L3b_MO_CHL.main|OFILE = S19972441997273.L3m_MO_CHL_chlor_a_9km|PFILE = |PROD = chlor_a|PALFILE = DEFAULT|PROCESSING VERSION = 2010.0|MEAS = 1|STYPE = 0|DATAMIN = 0.000000|DATAMAX = 0.000000|LONWEST = -180.000000|LONEAST = 180.000000|LATSOUTH = -90.000000|LATNORTH = 90.000000|RESOLUTION = 9km|PROJECTION = RECT|GAP_FILL = 0|SEAM_LON = -180.000000|PRECISION=F"
        : L2_Flag_Names = "ATMFAIL,LAND,HILT,HISATZEN,STRAYLIGHT,CLDICE,COCCOLITH,LOWLW,CHLWARN,CHLFAIL,NAVWARN,MAXAERITER,ATMWARN,HISOLZEN,NAVFAIL,FILTER,HIGLINT"
        : Period_Start_Year = 1997S
        : Period_Start_Day = 247S
        : Period_End_Year = 1997S
        : Period_End_Day = 273S
        : Start_Time = "1997247162633871"
        : End_Time = "1997274055157036"
        : Start_Year = 1997S
        : Start_Day = 247S
        : Start_Millisec = 59193871
        : End_Year = 1997S
        : End_Day = 274S
        : End_Millisec = 21117036
        : Start_Orbit = 519
        : End_Orbit = 905
        : Orbit = 712
        : Map_Projection = "Equidistant Cylindrical"
        : Latitude_Units = "degrees North"
        : Longitude_Units = "degrees East"
        : Northernmost_Latitude = 90.0f
        : Southernmost_Latitude = -90.0f
        : Westernmost_Longitude = -180.0f
        : Easternmost_Longitude = 180.0f
        : Latitude_Step = 0.083333336f
        : Longitude_Step = 0.083333336f
        : SW_Point_Latitude = -89.958336f
        : SW_Point_Longitude = -179.95833f
        : Data_Bins = 2750226
        : Number_of_Lines = 2160
        : Number_of_Columns = 4320
        : Parameter = "Chlorophyll a concentration"
        : Measure = "Mean"
        : Units = "mg m^-3"
        : Scaling = "linear"
        : Scaling_Equation = "(Slope*l3m_data) + Intercept = Parameter value"
        : Slope = 1.0f
        : Intercept = 0.0f
        : Data_Minimum = 0.00149f
        : Data_Maximum = 97.93204f
        : Suggested_Image_Scaling_Minimum = 0.01f
        : Suggested_Image_Scaling_Maximum = 20.0f
        : Suggested_Image_Scaling_Type = "LOG"
        : Suggested_Image_Scaling_Applied = "No"
        : _History = "Direct read of HDF4 file through CDM library"
        : HDF4_Version = "4.2.5 (HDF Version 4.2 Release 5, February 17, 2010)"
Variations: 1
        float l3m_data(fakeDim0,fakeDim1);
                l3m_data: Fill = -32767.0f
                l3m_data: Scaling = "linear"
                l3m_data: Scaling_Equation = "(Slope*l3m_data) + Intercept = Parameter value"
                l3m_data: Slope = 1.0f
                l3m_data: Intercept = 0.0f


文件信息里有格点投影、起始经纬度、经纬度间隔等信息，利用这些信息在MeteoInfo的脚本编辑器中写一个Jython脚本即可正确打开数据并生成数据图层。修改图例后显示如下：


脚本代码：

1. #-----------------------------------------------------
   
2. # Author: Yaqiang Wang
   
3. # Date: 2013-10-29
   
4. # Purpose: Read 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.projection import ProjectionInfo
    
12. from org.meteoinfo.legend import LegendManage
    
13. from org.meteoinfo.legend import LegendType
    
14. from org.meteoinfo.shape import ShapeTypes
    
15. import os.path
    
16. import jarray
    
17. 
    
18. #Set data directory
    
19. dataDir = 'D:/Temp/hdf/'
    
20. 
    
21. #Create MeteoDataInfo object
    
22. mdi = MeteoDataInfo()
    
23. 
    
24. #Read hdf data file
    
25. fn = dataDir + 'S19972441997273.L3m_MO_CHL_chlor_a_9km'
    
26. if os.path.isfile(fn):
    
27.         print fn
    
28.         mdi.openNetCDFData(fn)
    
29.         dataInfo = mdi.getDataInfo()
    
30.         print mdi.getInfoText()
    
31.         xmin = -180.0
    
32.         ymin = -90.0
    
33.         xnum = 4320
    
34.         ynum = 2160
    
35.         xdelt = 0.083333336
    
36.         ydelt = 0.083333336
    
37.         xlist = []
    
38.         ylist = []
    
39.         for i in range(0,xnum):
    
40.                 xlist.append(xmin + xdelt * i)
    
41.         for i in range(0,ynum):
    
42.                 ylist.append(ymin + ydelt * i)
    
43. 
    
44.         X = jarray.array(xlist, 'd')
    
45.         Y = jarray.array(ylist, 'd')
    
46.         xDim = Dimension(DimensionType.X)
    
47.         xDim.setValues(X)
    
48.         dataInfo.setXDimension(xDim)
    
49.         yDim = Dimension(DimensionType.Y)
    
50.         yDim.setValues(Y)
    
51.         dataInfo.setYDimension(yDim)
    
52.         var = dataInfo.getVariable('l3m_data')
    
53.         dimList = [yDim, xDim]
    
54.         var.setDimensions(dimList)
    
55.         gData = mdi.getGridData(var.getName())
    
56.         gData.yReverse()
    
57.         gData.missingValue = -32767.0
    
58.         aLS = LegendManage.createLegendSchemeFromGridData(gData, LegendType.GraduatedColor, ShapeTypes.Polygon)
    
59.         aLayer = DrawMeteoData.createRasterLayer(gData, "Test_HDF", aLS)
    
60.         aLayer.setProjInfo(dataInfo.getProjectionInfo())
    
61.         mf = miapp.getMapDocument().getActiveMapFrame()
    
62.         mf.addLayer(aLayer)
    
63.         mf.moveLayer(aLayer, 0)
    
64. 
    
65. print 'Finished!'

#-----------------------------------------------------<br /> # Author: Yaqiang Wang<br /> # Date: 2013-10-29<br /> # Purpose: Read 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.projection import ProjectionInfo<br /> from org.meteoinfo.legend import LegendManage<br /> from org.meteoinfo.legend import LegendType<br /> from org.meteoinfo.shape import ShapeTypes<br /> import os.path<br /> import jarray<br /> <br /> #Set data directory<br /> dataDir = 'D:/Temp/hdf/'<br /> <br /> #Create MeteoDataInfo object<br /> mdi = MeteoDataInfo()<br /> <br /> #Read hdf data file<br /> fn = dataDir + 'S19972441997273.L3m_MO_CHL_chlor_a_9km'<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 = 4320<br />         ynum = 2160<br />         xdelt = 0.083333336<br />         ydelt = 0.083333336<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('l3m_data')<br />         dimList = [yDim, xDim]<br />         var.setDimensions(dimList)<br />         gData = mdi.getGridData(var.getName())<br />         gData.yReverse()<br />         gData.missingValue = -32767.0<br />         aLS = LegendManage.createLegendSchemeFromGridData(gData, LegendType.GraduatedColor, ShapeTypes.Polygon)<br />         aLayer = DrawMeteoData.createRasterLayer(gData, "Test_HDF", aLS)<br />         aLayer.setProjInfo(dataInfo.getProjectionInfo())<br />         mf = miapp.getMapDocument().getActiveMapFrame()<br />         mf.addLayer(aLayer)<br />         mf.moveLayer(aLayer, 0)<br /> <br /> print 'Finished!'

Aires

膜拜王老师
小A的帖子能让王老师看到好开心