用 我 所 能,精 彩 气 象。
在实际业务中经常需要对指定经纬度点进行一个相关气象数据的分析和研究,需要将格点数据插值到站点上面。本文介绍了三种在MeteoInfoLab中如何将格点数据插值到站点上面的方法。【本文参考了王老师的书和代码】
格点数据插值到站点主要有两种方法:双线性插值和最近距离,算法都很简单,MeteoInfoLab中插值到站点有几种方法:
(a)利用DimDataFile的tostation方法
(b)利用DimArray的tostation方法
(c)利用interp2d插值函数。推荐使用interp2d方法,该方法中的kind参数缺省为'linear'双线性插值,也可以设置为kind='neareast'最近距离插值(其实就是找离站点最近的格点将其值赋给站点)
总结:其实这几种方法插值出来的结果都差不多,王老师也推荐使用interp2d。
f = addfile('model.ctl')
ps = f['PS'][0,'10:60','60:140']
#Interpolate to a point
x = 123
y = 44.5
z = None
t = 0
d = f.tostation('PS', x, y, z, t)
d1 = ps.tostation(x, y)
d2 = interp2d(ps, x, y)
d3 = interp2d(ps, x, y, kind='neareast')
print d
print d1
print d2
print d3
#942.347668457
#942.347668457
#942.347668457
#980.492126465
#Plot
figure(figsize=[600,600], newfig=False,bgcolor=[255,255,255])
axesm()
geoshow('cn_province.shp',labelfield='NAME', fontname=u'黑体', fontsize=12, yoffset=15,avoidcoll=False)
scatterm(x, y, d, size=10, color='g')
text(x + 0.3, y+1, 'DimDataFile.tostation:%.4f' % d, color='r', fontsize=12)
text(x + 0.3, y-1, 'DimArray.tostation:%.4f' % d, color='r', fontsize=12)
text(x + 0.3, y-2, 'interp2d.linear:%.4f' % d, color='r', fontsize=12)
text(x + 0.3, y-3, 'interp2d.neareast:%.4f' % d, color='r', fontsize=12)
xlim(110, 132)
ylim(32, 54)
title('Grid to point interpolation')