处理FY2E圆盘标称投影数据

MeteoInfo

leikunjiang 发表于 2014-11-18 20:10
各个通道的分辨率确实是5km，不过2288乘以5000等于11440000，我不知道程序中的起始点-5717500是怎么得来 ...

应该是(2288 - 1) * 5000

leikunjiang

MeteoInfo 发表于 2014-11-18 21:52
应该是(2288 - 1) * 5000

王老师，我又发现另外一个问题，图片上的值与实际值(用micaps打开显示)是不正确的，var = dataInfo.getVariable('NOMChannelVIS')，这句话得到的数据应该是索引值（正整数（0到63的任意数）），要结合HDF里面的另外一个数据(CALChannelVIS)才是通道数据的值，比如读出'NOMChannelVIS里面的一个数据是5，则要在'CALChannelVIS'里找到5所对应的数据，这个5所对应的数据既为该点可见光的反射率值。这里'CALChannelVIS'只有64个值(0~63所对应的可将光通道数据)，就是要定标数据。我不知道有没有讲清楚。另外我想要那个由标称投影转换而成的等经纬度数据，脚本程序中有两句话outfn = 'E:\MeteInfo Java\MeteoInfo/NOMChannelVIS.grd'      和    gData.saveAsSurferASCIIFile(outfn)    这里的NOMChannelVIS.grd是不是等经纬度数据，还是那个生成的栅格数据才是等经纬度数据啊。

leikunjiang

刚刚经过检验，NOMChannelVIS.grd里的数据确实索引值，那么它就不是等经纬度投影......那转换的等经纬度数据如何得来呢？通常是由数据出图，那能不能反过来，由图得到那个等经纬度数据。主要是还要对那个等经纬度数据进行处理

MeteoInfo

leikunjiang 发表于 2014-11-20 21:34
刚刚经过检验，NOMChannelVIS.grd里的数据确实索引值，那么它就不是等经纬度投影......那转换的等经纬度数 ...

获取真实值可以用此脚本，具体看脚本里的代码吧，我就不详细说了：
代码贴出来有点小问题，把代码文件也附上： read_fy2e_hdf.py (2.95 KB, 下载次数: 35)

2014-11-21 14:28 上传

点击文件名下载附件

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

#-----------------------------------------------------<br /> # 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 />         #Get CALChannelVIS data<br />         calvis = mdi.getDataInfo().read('CALChannelVIS')<br />         cvs = []<br />         for i in range(0, calvis.getSize()):<br />                 cvs.append(calvis.getFloat(i))<br />         print cvs<br />         for i in range(0, ynum):<br />                 for j in range(0, xnum):<br />                         v = int(gData.data[j]<i>)<br />                         if v >= 0 and v < len(cvs):<br />                                 gData.data[j]<i> = cvs[v]<br />                                 <br />         print 'Ouput grid data to Sufer ASCII grid data file...'<br />         outfn = 'D:/Temp/test/NOMChannelVIS-1.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.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!'

leikunjiang

MeteoInfo 发表于 2014-11-21 14:27
获取真实值可以用此脚本，具体看脚本里的代码吧，我就不详细说了：
代码贴出来有点小问题，把代码文件也 ...

灰常灰常感谢王老师耐心细致的回答。我还想问问那个等经纬度数据能不能把它提取出来啊，因为还要对这个数据进行处理。程序代码中有两句话outfn = 'D:/Temp/test/NOMChannelVIS-1.grd'
                                             gData.saveAsSurferASCIIFile(outfn)
这两句话的意思是不是把它转换成了等经纬度了的啊

MeteoInfo

leikunjiang 发表于 2014-11-21 14:38
灰常灰常感谢王老师耐心细致的回答。我还想问问那个等经纬度数据能不能把它提取出来啊，因为还要对这个数 ...

已经转成了等经纬度数据了

leikunjiang

MeteoInfo 发表于 2014-11-21 14:44
已经转成了等经纬度数据了

王老师，FY2E轨道高度从HDF数据里面可以看到是3.57859e+007。就是那个h应该弄成这个数据。转换后的等经纬度数据的分辨率是0.05度吗？起始经纬度又该如何确定呢？

leikunjiang

MeteoInfo 发表于 2014-11-21 14:44
已经转成了等经纬度数据了

王老师，我回想起昨天检验NOMChannelVIS.grd里面的数据的时候，发现里面数据都是正整数，和HDF里面的数据是一样的。这说明这个文件里面的数据就不是等经纬度数据了撒，只是把HDF里的数据读出来保存在NOMChannelVIS.grd文件里了，它的数据部分也是2288行2288列。那经过转换成的等经纬度的数据到底是那个呢？

MeteoInfo

leikunjiang 发表于 2014-11-21 15:14
王老师，FY2E轨道高度从HDF数据里面可以看到是3.57859e+007。就是那个h应该弄成这个数据。转换后的等经纬 ...

输出的数据文件是Surfer ASCII Grid data数据格式，头几行有坐标参数，具体请自己网上搜索相关格式说明。

MeteoInfo

leikunjiang 发表于 2014-11-21 15:43
王老师，我回想起昨天检验NOMChannelVIS.grd里面的数据的时候，发现里面数据都是正整数，和HDF里面的数据 ...

我再强调一遍，脚本里下面的语句功能就是将Geostationary投影数据转换为等经纬度投影数据，如果脚本里gData = gData.project(toProj)语句未被注释掉（注释指的语句前面有#，该语言不会运行），那么输出的数据就是等经纬度数据。投影只改变坐标，并不改变数据本身。

        print 'Reproject grid data to Lon/Lat...'
        toProj = KnownCoordinateSystems.geographic.world.WGS1984
        gData = gData.project(toProj)