处理FY2E圆盘标称投影数据

MeteoInfo

圆盘标称投影数据时静止气象卫星常见的数据产品，比如FY2E静止气象卫星就有很多这样的产品（可以从国家卫星气象中心网站上下载）。所谓的圆盘标称投影就是Geostationary投影，主要的投影参数有中央经度和卫星高度。这里给出一个示例脚本程序可以读取FY2E HDF格式数据、将数据投影为等经纬度并生成数据图层。需要最新的MeteoInfo Java版本（MeteoInfo Java 1.1.6R1，http://www.meteothinker.com）。

FY2E的HDF数据自描述得很粗糙，数据的范围（脚本中的xmin等）是自己调试出来的，如果有人知道确切的数据也请指出。

脚本程序：

1. #-----------------------------------------------------# Author: Yaqiang Wang
   
2. # Date: 2014-11-16
   
3. # Purpose: Read FY2E hdf5 data
   
4. # Note: Sample
   
5. #-----------------------------------------------------
   
6. from org.meteoinfo.data.meteodata import MeteoDataInfo
   
7. from org.meteoinfo.data.meteodata import Dimension
   
8. from org.meteoinfo.data.meteodata import DimensionType
   
9. from org.meteoinfo.data.meteodata import DrawMeteoData
   
10. from org.meteoinfo.projection import ProjectionInfo
    
11. from org.meteoinfo.projection import KnownCoordinateSystems
    
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. projInfo = ProjectionInfo("+proj=geos +lon_0=104.5 +h=35785864")
    
25. 
    
26. #Read hdf data file
    
27. #fn = dataDir + 'FY2E_CTA_MLT_NOM_20141116_0200.hdf'
    
28. fn = dataDir + 'FY2E_FDI_ALL_NOM_20110713_0100.hdf'
    
29. if os.path.isfile(fn):
    
30.         print fn
    
31.         mdi.openNetCDFData(fn)
    
32.         dataInfo = mdi.getDataInfo()
    
33.         print mdi.getInfoText()
    
34.         xmin = -5750000.0
    
35.         ymin = -5750000.0
    
36.         xnum = 2288
    
37.         ynum = 2288
    
38.         xmax = 5750000.0
    
39.         ymax = 5750000.0
    
40.         xdelt = (xmax - xmin) / xnum
    
41.         ydelt = (ymax - ymin) / ynum
    
42.         xlist = []
    
43.         ylist = []
    
44.         for i in range(0,xnum):
    
45.                 xlist.append(xmin + xdelt * i)
    
46.         for i in range(0,ynum):
    
47.                 ylist.append(ymin + ydelt * i)
    
48. 
    
49.         X = jarray.array(xlist, 'd')
    
50.         Y = jarray.array(ylist, 'd')
    
51.         xDim = Dimension(DimensionType.X)
    
52.         xDim.setValues(X)
    
53.         dataInfo.setXDimension(xDim)
    
54.         yDim = Dimension(DimensionType.Y)
    
55.         yDim.setValues(Y)
    
56.         dataInfo.setYDimension(yDim)
    
57.         #var = dataInfo.getVariable('FY2E_CLM_Hourly_Cloud_Amount')
    
58.         var = dataInfo.getVariable('NOMChannelVIS')
    
59.         dimList = [yDim, xDim]
    
60.         var.setDimensions(dimList)
    
61.         print 'Get grid data...'
    
62.         gData = mdi.getGridData(var.getName())        
    
63.         gData.yReverse()
    
64.         gData.projInfo = projInfo
    
65.         #gData.missingValue = 0.0
    
66.         print 'Reproject grid data to Lon/Lat...'
    
67.         toProj = KnownCoordinateSystems.geographic.world.WGS1984
    
68.         gData = gData.project(toProj)
    
69.         print 'Ouput grid data to Sufer ASCII grid data file...'
    
70.         outfn = 'D:/Temp/test/NOMChannelVIS.grd'
    
71.         gData.saveAsSurferASCIIFile(outfn)
    
72.         print 'Create raster layer from the grid data...'
    
73.         aLS = LegendManage.createLegendSchemeFromGridData(gData, LegendType.GraduatedColor, ShapeTypes.Polygon)
    
74.         aLayer = DrawMeteoData.createRasterLayer(gData, "Test_FY2E", aLS)
    
75.         #aLayer = DrawMeteoData.createShadedLayer(gData, aLS, 'Test_HDF_Shaded', 'AOD', True)
    
76.         aLayer.setProjInfo(toProj)
    
77.         #aLayer.setProjInfo(projInfo)
    
78.         print 'Add the layer...'
    
79.         mf = miapp.getMapDocument().getActiveMapFrame()
    
80.         mf.addLayer(aLayer)
    
81.         mf.moveLayer(aLayer, 0)
    
82. 
    
83. print 'Finished!'

#-----------------------------------------------------# Author: Yaqiang Wang<br /> # Date: 2014-11-16<br /> # Purpose: Read FY2E hdf5 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.projection import KnownCoordinateSystems<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 /> projInfo = ProjectionInfo("+proj=geos +lon_0=104.5 +h=35785864")<br /> <br /> #Read hdf data file<br /> #fn = dataDir + 'FY2E_CTA_MLT_NOM_20141116_0200.hdf'<br /> fn = dataDir + 'FY2E_FDI_ALL_NOM_20110713_0100.hdf'<br /> if os.path.isfile(fn):<br />         print fn<br />         mdi.openNetCDFData(fn)<br />         dataInfo = mdi.getDataInfo()<br />         print mdi.getInfoText()<br />         xmin = -5750000.0<br />         ymin = -5750000.0<br />         xnum = 2288<br />         ynum = 2288<br />         xmax = 5750000.0<br />         ymax = 5750000.0<br />         xdelt = (xmax - xmin) / xnum<br />         ydelt = (ymax - ymin) / ynum<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('FY2E_CLM_Hourly_Cloud_Amount')<br />         var = dataInfo.getVariable('NOMChannelVIS')<br />         dimList = [yDim, xDim]<br />         var.setDimensions(dimList)<br />         print 'Get grid data...'<br />         gData = mdi.getGridData(var.getName())        <br />         gData.yReverse()<br />         gData.projInfo = projInfo<br />         #gData.missingValue = 0.0<br />         print 'Reproject grid data to Lon/Lat...'<br />         toProj = KnownCoordinateSystems.geographic.world.WGS1984<br />         gData = gData.project(toProj)<br />         print 'Ouput grid data to Sufer ASCII grid data file...'<br />         outfn = 'D:/Temp/test/NOMChannelVIS.grd'<br />         gData.saveAsSurferASCIIFile(outfn)<br />         print 'Create raster layer from the grid data...'<br />         aLS = LegendManage.createLegendSchemeFromGridData(gData, LegendType.GraduatedColor, ShapeTypes.Polygon)<br />         aLayer = DrawMeteoData.createRasterLayer(gData, "Test_FY2E", aLS)<br />         #aLayer = DrawMeteoData.createShadedLayer(gData, aLS, 'Test_HDF_Shaded', 'AOD', True)<br />         aLayer.setProjInfo(toProj)<br />         #aLayer.setProjInfo(projInfo)<br />         print 'Add the layer...'<br />         mf = miapp.getMapDocument().getActiveMapFrame()<br />         mf.addLayer(aLayer)<br />         mf.moveLayer(aLayer, 0)<br /> <br /> print 'Finished!'

数据本身投影的图层：




投影为等经纬度数据：

river

高端 大气 上档次，好好学习一下

沙颖凯

脚本的语法看着像Python. 楼主的MeteoInfo有基于Python的部分吗？

MeteoInfo

沙颖凯 发表于 2014-11-17 03:52
脚本的语法看着像Python. 楼主的MeteoInfo有基于Python的部分吗？

脚本语言是Jython，Python的一个分支，用Java语言开发的，所以能方便的调用各类Java库。

沙颖凯

MeteoInfo 发表于 2014-11-16 18:20
脚本语言是Jython，Python的一个分支，用Java语言开发的，所以能方便的调用各类Java库。

楼主有CPython（我想说的是通常意义上的Python）的MeteoInfo模块吗，或者近期有做这件事的计划吗？我对其中数据I/O和一些数据处理的模块很感兴趣。

平流层的萝卜

Beautiful.希望有朝一日在全世界普及

MeteoInfo

沙颖凯 发表于 2014-11-17 14:56
楼主有CPython（我想说的是通常意义上的Python）的MeteoInfo模块吗，或者近期有做这件事的计划吗？我对其 ...

没有，这个不是想有就能有的，工作量巨大。

leikunjiang

王老师，这个数据有专门的经纬度文件 ，和数据的的大小长度是一样的，都是2288*2288个经度和纬度数据，这些经纬度数据就是数据集里每个数据在地球上的位置。我用micaps打开同一个卫星数据，发现你画的图和micaps上的是相反的。这个卫星数据是从西往东，从南往北的。也就是说有两个问题，一个就是数据范围的问题，标称投影有专门的经纬度数据；第二个就是数据图层弄反了

MeteoInfo

leikunjiang 发表于 2014-11-17 21:32
王老师，这个数据有专门的经纬度文件 ，和数据的的大小长度是一样的，都是2288*2288个经度和纬度数据，这些 ...

数据南北反向的问题已解决（见1楼帖子），在脚本中加一句：
gData.yReverse()

另一个问题“一个就是数据范围的问题，标称投影有专门的经纬度数据”，不知道你想表达什么意思？

沙颖凯

MeteoInfo 发表于 2014-11-17 02:03
没有，这个不是想有就能有的，工作量巨大。

我想我可以做一些这方面的事情，比如基于楼主的源代码写Python的模块之类，问下楼主的MeteoInfo是不是开源的，如果代码有托管在三方网站上就更好了。