python用RBF方法插值画降水落区，出现不可预知问题，求教

hut171

本人试过Python scipy.interpolate中各个插值方法，最后发现RBF方法还算差强人意，但是仍然存在一点问题

首先读取的txt文件格式为每一行“经度，纬度，降水”
当行数在386行（包括386行）以内的时候，可以正常做出落区图
当超过386行后，就成花图了。

请教坛里大神帮忙看下什么原因，怎么能够调好，毕竟区域站超过386个还是很正常的

1. import numpy as np
   
2. import matplotlib.pyplot as plt
   
3. from scipy.interpolate import Rbf
   
4. from mpl_toolkits.basemap import Basemap
   
5. 
   
6. 
   
7. def leves_colors():  # 设置色标方法2
   
8.     levels = [0, 0.000000001, 1.5, 7, 15, 40, 50]
   
9.     colors = ['#FFFFFF', '#A6F28F', '#3DBA3D', '#61B8FF', '#0000E1', '#FA00FA', '#800040']
   
10.     return levels, colors
    
11. 
    
12. 
    
13. filename = r'E:\qls\2019022604.txt'
    
14. file = np.loadtxt(filename, delimiter=',', dtype='str')
    
15. # lon = file[:, 0].astype(np.float).reshape(-1, 1)
    
16. lon = file[:, 0].astype(np.float)
    
17. lat = file[:, 1].astype(np.float)
    
18. rain = file[:, 2].astype(np.float)
    
19. result = np.ma.masked_greater(rain, 999989)
    
20. 
    
21. olon = np.linspace(104, 112, 300)
    
22. olat = np.linspace(31, 40, 300)
    
23. olon, olat = np.meshgrid(olon, olat)
    
24. 
    
25. func = Rbf(lon, lat, result, function='linear')
    
26. rain_data_new = func(olon, olat)
    
27. 
    
28. fig = plt.figure(figsize=(8, 8))
    
29. ax = fig.add_subplot(111)
    
30. m = Basemap(llcrnrlon=105, llcrnrlat=31, urcrnrlon=112, urcrnrlat=40, projection='cyl')
    
31. m.readshapefile('./Map/Shaanxi/Shaanxi_province', 'Shaanxi', default_encoding='utf-8')
    
32. xx, yy = m(olon, olat)
    
33. 
    
34. levels, colors = leves_colors()
    
35. c = m.contourf(xx, yy, rain_data_new, levels=levels, colors=colors, extend='max')
    
36. m.colorbar(c)
    
37. 
    
38. plt.show()
    
39. plt.close()

复制代码

数据文件示例，调试时可以在后面复制这些数据108.88,34.08,3.9
108.87,34.17,3.5
108.95,34.14,3.2
108.1,34.09,2.1
108.07,34.13,2.1
108.32,34.07,2.1
108.04,34.14,2
108.73,34.17,2
108.19,34.05,1.9
108.02,33.61,1.8
108.41,33.99,1.8
108.17,34.08,1.7
108.89,34.2,1.7

附件传错了，不知道怎么删除，大家千万不要下载，浪费金币

werewolf

站点太密了，这个函数针对少量站点比较好。站点密会出现病态矩阵。你可以考虑用LinearNDInterpolator或者natural neighbor。再者就自己写个反距离权重插值，我们学习一下。

hut171

LinearNDInterpolator这个我知道，natural neighbor这个方法在哪可以找到呀？

反距离权重插值，我在网上搜了下，挺复杂的，能力有限啊

835479131

hut171 发表于 2019-4-12 15:45
LinearNDInterpolator这个我知道，natural neighbor这个方法在哪可以找到呀？

反距离权重插值，我在网上 ...

你要的natural_neighbor
https://unidata.github.io/MetPy/ ... e.interpolate.html#

小其其格

用MI绘制，了解一下：

1. http://bbs.06climate.com/forum.php?mod=viewthread&tid=53210&extra=page%3D1

复制代码

hut171

好的，多谢上面各位专家，明天有空学习下

astiny

格点化的话，建议最好可以用cressman插值(对应ncl中的obj_anal_ic_deprecated方法）.metpy似乎有相关函数。

werewolf

astiny 发表于 2019-4-12 21:20
格点化的话，建议最好可以用cressman插值(对应ncl中的obj_anal_ic_deprecated方法）.metpy似乎有相关函数。

cressman反距离权重不推荐使用了，ncl官网推荐ESMF。不过python的esmfpy不会。

hut171

metpy下载了，准备先试试看cressman效果

astiny

werewolf 发表于 2019-4-13 09:48
cressman反距离权重不推荐使用了，ncl官网推荐ESMF。不过python的esmfpy不会。

您好。请教一下，具体是哪个函数？我只找到了从格点到格点的函数，没找到从站点到格点的函数。