python 画降水量色斑图问题

yangshuang

本人小白，通过网上学习的资料画降水量色斑图，并裁切成地图。但是最终的成果和我心中的理想值还是有落差，而且数据看起来不太对。
比如我有一个依兰、通河、方正三个区的降雨量分别为5mm、50mm、181mm，但图上渲染出来的无论是结果还是区域均对应不上。下边是我的代码,求解大神帮看下我哪里出了问题。其实我对整个绘制降雨量色斑图过程并不了解，都是网上查阅资料。

1. # 引用部分
   
2. import numpy as np
   
3. import pandas as pd
   
4. from scipy.interpolate import Rbf  # 径向基函数 ： 将站点信息插到格点上 用于绘制等值线
   
5. 
   
6. import matplotlib.pyplot as plt
   
7. import matplotlib.colors as colors
   
8. import matplotlib as mpl
   
9. 
   
10. import cartopy.crs as ccrs
    
11. import cartopy.io.shapereader as shpreader
    
12. from cartopy.mpl.ticker import LongitudeFormatter,LatitudeFormatter
    
13. 
    
14. import maskout  #只显示某个地区
    
15. 
    
16. plt.rcParams['font.sans-serif']=['SimHei'] #用来正常显示中文
    
17. plt.rcParams['axes.unicode_minus']=False #用来正常显示负号
    
18. 
    
19. # 数据准备
    
20. df = pd.read_excel('工作簿1.xlsx') #读取Excel
    
21. df.columns = ['stid','lon','lat','height','rain']
    
22. # df = df[(df['rain']>0)]  #筛选坐标在某个范围内，由于数据量少，直接在EXCEL中做处理  &(df['lat']>35)&(df['lon']>115)&(df['lat']<50)&(df['lon']<130)
    
23. lon = df['lon']
    
24. lat = df['lat']
    
25. rain = df['rain']
    
26. 
    
27. # 绘图准备
    
28. olon = np.linspace (125,131,120) #经纬坐标，0.05°分辨率 118°到126°0.05分辨率是有160个格点
    
29. olat = np.linspace (44,47,60) #纬度坐标，0.05°分辨率
    
30. olon,olat = np.meshgrid(olon,olat) #生成坐标网格 meshgrid网格化
    
31. func = Rbf(lon,lat,rain,function='linear') #插值函数 调用Rbf插值函数中的 cubic 插值法 linear
    
32. rain_data_new = func(olon,olat) #插值
    
33. rain_data_new[rain_data_new <0 ] = 0.
    
34. 
    
35. #画布及绘图声明
    
36. fig = plt.figure(figsize = (16,9.6),facecolor = '#666666',edgecolor = 'Blue',frameon = False)  # 画布
    
37. ax = fig.add_subplot(111,projection=ccrs.PlateCarree()) #绘图区
    
38. 
    
39. 
    
40. #色彩定制：24小时降水量级色标
    
41. clevs = [0.1,10.,25.,50.,100.,250.,500] #自定义颜色列表
    
42. cdict = ['#A9F090','#40B73F','#63B7FF','#0000FE','#FF00FC','#850042'] #自定义颜色列表 '#A9F090','#40B73F','#63B7FF','#0000FE','#FF00FC','#850042'
    
43. my_cmap = colors.ListedColormap(cdict) # 自定义颜色映射 color-map
    
44. norm = mpl.colors.BoundaryNorm(clevs,my_cmap.N) # 基于离散区间生成颜色映射索引
    
45. 
    
46. #  绘制等值线、等值线填色
    
47. cf = ax.contourf(olon,olat,rain_data_new,clevs,transform = ccrs.PlateCarree(),cmap=my_cmap,norm = norm)
    
48. ct = ax.contour(olon,olat,rain_data_new,clevs)   # 绘制等值线
    
49. clabel = ax.clabel(ct,fmt = '%i')
    
50. position = fig.add_axes([0.82,0.2,0.05,0.2]) #位置【左，下，宽。高】
    
51. plt.colorbar(cf,cax=position)     # 颜色参照表
    
52. position.set_yticklabels((0,10,25,50,100,250,500))
    
53. ax.xaxis.set_major_formatter(LongitudeFormatter(zero_direction_label=True))
    
54. ax.yaxis.set_major_formatter(LatitudeFormatter())
    
55. ax.set_xticks(np.arange(125,131,2),crs = ccrs.PlateCarree())    # x轴
    
56. ax.set_yticks(np.arange(44,47,2),crs = ccrs.PlateCarree())      # y轴
    
57. ax.gridlines()  #显示背景线
    
58. 
    
59. #裁切
    
60. # clip = maskout.shp2clip(cf,ax,'haerbin/hrb.shp')
    
61. clip = maskout.shp2clip(cf,ax,'haerbin/hrb.shp','hrb') #haerbin/shijie.shp  shp/黑龙江省/市界.shp
    
62. plt.show()
    
63. #从全国实际地图中获取辽宁省的市级边界并加载
    
64. # shpname = r'haerbin/hrb.shp'  # shp/country1.shp
    
65. # adm1_shapes=list(shpreader.Reader(shpname).geometries())
    
66. # ax.add_geometries(adm1_shapes[:],ccrs.PlateCarree(),edgecolor='k',facecolor='') #36:72东三省
    
67. # plt.savefig('111.png')

复制代码

笺疏

请问maskout是怎么定义的

yangshuang

笺疏 发表于 2020-3-11 15:20
请问maskout是怎么定义的

1. #coding=utf-8
   
2. ###################################################################################################################################
   
3. #####This module enables you to maskout the unneccessary data outside the interest region on a matplotlib-plotted output instance
   
4. ####################in an effecient way,You can use this script for free     ########################################################
   
5. #####################################################################################################################################
   
6. #####USAGE: INPUT  include           'originfig':the matplotlib instance##
   
7. #                                    'ax': the Axes instance
   
8. #                                    'shapefile': the shape file used for generating a basemap A
   
9. #                                    'region':the name of a region of on the basemap A,outside the region the data is to be maskout
   
10. #           OUTPUT    is             'clip' :the the masked-out or clipped matplotlib instance.
    
11. import shapefile
    
12. from matplotlib.path import Path
    
13. from matplotlib.patches import PathPatch
    
14. def shp2clip(originfig,ax,shpfile,region):
    
15.    
    
16.     sf = shapefile.Reader(shpfile,encoding='utf-8')
    
17.     for shape_rec in sf.shapeRecords():
    
18.         if shape_rec.record[2] == region:  ####这里需要找到和region匹配的唯一标识符，record[]中必有一项是对应的。
    
19.             vertices = []
    
20.             codes = []
    
21.             pts = shape_rec.shape.points
    
22.             prt = list(shape_rec.shape.parts) + [len(pts)]
    
23.             for i in range(len(prt) - 1):
    
24.                 for j in range(prt[i], prt[i+1]):
    
25.                     vertices.append((pts[j][0], pts[j][1]))
    
26.                 codes += [Path.MOVETO]
    
27.                 codes += [Path.LINETO] * (prt[i+1] - prt[i] -2)
    
28.                 codes += [Path.CLOSEPOLY]
    
29.             clip = Path(vertices, codes)
    
30.             clip = PathPatch(clip, transform=ax.transData)
    
31.     for contour in originfig.collections:
    
32.         contour.set_clip_path(clip)
    
33.     return clip
    

复制代码

Eegle

maskout中有误，增加对clabel的裁切，具体的可以在家园里搜索“白化”找到“晋陵小生”的那篇帖子，里面有讲解

hxylry

图画得很好啊，我还是小白，向你学习。

haihai123

我来学习一下。

Rainch

太漂亮了，学习一下，希望多多指教啊

snoopy-_-99

刚学python，出现问题都不知咋解决咋办