MeteoInfoLab脚本示例：格点插值到站点

MeteoInfo

格点数据插值到站点主要有两种方法：双线性插值和最近距离，算法都很简单，MeteoInfoLab中插值到站点有几种方法：（a）利用DimDataFile的tostation方法，（b）利用DimArray的tostation方法，（c）利用interp2d插值函数。推荐使用interp2d方法，该方法中的kind参数缺省为'linear'双线性插值，也可以设置为kind='neareast'最近距离插值（其实就是找离站点最近的格点将其值赋给站点）。

示例脚本：

1. f = addfile('D:/Temp/GrADS/model.ctl')
   
2. ps = f['PS'][0,'10:60','60:140']
   
3. 
   
4. #Interpolate to a point
   
5. x = 123
   
6. y = 44.5
   
7. z = None
   
8. t = 0
   
9. d = f.tostation('PS', x, y, z, t)
   
10. d1 = ps.tostation(x, y)
    
11. d2 = interp2d(ps, x, y)
    
12. d3 = interp2d(ps, x, y, kind='neareast')
    
13. print d
    
14. print d1
    
15. print d2
    
16. print d3
    
17. 
    
18. #Plot
    
19. axesm()
    
20. geoshow('country', edgecolor=(0,0,255))
    
21. layer = scatterm(ps, fill=False, edgecolor='r')
    
22. layer.addlabels('data', yoffset=15, decimals=2)
    
23. scatterm(x, y, d, size=10, color='g')
    
24. text(x + 0.3, y, '%.2f' % d)
    
25. xlim(x - 6, x + 6)
    
26. ylim(y - 5, y + 5)
    
27. title('Grid to point interpolation')

f = addfile('D:/Temp/GrADS/model.ctl')<br /> ps = f['PS'][0,'10:60','60:140']<br /> <br /> #Interpolate to a point<br /> x = 123<br /> y = 44.5<br /> z = None<br /> t = 0<br /> d = f.tostation('PS', x, y, z, t)<br /> d1 = ps.tostation(x, y)<br /> d2 = interp2d(ps, x, y)<br /> d3 = interp2d(ps, x, y, kind='neareast')<br /> print d<br /> print d1<br /> print d2<br /> print d3<br /> <br /> #Plot<br /> axesm()<br /> geoshow('country', edgecolor=(0,0,255))<br /> layer = scatterm(ps, fill=False, edgecolor='r')<br /> layer.addlabels('data', yoffset=15, decimals=2)<br /> scatterm(x, y, d, size=10, color='g')<br /> text(x + 0.3, y, '%.2f' % d)<br /> xlim(x - 6, x + 6)<br /> ylim(y - 5, y + 5)<br /> title('Grid to point interpolation')

又是那隻貓

weinihou

{:eb502:}{:eb502:}{:eb502:}{:eb502:}

头两天才见到有人提如何使用最邻近法，王老师怎么这么快就附上啦！！

lovechang1314

想知道这个是否考虑了投影的问题？
给的站点经纬度需要转换到ctl的投影下再插值，不知道函数tostation内部是怎么考虑的。
此外interp2d显然是不会考虑投影的，直接将站点经纬度当成了ctl里的经纬度

MeteoInfo

lovechang1314 发表于 2017-6-30 15:46
想知道这个是否考虑了投影的问题？
给的站点经纬度需要转换到ctl的投影下再插值，不知道函数tostation内部 ...

如果格点数据不是经纬度投影，需要将站点经纬度投影为和格点一样的投影坐标系，MeteoInfoLab有投影函数可以用：http://www.meteothinker.com/docs ... ctions/project.html

lovechang1314

MeteoInfo 发表于 2017-6-30 16:05
如果格点数据不是经纬度投影，需要将站点经纬度投影为和格点一样的投影坐标系，MeteoInfoLab有投影函数可 ...

谢谢王老师！
可以在示例中加以说明，我走过了这个坑，没注意，后来看不太对，想起来问一下

Honlenn

感谢分享，非常有用

孤蓝et

王老师，咨询下，假如我指定一个站点的经纬度，能直接获得他周围四个网格点的数据么？

MeteoInfo

孤蓝et 发表于 2018-11-29 09:45
王老师，咨询下，假如我指定一个站点的经纬度，能直接获得他周围四个网格点的数据么？

知道格点坐标当然就可以

孤蓝et

嗯嗯，理解了，大致思路之前用meteoinfo类库实现，就是先找到站点周围的四个网格点经纬度信息，然后循环找到这四个点，就是在读取三年的数据效率有点低，比较耗时。