python完美白化，去除区域外标签,去除区域外风矢量

晋陵小生

@平流层的萝卜  给出的python完美白化脚本非常好用（详见：[经验总结] Python完美白化）
可以很好地去除目标区域以外的等值线和填色。但是区域外的线标签无法去除，如：

这里在平流层的萝卜的maskout程序中多添加了一个可选参数
def shp2clip(originfig,ax,m,shpfile,region, clabel = None)


加入对应的处理语句：

1. if  clabel:
   
2.     clip_map_shapely = ShapelyPolygon(vertices)
   
3.     for text_object in clabel:
   
4.         if not clip_map_shapely.contains(ShapelyPoint(text_object.get_position())):
   
5.         text_object.set_visible(False)  

复制代码

对应绘图语句：

1. ###绘制填色图和色标
   
2. cf = m.contourf(x,y,pcp_grid, levels=levels, cmap=cmaps.CBR_wet)
   
3. cbar = m.colorbar(cf,location='right',format='%d',size=0.3,
   
4.         ticks=np.linspace(200,1100,19),label='毫米')
   
5. 
   
6. ###绘制等值线图和标签
   
7. contours = m.contour(x,y,pcp_grid, levels=levels,colors="black")
   
8. cn_label = plt.clabel(contours,inline=1,fontsize=10,fmt="%i")
   
9. 
   
10. #####擦除区域外的填色和等值线
    
11. clip = maskout.shp2clip(cf,ax,m,china_border,["吉林省","黑龙江省","辽宁省","内蒙古自治区"])
    
12. clip = maskout.shp2clip(contours,ax,m,china_border,["吉林省","黑龙江省","辽宁省","内蒙古自治区"],clabel=cn_label)

复制代码

最终效果图：

程序中的判断标识我也做了修改，大家可以根据自己用的shp文件进行调整（if shape_rec.record[7] in region:）



##########分隔符：2019/3/17更新###########


感谢@janemars提出白化风矢量的需求。风场返回的是一个元组，并且元素中包含有“set_clip_path”方法，与填色图的白化方法类似。
在shp2clip函数中增加一个可选参数，判断是否为矢量图

1. def shp2clip(originfig,ax,m,shpfile,region, clabel = None, vcplot = False):

复制代码

修改contour白化部分的脚本

1. if vcplot:
   
2.    for ivec in originfig:
   
3.         ivec.set_clip_path(clip)
   
4. else:
   
5.     for contour in originfig.collections:
   
6.         contour.set_clip_path(clip)

复制代码

对应绘图脚本（修改自wrf-python官方文档例子 Horizontal Interpolation to a Pressure Level）

1. from netCDF4 import Dataset
   
2. import numpy as np
   
3. import matplotlib.pyplot as plt
   
4. from matplotlib.cm import get_cmap
   
5. import maskout,sys
   
6. from wrf import getvar, interplevel, to_np, get_basemap, latlon_coords
   
7. 
   
8. # Open the NetCDF file
   
9. ncfile = Dataset("wrfout_d01_2017-01-26_06_00_00")
   
10. 
    
11. # Extract the pressure, geopotential height, and wind variables
    
12. p = getvar(ncfile, "pressure")
    
13. z = getvar(ncfile, "z", units="dm")
    
14. ua = getvar(ncfile, "ua", units="kt")
    
15. va = getvar(ncfile, "va", units="kt")
    
16. wspd = getvar(ncfile, "wspd_wdir", units="kts")[0,:]
    
17. 
    
18. # Interpolate geopotential height, u, and v winds to 500 hPa
    
19. ht_500 = interplevel(z, p, 500)
    
20. u_500 = interplevel(ua, p, 500)
    
21. v_500 = interplevel(va, p, 500)
    
22. wspd_500 = interplevel(wspd, p, 500)
    
23. 
    
24. # Get the lat/lon coordinates
    
25. lats, lons = latlon_coords(ht_500)
    
26. 
    
27. # Get the basemap object
    
28. bm = get_basemap(ht_500)
    
29. 
    
30. # Create the figure
    
31. fig = plt.figure(figsize=(12,9))
    
32. ax = plt.axes()
    
33. 
    
34. # Convert the lat/lon coordinates to x/y coordinates in the projection space
    
35. x, y = bm(to_np(lons), to_np(lats))
    
36. 
    
37. # Add the 500 hPa wind barbs, only plotting every 125th data point.
    
38. fbar = bm.barbs(x[::100,::100], y[::100,::100], to_np(u_500[::100, ::100]),
    
39.              to_np(v_500[::100, ::100]), length=6)
    
40. 
    
41. china_border = r"xxx/Export_Output_2"
    
42. bm.readshapefile(china_border, name = "china"  ,drawbounds=True, color = 'k',linewidth=1.)
    
43. clip = maskout.shp2clip(fbar,ax,bm,china_border,["吉林省","黑龙江省","辽宁省","内蒙古自治区"],vcplot=True)
    
44. plt.title("500 MB Height (dm), Wind Speed (kt), Barbs (kt)")
    
45. 
    
46. plt.show()

复制代码

对比图：




######################

若风矢量绘图返回的不是元组（可以print(fbar)查看），则maskout中的对应代码可以改为

1. if vcplot:
   
2.     originfig.set_clip_path(clip)
   
3. else:
   
4.     for contour in originfig.collections:
   
5.         contour.set_clip_path(clip)

复制代码

刚接触python不久，不清楚为什么会出现不同的返回值，大家有知道的可以帮忙补充一下，谢谢。

##########################################
############2019/4/1更新################################################################

判断风矢量是否可迭代

1. if vcplot:
   
2.     if isinstance(originfig,Iterable):
   
3.         for ivec in originfig:
   
4.             ivec.set_clip_path(clip)
   
5.     else:
   
6.         originfig.set_clip_path(clip)
   
7. else:
   
8.     for contour in originfig.collections:
   
9.         contour.set_clip_path(clip)

复制代码

白化语句：

1. clip = maskout.shp2clip(fbar,ax,bm,china_border,["吉林省","黑龙江省","辽宁省","内蒙古自治区"], vcplot=True) #白化风矢量场
   
2. clip = maskout.shp2clip(contour1,ax,bm,china_border,["吉林省","黑龙江省","辽宁省","内蒙古自治区"]) #白化填色图（无标签）
   
3. clip = maskout.shp2clip(contour2,ax,bm,china_border,["吉林省","黑龙江省","辽宁省","内蒙古自治区"]，clabel=True) #白化填色图（有标签）

复制代码

最新白化脚本：

maskout.py (3.17 KB, 下载次数: 951)

2019-4-1 15:47 上传

点击文件名下载附件


注意：maskout脚本中的第27行 “if shape_rec.record[7] in region:”需要根据自己使用的shp文件进行调整，区域的选定（["吉林省","黑龙江省","辽宁省","内蒙古自治区"]）也需要进行相应的调整，具体方法移步http://bbs.06climate.com/forum.p ... A%C3%C0%B0%D7%BB%AF

chiqu296

楼主能否分享一下你用的shp文件？

tulalang

这就非常好了，把等值线的label给去掉就很舒服了

晋陵小生

chiqu296 发表于 2019-3-14 15:33
楼主能否分享一下你用的shp文件？

可以啊，你把邮箱给我，我发给你

勿忘草

楼主可以给我分享一下你用的shp文件不？邮箱2452912925@qq.com。谢谢啦

janemars

楼主，请问这个报错是什么地方的问题呢？AttributeError: 'str' object has no attribute ' get_position'

晋陵小生

janemars 发表于 2019-3-15 17:31
楼主，请问这个报错是什么地方的问题呢？AttributeError: 'str' object has no attribute ' get_position'
...

可以看下你的脚本吗？
我怀疑可能是画标签的那边出问题了

1. ###绘制等值线图和标签
   
2. contours = m.contour(x,y,pcp_grid, levels=levels,colors="black")
   
3. cn_label = plt.clabel(contours,inline=1,fontsize=10,fmt="%i")

复制代码

janemars

晋陵小生 发表于 2019-3-16 16:14
可以看下你的脚本吗？
我怀疑可能是画标签的那边出问题了

谢谢楼主，我解决了

janemars

晋陵小生 发表于 2019-3-16 16:14
可以看下你的脚本吗？
我怀疑可能是画标签的那边出问题了

楼主，我画风杆遇到掩膜区域风杆去不掉的问题，不知道你有心得吗C:\Users\nriet\Documents\Tencent Files\949736402\FileRecv\MobileFile\IMG_7331

晋陵小生

janemars 发表于 2019-3-16 17:08
楼主，我画风杆遇到掩膜区域风杆去不掉的问题，不知道你有心得吗

已更新，请下载最新附件