求助！站点降水插值为格点批量出图后，有正常有不正常的

nightelfxu

2022-07-23.xlsx (11.75 KB, 下载次数: 9)

2024-8-27 16:12 上传

点击文件名下载附件
不正常那天的数据并没有缺测 最近在做降水色斑图的批量出图，将站点数据通过调用Rbf插值函数插值为格点数据出图，结果有的降水落区就比较正常，有的非常奇怪，似乎插值有问题（检查过了数据里面没有缺测），有没有大神知道是怎么回事呢？附件Excel就是图片不正常那天的数据，并没有缺测，也找不到其他问题，求大神回复

1. # 引用部分
   
2. import numpy as np
   
3. import pandas as pd
   
4. import os
   
5. from datetime import datetime, timedelta
   
6. from scipy.interpolate import Rbf  # 径向基函数 ： 将站点信息插到格点上 用于绘制等值线
   
7. import matplotlib.pyplot as plt
   
8. import matplotlib.colors as colors
   
9. import matplotlib as mpl
   
10. from mpl_toolkits.basemap import Basemap
    
11. import matplotlib.patches as mpatches
    
12. import shapefile
    
13. #import maskout  #只显示某个地区
    
14. from matplotlib.path import Path
    
15. from matplotlib.patches import PathPatch
    
16. #m：坐标的投影对象；originfig：等值线对象；ax：绘制的句柄；shpfile：shp边界文件
    
17. def shp2clip(m,originfig,ax,shpfile):
    
18.     sf = shapefile.Reader(shpfile)
    
19.     clip=None
    
20.     vertices = []
    
21.     codes = []
    
22.     for shape_rec in sf.shapeRecords():
    
23.             pts = shape_rec.shape.points
    
24.             prt = list(shape_rec.shape.parts) + [len(pts)]
    
25.             # print(prt)
    
26.             for i in range(len(prt) - 1):
    
27.                 for j in range(prt[i], prt[i+1]):
    
28.                     vertices.append(m(pts[j][0], pts[j][1]))
    
29.                 codes += [Path.MOVETO]
    
30.                 codes += [Path.LINETO] * (prt[i+1] - prt[i] -2)
    
31.                 codes += [Path.CLOSEPOLY]
    
32.             clip = Path(vertices, codes)
    
33.             clip = PathPatch(clip, transform=ax.transData)
    
34.     for contour in originfig.collections:
    
35.         contour.set_clip_path(clip)
    
36.     return clip
    
37. plt.rcParams['font.sans-serif']=['SimHei'] #用来正常显示中文
    
38. plt.rcParams['axes.unicode_minus']=False #用来正常显示负号
    
39. # 数据准备
    
40. path = r"D:\\baoyu\\暴雨日\\2021-05-08.xlsx"
    
41. df  =  pd.read_excel(path,engine='openpyxl')    #  修改为读取Excel文件
    
42. name  =  df['id']
    
43. lons  =  df['lon']
    
44. lats  =  df['lat']
    
45. data  =  df['data0']
    
46. #绘制站点
    
47. sta =pd.read_excel("D:\\baoyu\\sta.xlsx", engine='openpyxl')
    
48. name = sta['Station_Name']
    
49. llon = sta['Lon']
    
50. llat = sta['Lat']
    
51. # 绘图准备
    
52. olon = np.linspace (115.0,126.0,550) #经纬坐标，有50个格点
    
53. olat = np.linspace (40.0,50.0,500) #纬度坐标，，有50个格点
    
54. olon,olat = np.meshgrid(olon,olat) #生成坐标网格 meshgrid网格化
    
55. func = Rbf(lons,lats,data,function='linear') #插值函数 调用Rbf插值函数中的插值法linear
    
56. rain_data_new = func(olon,olat) #插值
    
57. rain_data_new[rain_data_new <0 ] = 0
    
58. print(rain_data_new.shape)  # 打印形状
    
59. print(rain_data_new)  # 打印数据
    
60. #新建figure
    
61. fig = plt.figure(figsize=(19.6, 17.4))
    
62. #新建区域ax1
    
63. #figure的百分比,从figure 10%的位置开始绘制, 宽高是figure的80%
    
64. left, bottom, width, height = 0.01, 0.00, 0.98, 0.98
    
65. #获得绘制的句柄
    
66. ax1 = fig.add_axes([left, bottom, width, height])
    
67. m= Basemap(llcrnrlon=118.2, llcrnrlat=44.7, urcrnrlon=124.5, urcrnrlat=47.5, projection='lcc', lat_1=33, lat_2=45, lon_0=100)
    
68. m.readshapefile("D:/baoyu/xam/x",'xam', drawbounds=True)
    
69. #色彩定制：24小时降水量级色标
    
70. clevs = [0.0,0.1,10.,25.,50.,100.,250.,500] #自定义颜色列表
    
71. cdict = ['#FFFFFF','#A9F090','#40B73F','#63B7FF','#0000FE','#FF00FC','#850042'] #自定义颜色列表 '#A9F090','#40B73F','#63B7FF','#0000FE','#FF00FC','#850042'
    
72. my_cmap = colors.ListedColormap(cdict) # 自定义颜色映射 color-map
    
73. norm = mpl.colors.BoundaryNorm(clevs,my_cmap.N) # 基于离散区间生成颜色映射索引
    
74. x, y = m(olon, olat) # compute map proj coordinates.
    
75. xx,yy = m(llon,llat)
    
76. #  绘制等值线、等值线填色
    
77. cs = m.contourf(x,y,rain_data_new,clevs,cmap=my_cmap,norm=norm,extend='both')
    
78. st = m.scatter(xx-0.1,yy,c='k',s=10,marker='o') #绘制站点
    
79. for i in range(0,len(xx)):
    
80.     plt.text(xx[i],yy[i],sta['Station_Name'][i],va='center',fontsize=25)
    
81.     # 从数据路径中提取日期部分
    
82.     date_str = os.path.basename(path).split("\")[-1].split(".")[0]
    
83.     date_str1 = datetime.strptime(date_str, "%Y-%m-%d").strftime("%Y%m%d")  # 格式化日期字符串
    
84.     ymd = datetime.strptime(date_str, "%Y-%m-%d").strftime("%Y年%m月%d日")  # 格式化日期字符串
    
85.     previous_day = (datetime.strptime(date_str, "%Y-%m-%d") - timedelta(days=1)).strftime("%Y年%m月%d日")
    
86. #裁切
    
87. clip = shp2clip(m,cs,ax1,"D:/baoyu/xam/xamwbj")
    
88. xx,yy =m(121.8, 47.5)
    
89. plt.text(xx,yy,f'兴安盟降水实况图\n{previous_day}20时-{ymd}20时\n兴安盟气象台',fontsize=40,ha='center')
    
90. 
    
91. # 自定义图例标记
    
92. patch0 = mpatches.Patch(color=cdict[0], label='无降水')
    
93. patch1 = mpatches.Patch(color=cdict[1], label='0~10')
    
94. patch2 = mpatches.Patch(color=cdict[2], label='10~25')
    
95. patch3 = mpatches.Patch(color=cdict[3], label='25~50')
    
96. patch4 = mpatches.Patch(color=cdict[4], label='50~100')
    
97. patch5 = mpatches.Patch(color=cdict[5], label='100~250')
    
98. patch6 = mpatches.Patch(color=cdict[6], label='≥250')
    
99. # 构造图例
    
100. handles = [patch0, patch1, patch2, patch3, patch4, patch5, patch6]
     
101. ax1.legend(loc='lower left',handles=handles, title='图例(毫米)', title_fontsize=35, fontsize=30,
     
102.            shadow=False, frameon=True, edgecolor='#000000'  )
     
103. 
     
104. plt.savefig(f"D:\\baoyu\\暴雨日\\output_image_{date_str1}.png", dpi=300, bbox_inches='tight')
     
105. plt.show()
     
106. 
     

复制代码

付亚男

我做了一下测试
应该是插值的时候出现了问题
Rbf中要给出smooth，不能使用缺省的0





按道理讲，设为缺省的0也不应该出现问题
就好像你说的，有时候也能正常出图
不知道为什么会出现这种现象

nightelfxu

抱歉，不会设附件免费啊

付亚男

原因在这里

C0911和50833这两个站点的经纬度是相同的

nightelfxu

付亚男 发表于 2024-8-27 17:58
我做了一下测试
应该是插值的时候出现了问题
Rbf中要给出smooth，不能使用缺省的0

感谢，就是这个参数的问题，改了一下func = Rbf(lons, lats, data, function='linear', smooth=0.1) ，都能正常出图了。

nightelfxu

file:///C:/Users/Administrator/Desktop/20220723.png