对站点资料做相关分析后画色斑图的问题

veir

按照站点资料画色斑图的思路，我对站点资料进行相关分析，分析后的结果画色斑图，结果色板图是下面这样。
相关分析的结果应该是没有问题，就是画图不知道怎么回事，但是我用java的MeteoInfo画出来了，不知道是不是插值算法问题，麻烦各位帮忙看看是画图哪里有问题吗？先谢谢各位了。
代码参考自
http://bbs.06climate.com/forum.php?mod=viewthread&tid=42437
以及http://bbs.06climate.com/forum.php?mod=viewthread&tid=57932&extra=page%3D1

附代码：

1. import os
   
2. import maskout2 as maskout
   
3. import numpy as np
   
4. import xarray as xr
   
5. from scipy.interpolate import Rbf
   
6. import matplotlib as mpl
   
7. import matplotlib.pyplot as plt
   
8. import cartopy.crs as ccrs
   
9. import cartopy.feature as cfeature
   
10. from cartopy.io.shapereader import Reader
    
11. import os
    
12. import time
    
13. import datetime
    
14. import pandas as pd
    
15. 
    
16. last_time = time.time()
    
17. 
    
18. plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文
    
19. plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号
    
20. 
    
21. SHP = './shapefile/china_shp'
    
22. root = os.path.dirname(__file__)
    
23. SHP = os.path.join(root, SHP)
    
24. 
    
25. # 数据处理
    
26. data = pd.read_csv('corr1.csv', header=0)
    
27. lon = data['Lon']
    
28. lat = data['Lat']
    
29. rain_data = data['corr']
    
30. 
    
31. 
    
32. # 插值
    
33. olon = np.linspace(71, 138, 88)
    
34. olat = np.linspace(16, 56, 88)
    
35. olon, olat = np.meshgrid(olon, olat)
    
36. func = Rbf(lon, lat, rain_data, function='cubic')
    
37. pre_grid = func(olon, olat)
    
38. 
    
39. # 字体
    
40. mpl.rcParams['font.sans-serif'] = ['Times New Roman']
    
41. mpl.rc('font', size=15, weight='normal')
    
42. font = {'family': 'SimHei', 'weight': 'normal', 'size': 15}
    
43. 
    
44. # 颜色
    
45. # 颜色
    
46. levels = [-1.0,
    
47.           -0.872,
    
48.           -0.764,
    
49.           -0.631,
    
50.           -0.549,
    
51.           0.0,
    
52.           0.549,
    
53.           0.631,
    
54.           0.764,
    
55.           0.872
    
56.           ]
    
57. crgb = [
    
58.     '#200080',
    
59.     '#1000bf',
    
60.     '#0000ff',
    
61.     '#007fff',
    
62.     '#00ffff',
    
63.     '#ffff00',
    
64.     '#ff7f00',
    
65.     '#ff0000',
    
66.     '#800020',
    
67. ]
    
68. cmaps = mpl.colors.ListedColormap(crgb)
    
69. norm = mpl.colors.BoundaryNorm(levels, 9)
    
70. 
    
71. # 画图
    
72. proj = ccrs.PlateCarree(central_longitude=105)
    
73. fig = plt.figure(figsize=(16, 9), facecolor='none')
    
74. ax = fig.add_subplot(1, 1, 1, projection=proj)
    
75. ax.set_extent([73, 136, 15, 55], ccrs.PlateCarree())
    
76. 
    
77. ax.add_geometries(Reader(os.path.join(SHP, 'cnhimap.shp')).geometries(),
    
78.                   ccrs.PlateCarree(),
    
79.                   facecolor='none',
    
80.                   edgecolor='#737373',
    
81.                   linewidth=1,
    
82.                   )
    
83. 
    
84. cf = ax.contourf(olon, olat, pre_grid, levels=levels, cmap=cmaps, norm=norm, zorder=0, transform=ccrs.PlateCarree())
    
85. 
    
86. # 位置[左,下,右,上]
    
87. pos_cb = []
    
88. cb_ax = fig.add_axes()
    
89. cbar = fig.colorbar(cf, ax=ax, orientation="vertical", aspect=25, pad=0.01, shrink=0.7)
    
90. cbar.set_label('', fontdict=font)
    
91. cbar.set_ticklabels(levels)
    
92. ax.set_title('站点资料温度-降水相关分析', fontdict=font)
    
93. 
    
94. #  南海
    
95. pos1 = ax.get_position()
    
96. # 位置[左,下,右,上]
    
97. pos2 = [pos1.x1 - ((pos1.x1 - pos1.x0) / 7), pos1.y0, pos1.width / 7, pos1.height / 5]
    
98. 
    
99. proj_n = ccrs.LambertConformal(central_latitude=90, central_longitude=115)
    
100. ax_n = fig.add_axes(pos2, projection=proj_n)
     
101. ax_n.set_extent([105, 125, 0, 25], ccrs.PlateCarree())
     
102. ax_n.add_geometries(Reader(os.path.join(SHP, 'cnhimap.shp')).geometries(),
     
103.                     ccrs.PlateCarree(),
     
104.                     facecolor='none',
     
105.                     edgecolor='k',
     
106.                     linewidth=1)
     
107. 
     
108. # 白化
     
109. clip = maskout.shp2clip(cf, ax, shpfile=os.path.join(SHP, 'country1.shp'), region='China', proj=proj)
     
110. 
     
111. # plt.show()
     
112. plt.savefig('out/corr1.jpg',
     
113.             format='jpg',
     
114.             # bbox_inches='tight',
     
115.             # transparent=True,
     
116.             dpi=600) # bbox_inches='tight' 图片边界空白紧致, 背景透明
     
117. 
     
118. print(f"cost: {datetime.timedelta(seconds=time.time() - last_time)} s")
     

复制代码

代码中使用到的相关分析结果是下面这样：


相关分析结果文件：
链接: https://pan.baidu.com/s/1GRLHfKrjuNPNzPIXJp-edA 提取码: 43ru

duoci01

请问python库中怎么安装maskout2?

dongyi

这个需要学习一波

来自安卓APP客户端

雪山笑侠

谢谢分享{:eb302:}{:eb302:}{:eb328:}