请选择 进入手机版 | 继续访问电脑版
爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 132782|回复: 103

[经验总结] python完美白化,去除区域外标签,去除区域外风矢量

  [复制链接]

新浪微博达人勋

发表于 2019-3-14 10:48:41 | 显示全部楼层 |阅读模式

登录后查看更多精彩内容~

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

x
本帖最后由 晋陵小生 于 2019-4-3 11:02 编辑

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


这里在平流层的萝卜的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. contours = m.contour(x,y,pcp_grid, levels=levels,colors="black")
  7. cn_label = plt.clabel(contours,inline=1,fontsize=10,fmt="%i")

  8. #####擦除区域外的填色和等值线
  9. clip = maskout.shp2clip(cf,ax,m,china_border,["吉林省","黑龙江省","辽宁省","内蒙古自治区"])
  10. clip = maskout.shp2clip(contours,ax,m,china_border,["吉林省","黑龙江省","辽宁省","内蒙古自治区"],clabel=cn_label)
复制代码


最终效果图:

白化标签.png


程序中的判断标识我也做了修改,大家可以根据自己用的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. # Open the NetCDF file
  8. ncfile = Dataset("wrfout_d01_2017-01-26_06_00_00")

  9. # Extract the pressure, geopotential height, and wind variables
  10. p = getvar(ncfile, "pressure")
  11. z = getvar(ncfile, "z", units="dm")
  12. ua = getvar(ncfile, "ua", units="kt")
  13. va = getvar(ncfile, "va", units="kt")
  14. wspd = getvar(ncfile, "wspd_wdir", units="kts")[0,:]

  15. # Interpolate geopotential height, u, and v winds to 500 hPa
  16. ht_500 = interplevel(z, p, 500)
  17. u_500 = interplevel(ua, p, 500)
  18. v_500 = interplevel(va, p, 500)
  19. wspd_500 = interplevel(wspd, p, 500)

  20. # Get the lat/lon coordinates
  21. lats, lons = latlon_coords(ht_500)

  22. # Get the basemap object
  23. bm = get_basemap(ht_500)

  24. # Create the figure
  25. fig = plt.figure(figsize=(12,9))
  26. ax = plt.axes()

  27. # Convert the lat/lon coordinates to x/y coordinates in the projection space
  28. x, y = bm(to_np(lons), to_np(lats))

  29. # Add the 500 hPa wind barbs, only plotting every 125th data point.
  30. fbar = bm.barbs(x[::100,::100], y[::100,::100], to_np(u_500[::100, ::100]),
  31.              to_np(v_500[::100, ::100]), length=6)

  32. china_border = r"xxx/Export_Output_2"
  33. bm.readshapefile(china_border, name = "china"  ,drawbounds=True, color = 'k',linewidth=1.)
  34. clip = maskout.shp2clip(fbar,ax,bm,china_border,["吉林省","黑龙江省","辽宁省","内蒙古自治区"],vcplot=True)
  35. plt.title("500 MB Height (dm), Wind Speed (kt), Barbs (kt)")

  36. plt.show()
复制代码



对比图:
风场-未白化.png 风场-白化.png



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

若风矢量绘图返回的不是元组(可以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, 下载次数: 911)

评分

参与人数 3金钱 +30 贡献 +2 收起 理由
Om_MyGad + 10 赞一个!
chongzika + 10
mofangbao + 10 + 2

查看全部评分

密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2019-3-14 15:33:28 | 显示全部楼层
楼主能否分享一下你用的shp文件?
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2019-3-14 22:46:37 来自手机 | 显示全部楼层
这就非常好了,把等值线的label给去掉就很舒服了
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2019-3-15 07:28:52 | 显示全部楼层
chiqu296 发表于 2019-3-14 15:33
楼主能否分享一下你用的shp文件?

可以啊,你把邮箱给我,我发给你
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2019-3-15 09:05:53 | 显示全部楼层
楼主可以给我分享一下你用的shp文件不?邮箱2452912925@qq.com。谢谢啦
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2019-3-15 17:31:39 | 显示全部楼层
本帖最后由 janemars 于 2019-3-16 10:44 编辑

楼主,请问这个报错是什么地方的问题呢?AttributeError: 'str' object has no attribute ' get_position'
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2019-3-16 16:14:40 | 显示全部楼层
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")
复制代码
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2019-3-16 17:05:55 | 显示全部楼层
晋陵小生 发表于 2019-3-16 16:14
可以看下你的脚本吗?
我怀疑可能是画标签的那边出问题了

谢谢楼主,我解决了
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2019-3-16 17:08:24 | 显示全部楼层
本帖最后由 janemars 于 2019-3-16 17:12 编辑
晋陵小生 发表于 2019-3-16 16:14
可以看下你的脚本吗?
我怀疑可能是画标签的那边出问题了

楼主,我画风杆遇到掩膜区域风杆去不掉的问题,不知道你有心得吗C:\Users\nriet\Documents\Tencent Files\949736402\FileRecv\MobileFile\IMG_7331
QQ截图20190316171049.png
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2019-3-17 10:16:43 | 显示全部楼层
janemars 发表于 2019-3-16 17:08
楼主,我画风杆遇到掩膜区域风杆去不掉的问题,不知道你有心得吗

已更新,请下载最新附件
密码修改失败请联系微信:mofangbao
您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表