爱气象，爱气象家园！ 
·设为首页 开启辅助访问 新帖子新回复热门贴精华帖

气象家园

 找回密码
 立即注册

QQ登录

只需一步，快速开始

新浪微博登陆

只需一步, 快速开始

搜索
气象家园»论坛 › 精品推荐 › MeteoInfo › 处理FY2E圆盘标称投影数据
123456789下一页
返回列表 发新帖
楼主: MeteoInfo

处理FY2E圆盘标称投影数据

[复制链接]
MeteoInfo

新浪微博达人勋
 楼主| 发表于 2014-11-18 09:15:27 | 显示全部楼层
沙颖凯 发表于 2014-11-18 04:10
我想我可以做一些这方面的事情，比如基于楼主的源代码写Python的模块之类，问下楼主的MeteoInfo是不是开 ...

MeteoInfo目前还没有开源。
密码修改失败请联系微信：mofangbao
回复 支持 反对

使用道具 举报

leikunjiang

新浪微博达人勋
发表于 2014-11-18 10:32:49 | 显示全部楼层
MeteoInfo 发表于 2014-11-17 22:43
数据南北反向的问题已解决（见1楼帖子），在脚本中加一句：
gData.yReverse()

就是说，这个FY2E全圆盘数据有2288行2288列，光有这个数据是不能得到它在地球上的经纬度位置的。因为它有个专门的经纬度文件，纬度有2288行2288列，经度也是有2288行2288列数据，全圆盘里每个数据的位置对应着这个经纬度文件。
密码修改失败请联系微信：mofangbao
回复 支持 反对

使用道具 举报

leikunjiang

新浪微博达人勋
发表于 2014-11-18 10:38:26 | 显示全部楼层
MeteoInfo 发表于 2014-11-17 22:43
数据南北反向的问题已解决（见1楼帖子），在脚本中加一句：
gData.yReverse()

那个经纬度文件我通过邮箱传给你了的
密码修改失败请联系微信：mofangbao
回复 支持 反对

使用道具 举报

MeteoInfo

新浪微博达人勋
 楼主| 发表于 2014-11-18 10:42:56 | 显示全部楼层
leikunjiang 发表于 2014-11-18 10:32
就是说，这个FY2E全圆盘数据有2288行2288列，光有这个数据是不能得到它在地球上的经纬度位置的。因为它有 ...

你到底要做什么？能不能说清楚些？不是要转为等经纬度数据吗？
密码修改失败请联系微信：mofangbao
回复 支持 反对

使用道具 举报

leikunjiang

新浪微博达人勋
发表于 2014-11-18 11:31:19 | 显示全部楼层
MeteoInfo 发表于 2014-11-18 10:42
你到底要做什么？能不能说清楚些？不是要转为等经纬度数据吗？

sorry，可能是表达的问题哈，就是说我觉得FY2E全圆盘数据要依靠那个专门的经纬度数据才能定位在地球上，你的第一张图（可见光通道），是不是MeteoInfo软件像micaps软件一样可以不依靠那个经纬度数据就能打开或者直接生成填色图？因为我安装的meteinfo包括室友安装的meteinfo都打不开FY2E全圆盘数据。我也问了其他同学，他们也不知道第一张图的底图3D模型是怎么弄出来的。
密码修改失败请联系微信：mofangbao
回复 支持 反对

使用道具 举报

MeteoInfo

新浪微博达人勋
 楼主| 发表于 2014-11-18 12:52:57 | 显示全部楼层
本帖最后由 MeteoInfo 于 2014-11-18 12:54 编辑
leikunjiang 发表于 2014-11-18 11:31
sorry，可能是表达的问题哈，就是说我觉得FY2E全圆盘数据要依靠那个专门的经纬度数据才能定位在地球上， ...

通常卫星数据无论是什么投影，数据在该投影中的坐标是等间隔的，也就是格点数据。只有在等经纬度投影中坐标的单位才是度，其它投影中的坐标单位通常是米。所以想要定位Geostationary投影下的数据，需要知道投影参数及该投影下数据坐标的起始点、x/y方向上的格距和格点数，也就是脚本程序中的xmin, ymin, xdelt, ydelt, xnum, ynum。遗憾的是FY2E的hdf数据中并没有给出这样的信息（xnum和ynum是有的2288），好在通过调试也能大概给出（脚本中的值）。你说的经纬度数据在这里用不上。micaps正确打开该数据应该是软件里给出了此类数据的固定参数。其实如果FY系列的hdf数据遵循EOS约定就好办了，hdf格式本身保罗万象，不遵循某一约定的话无法用软件自动读取定位。
首先下载MeteoInfo Java最新版本。要生成1楼帖子里第一张图，先将地图投影成数据的投影（“显示 -> 投影”菜单），具体参考第二张图上的参数。然后运行下面的脚本程序即可。这个脚本和1楼帖子里的基本一样，只是无需再进行投影转换。

  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

  19. #Set data directory
  20. dataDir = 'D:/Temp/hdf/'

  22. #Create MeteoDataInfo object
  23. mdi = MeteoDataInfo()

  25. projInfo = ProjectionInfo("+proj=geos +lon_0=104.5 +h=35785864")

  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)

  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.         print 'Ouput grid data to Sufer ASCII grid data file...'
  71.         outfn = 'D:/Temp/test/NOMChannelVIS.grd'
  72.         #gData.saveAsSurferASCIIFile(outfn)
  73.         print 'Create raster layer from the grid data...'
  74.         aLS = LegendManage.createLegendSchemeFromGridData(gData, LegendType.GraduatedColor, ShapeTypes.Polygon)
  75.         aLayer = DrawMeteoData.createRasterLayer(gData, "Test_FY2E", aLS)
  76.         #aLayer = DrawMeteoData.createShadedLayer(gData, aLS, 'Test_HDF_Shaded', 'AOD', True)
  77.         #aLayer.setProjInfo(toProj)
  78.         aLayer.setProjInfo(projInfo)
  79.         print 'Add the layer...'
  80.         mf = miapp.getMapDocument().getActiveMapFrame()
  81.         mf.addLayer(aLayer)
  82.         mf.moveLayer(aLayer, 0)

  84. print 'Finished!'



密码修改失败请联系微信：mofangbao
回复 支持 反对

使用道具 举报

小预报员

新浪微博达人勋
发表于 2014-11-18 17:22:46 | 显示全部楼层
真漂亮啊镇漂亮
密码修改失败请联系微信：mofangbao
回复 支持 反对

使用道具 举报

leikunjiang

新浪微博达人勋
发表于 2014-11-18 18:52:23 | 显示全部楼层
MeteoInfo 发表于 2014-11-18 12:52
通常卫星数据无论是什么投影，数据在该投影中的坐标是等间隔的，也就是格点数据。只有在等经纬度投影中坐 ...

让我受教良多啊，我用这个数据的最终目的是用它反演降水，所以对数据定位要求比较高，通常所使用的投影是等经纬度，所以思路绕不开这个。就FY2E全圆盘数据，文档里讲的各个通道的分辨率是5KM(0.05度)，我算是明白了
密码修改失败请联系微信：mofangbao
回复 支持 反对

使用道具 举报

MeteoInfo

新浪微博达人勋
 楼主| 发表于 2014-11-18 19:51:01 | 显示全部楼层
leikunjiang 发表于 2014-11-18 18:52
让我受教良多啊，我用这个数据的最终目的是用它反演降水，所以对数据定位要求比较高，通常所使用的投影是 ...

如果知道确切的数据分辨率（5000 m），就可以算出来确切的参数，详见下面的脚本：
  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

  19. #Set data directory
  20. dataDir = 'D:/Temp/hdf/'

  22. #Create MeteoDataInfo object
  23. mdi = MeteoDataInfo()

  25. projInfo = ProjectionInfo("+proj=geos +lon_0=104.5 +h=35785864")

  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 = -5717500.0
  36.         ymin = -5717500.0
  37.         xnum = 2288
  38.         ynum = 2288
  39.         xdelt = 5000
  40.         ydelt = 5000
  41.         xlist = []
  42.         ylist = []
  43.         for i in range(0,xnum):
  44.                 xlist.append(xmin + xdelt * i)
  45.         for i in range(0,ynum):
  46.                 ylist.append(ymin + ydelt * i)

  48.         X = jarray.array(xlist, 'd')
  49.         Y = jarray.array(ylist, 'd')
  50.         xDim = Dimension(DimensionType.X)
  51.         xDim.setValues(X)
  52.         dataInfo.setXDimension(xDim)
  53.         yDim = Dimension(DimensionType.Y)
  54.         yDim.setValues(Y)
  55.         dataInfo.setYDimension(yDim)
  56.         #var = dataInfo.getVariable('FY2E_CLM_Hourly_Cloud_Amount')
  57.         var = dataInfo.getVariable('NOMChannelVIS')
  58.         dimList = [yDim, xDim]
  59.         var.setDimensions(dimList)
  60.         print 'Get grid data...'
  61.         gData = mdi.getGridData(var.getName())       
  62.         gData.yReverse()
  63.         gData.projInfo = projInfo
  64.         #gData.missingValue = 0.0
  65.         print 'Reproject grid data to Lon/Lat...'
  66.         toProj = KnownCoordinateSystems.geographic.world.WGS1984
  67.         #gData = gData.project(toProj)
  68.         print 'Ouput grid data to Sufer ASCII grid data file...'
  69.         outfn = 'D:/Temp/test/NOMChannelVIS.grd'
  70.         #gData.saveAsSurferASCIIFile(outfn)
  71.         print 'Create raster layer from the grid data...'
  72.         aLS = LegendManage.createLegendSchemeFromGridData(gData, LegendType.GraduatedColor, ShapeTypes.Polygon)
  73.         aLayer = DrawMeteoData.createRasterLayer(gData, "Test_FY2E", aLS)
  74.         #aLayer = DrawMeteoData.createShadedLayer(gData, aLS, 'Test_HDF_Shaded', 'AOD', True)
  75.         #aLayer.setProjInfo(toProj)
  76.         aLayer.setProjInfo(projInfo)
  77.         print 'Add the layer...'
  78.         mf = miapp.getMapDocument().getActiveMapFrame()
  79.         mf.addLayer(aLayer)
  80.         mf.moveLayer(aLayer, 0)

  82. print 'Finished!'


密码修改失败请联系微信：mofangbao
回复 支持 反对

使用道具 举报

leikunjiang

新浪微博达人勋
发表于 2014-11-18 20:10:51 | 显示全部楼层
MeteoInfo 发表于 2014-11-18 19:51
如果知道确切的数据分辨率（5000 m），就可以算出来确切的参数，详见下面的脚本：

各个通道的分辨率确实是5km，不过2288乘以5000等于11440000，我不知道程序中的起始点-5717500是怎么得来的
密码修改失败请联系微信：mofangbao
回复 支持 反对

使用道具 举报

下一页 »
123456789下一页
返回列表 发新帖
高级模式
B Color Image Link Quote Code Smilies
您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

小黑屋|帮助中心|清除痕迹

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表，不代表气象家园立场，禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表