MeteoInfoLab脚本示例：MODIS Sinusoidal投影HDF数据

MeteoInfo

837078493 发表于 2015-10-29 22:10
http://pan.baidu.com/s/1o6irQJ4

需要知道你的模式的格点坐标，比如每个格点的经纬度，或者在模式投影坐标系下的坐标信息。

其实就是要将Sinusoida投影的数据投影到模式的格点上。

837078493

MeteoInfo 发表于 2015-10-30 13:31
需要知道你的模式的格点坐标，比如每个格点的经纬度，或者在模式投影坐标系下的坐标信息。

其实就是要 ...

模式的左上角经纬度是103.675,36.308右下角是104.583,35.5731  分辨率是1km

MeteoInfo

837078493 发表于 2015-10-30 16:07
模式的左上角经纬度是103.675,36.308右下角是104.583,35.5731  分辨率是1km

根据你提供的格点设置，获取格点在Lambert投影下的坐标（米为单位），用下面的脚本：

1. fromproj = projinfo(proj='longlat')
   
2. toproj = projinfo(proj='lcc', lon_0=104.137, lat_0=35.946, lat_1=30.0, lat_2=60.0)
   
3. lon = 104.583
   
4. lat = 35.5731
   
5. x, y = project(lon, lat, fromproj=fromproj, toproj=toproj)

fromproj = projinfo(proj='longlat')<br /> toproj = projinfo(proj='lcc', lon_0=104.137, lat_0=35.946, lat_1=30.0, lat_2=60.0)<br /> lon = 104.583<br /> lat = 35.5731<br /> x, y = project(lon, lat, fromproj=fromproj, toproj=toproj)

然后用下面的脚本将Sinusoidal投影数据转为Lambert投影数据。

1. f = addfile('D:/Temp/Hdf/MOD09A1.A2008217.h26v05.006.2015177090630.hdf')
   
2. vname = 'sur_refl_b01'
   
3. data = f[vname][:,:]
   
4. #Set projection
   
5. toproj = projinfo(proj='lcc', lon_0=104.137, lat_0=35.946, lat_1=30.0, lat_2=60.0)
   
6. #Set x/y
   
7. x = arange(-40586.90607311319, 39626.10751807479, 1000)
   
8. y = arange(-40422.11861349553, 39427.31922978533, 1000)
   
9. #Project
   
10. data = data.project(x, y, toproj)
    
11. #Write to a binary file
    
12. outfn = 'D:/Temp/test/' + vname + '.bin'
    
13. binwrite(outfn, data)
    
14. #Plot
    
15. axesm(projinfo=toproj, griddx=0.2, griddy=0.2)
    
16. mlayer = shaperead('D:/Temp/map/country1.shp')
    
17. geoshow(mlayer)
    
18. levs = arange(0, 100, 5)
    
19. layer = imshowm(x, y, data, 20, proj=toproj)
    
20. title(vname)
    
21. colorbar(layer)
    
22. #axism()

f = addfile('D:/Temp/Hdf/MOD09A1.A2008217.h26v05.006.2015177090630.hdf')<br /> vname = 'sur_refl_b01'<br /> data = f[vname][:,:]<br /> #Set projection<br /> toproj = projinfo(proj='lcc', lon_0=104.137, lat_0=35.946, lat_1=30.0, lat_2=60.0)<br /> #Set x/y<br /> x = arange(-40586.90607311319, 39626.10751807479, 1000)<br /> y = arange(-40422.11861349553, 39427.31922978533, 1000)<br /> #Project<br /> data = data.project(x, y, toproj)<br /> #Write to a binary file<br /> outfn = 'D:/Temp/test/' + vname + '.bin'<br /> binwrite(outfn, data)<br /> #Plot<br /> axesm(projinfo=toproj, griddx=0.2, griddy=0.2)<br /> mlayer = shaperead('D:/Temp/map/country1.shp')<br /> geoshow(mlayer)<br /> levs = arange(0, 100, 5)<br /> layer = imshowm(x, y, data, 20, proj=toproj)<br /> title(vname)<br /> colorbar(layer)<br /> #axism()

yysunnyboy

不明觉厉呀！   

837078493

MeteoInfo 发表于 2015-11-3 00:06
根据你提供的格点设置，获取格点在Lambert投影下的坐标（米为单位），用下面的脚本：

谢谢您啊   我有个关于投影的问题想请教一下  对于同一个经纬度坐标 该点对应的数据 如果用不同的投影得到的该点的数据值是一样的吗

MeteoInfo

837078493 发表于 2015-11-3 10:27
谢谢您啊   我有个关于投影的问题想请教一下  对于同一个经纬度坐标 该点对应的数据 如果用不同的投影得 ...

理论上应该是一样的，实际上由于浮点计算误差可能会有一点小的差别。

837078493

MeteoInfo 发表于 2015-11-3 10:29
理论上应该是一样的，实际上由于浮点计算误差可能会有一点小的差别。

哦哦  好的  太感谢您啦

gisorange

同志你好，我想请问如何将1km的modis影像与30米的landsat影像相匹配？还有在matlab里面如何实现影像的拼接，谢谢了

novelbean

不错，学习学习！！！

孤蓝et

MeteoInfo 发表于 2015-11-3 00:06
根据你提供的格点设置，获取格点在Lambert投影下的坐标（米为单位），用下面的脚本：

王老师，打算将卫星云图的圆盘图（Geostationary Satellite View）转换为其他投影方式进行绘图，
看了您的这个帖子，用您的demo和数据，但是其中有些参数设置不太理解，我怎样配合自己的数据来实现投影转换？

1、基本明白project函数的使用了，我也会输入一个x和y，得到另外一个投影方式下的x和y的数值。也明白demo中x = arange(-40586.90607311319, 39626.10751807479, 1000)的数值是通过这个方法转换而来的。

1. fromproj = projinfo(proj='longlat')
   
2. toproj = projinfo(proj='lcc', lon_0=104.137, lat_0=35.946, lat_1=30.0, lat_2=60.0)
   
3. lon = 104.583
   
4. lat = 35.5731
   
5. x, y = project(lon, lat, fromproj=fromproj, toproj=toproj)

复制代码

2、假如我针对自己手头的FY4数据，里面的参数怎么设置？利用下面的代码，输出data后发现里面都是0.

1. #Add data file
   
2. fn = 'D:/MeteoInfo/MeteoInfo.py/FY/FY4A-_AGRI--_N_REGC_1047E_L1-_FDI-_MULT_NOM_20190123043000_20190123043416_4000M_V0001.HDF'
   
3. f = addfile(fn)
   
4. print(f)
   
5. 
   
6. #Get data
   
7. data = f['NOMChannel02'][::-1,:]
   
8. 
   
9. #Set projection
   
10. #fromproj = projinfo(proj='geos', lon_0=104.74425, h=35781793)
    
11. #toproj = projinfo(proj='longlat')
    
12. toproj = projinfo() #此处为default，则默认是longlat经纬网格投影
    
13. 
    
14. #Set x/y
    
15. #x = arange(-40586.90607311319, 39626.10751807479, 1000)
    
16. #y = arange(-40422.11861349553, 39427.31922978533, 1000)
    
17. x=linspace(-5496000.0,5496000.0,2748)
    
18. y=linspace(-5496000.0,5496000.0,2748)
    
19. 
    
20. #Project
    
21. data = data.project(x, y, toproj)#此处仅包含toproj，那这个函数怎么知道原先这个数据是geos投影呢？
    
22. 
    
23. #保存数据至文本文件
    
24. f1 = open("out.txt",'w')
    
25. for tmp in data:
    
26.     f1.write(str(tmp))
    
27.     f1.write("\r")
    
28. f1.close()
    
29. 
    

复制代码

另外，我搜索了meteothink.org文档，没有找到data.project(x, y, toproj)的说明啊，不知道这个函数怎么使用？麻烦您指导一下如何将geos投影转换成其他类型投影的绘图方式，万分感谢。