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

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 11836|回复: 4

[经验总结] 关于等值线叠加地图,裁剪后地图区域外遗留等值线标签和重新设置标签位置的方法

[复制链接]

新浪微博达人勋

发表于 2020-12-22 11:15:05 | 显示全部楼层 |阅读模式

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

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

x
本帖最后由 heidicat 于 2020-12-22 11:15 编辑

一直以来,从论坛收获了很多,非常感谢大神们不吝惜的分享!!!!!
潜水多年,感觉是时候回报下社会了!因此,想把自己最近解决的一个问题跟大家分享一下。
问题来源是: 画等值线叠加地图时,对地图进行裁剪后发现,保留的地图区域外的等值线确实是裁掉了,但是保留地图区域外的等值线标签还在。
论坛有个类似的贴子也提出了这个问题:
http://bbs.06climate.com/forum.php?mod=viewthread&tid=47954&extra=page%3D1

下面举个栗子说明下解决办法(大神们有其他更简洁的办法请在后面留言哦~~~)

# -*- coding: utf-8 -*-
import netCDF4 as nc
import numpy as np
from shapely.geometry import Polygon as ShapelyPolygon  #先引入shapely,后引入matplotlib,否则可能会报错!
from shapely.geometry import Point as ShapelyPoint
import shapefile
import matplotlib.pyplot as plt
import matplotlib as mpl
from mpl_toolkits.basemap import Basemap
from matplotlib.patches import Polygon
from matplotlib.path import Path
from matplotlib.patches import PathPatch
from scipy.interpolate import Rbf
mpl.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus']=False
plt.rc('font',size=20,weight='bold')
###读离散点数据
with open(r'files\micaps3_threshold.000',"r") as f:   #读取台风极端降水阈值文件
context=f.readlines()
lat=[float(line.split()[2]) for line in context[5:] ]
lon=[float(line.split()[1]) for line in context[5:] ]
rainfall_threshold=[ float(line.split()[4]) for line in context[5:] ]
###设置网格
olon = np.linspace(115,121,50)      
olat = np.linspace(23,29,50)
olon,olat = np.meshgrid(olon,olat)
###离散点插值到网格
func = Rbf(lon,lat,rainfall_threshold,function='linear')  
rainfall_threshold_new = func(olon, olat)
###画等值线
levels =np.linspace(0,240,13,endpoint=True)
fig = plt.figure(figsize=(10,10))
ax = fig.add_subplot(111)
m = Basemap(projection='cyl',llcrnrlat=23.5,llcrnrlon=115.5,urcrnrlat=28.5,urcrnrlon=121)
cn =m.contour(olon,olat,rainfall_threshold_new,levels=levels,colors='k')
cn_label=plt.clabel(cn,cn.levels,inline=True,fmt="%.0f",fontsize=10)
###裁剪地图
shp_info2 = m.readshapefile(r'files\CHN_adm2','states',drawbounds=False,linewidth = 0.4,zorder=10)
for info, shp in zip(m.states_info, m.states):
proid = info['NAME_1']
if proid == 'Fujian':
  poly = Polygon(shp,facecolor='None',edgecolor='black', lw=0.5)
  ax.add_patch(poly)
vertices = []
codes = []   
shp_info = m.readshapefile(r'files\CHN_adm1','states',drawbounds=False,linewidth = 0.4,zorder=10)
for info, shp in zip(m.states_info, m.states):
proid = info['NAME_1']      
if proid == 'Fujian':        
  for j in np.arange(0, len(shp)):
   vertices.append((shp[j][0], shp[j][1]))
  codes += [Path.MOVETO]
  codes += [Path.LINETO] * (len(shp)-2)
  codes += [Path.CLOSEPOLY]
  clip = Path(vertices, codes)
  clip = PathPatch(clip, transform=ax.transData)
for contour in cn.collections:
contour.set_clip_path(clip)
###经纬度标注
font={'size' : 15 }
yticks=range(24,29,1)
xticks=range(116,122,2)
ax.set_xticks(xticks)
ax.set_yticks(yticks)
xlabels=[r'$'+str(i)+'^\circ E$' for i in xticks]
ylabels=[r'$'+str(i)+'^\circ N$' for i in yticks]
ax.set_yticklabels(ylabels,fontdict=font)
ax.set_xticklabels(xlabels,fontdict=font)
plt.show()



见图1(with contour labels outside the map area),图上有些标注不在福建省区域内!
由于contour等值线一条只有一个label,有的label在福建省外,所以等值线裁剪掉以后,标签留下了。
在前面那个帖子里我回复的是直接把地图外的标签设为invisible,但是这样就导致有些等值线没有label了。
因此,这里给出的方法是 采集每条等值线的经纬度位置,把label的位置设置在等值线上位于福建省内的经纬度(好像有点绕~~~)
cn_label=plt.clabel(cn,cn.levels,inline=True,fmt="%.0f",fontsize=10) # 从这里开始添加程序块
###############添加开始#######################
sf = shapefile.Reader(r'files\CHN_adm1.shp')  #读取福建省的经纬度范围
recs = sf.records()
shapes = sf.shapes()
for i in range(len(recs)):
if recs[4]=='Fujian':
  fj_area=ShapelyPolygon(shapes.points)
  break
# 重新计算contour labels的位置
'''
获取每条等值线的经纬度信息,levels有13个等级,就有13个line,每个line
可以对于多条等值线,比如说 40mm的等值线有3根,则line=2,len(path.vertices)=3
'''
keep_labels=[]
for i,line in enumerate(cn.collections):
for path in line.get_paths():         
  for kk in range(len(path.vertices)):
   if fj_area.contains(ShapelyPoint(path.vertices[kk]))==True:
   #判断等值线上的点是否在福建省区域内
    keep_labels.append(path.vertices[kk])  #如果等值线的点在福建省内就保留作为重新绘制label的位置
    break           
for label in cn_label:  #删掉之前的contour labels
label.remove()
# 重新画labels
cn_label=plt.clabel(cn,cn.levels,inline=True,fmt="%.0f",fontsize=10,manual=keep_labels)
################添加结束#######################
###裁剪地图

~~~~~~~~~~~~~~~~~~~~~~~~~
见图2(without contour labels outside the map area),从图上可以看到,地图区域外的label都没了,但是有个问题就是很多等值线的label都在等值线与地图的交界处,那是因为上面的这段程序我只判断等值线上的点在福建省内就作为label的位置,所以基本上都是起始点。
如果需要标注在中间,可以把在等值线上在福建省内的点采集出来,再选择中间位置的点作为label的位置,也就是说你可以任意操作了!!!

第一次在论坛发帖,希望对大家有所帮助!
最后感谢StackOverflow网站的大神们,从他们的回复里得到很多启发,以及死胖纸topmad




with conotur labels outside the map area

with conotur labels outside the map area

without contour labels outside the map area

without contour labels outside the map area

评分

参与人数 1威望 +10 金钱 +200 贡献 +5 体力 +100 收起 理由
topmad + 10 + 200 + 5 + 100 很给力!

查看全部评分

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

新浪微博达人勋

发表于 2020-12-22 13:30:29 | 显示全部楼层
如果底图是确定的那几个  可以预设一个字典 把label的location都定义好..
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2020-12-22 14:36:41 来自手机 | 显示全部楼层
topmad 发表于 2020-12-22 13:30
如果底图是确定的那几个  可以预设一个字典 把label的location都定义好..

那也要有等值线从那里经过吧
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2020-12-22 19:09:44 | 显示全部楼层
{:5_235:}
密码修改失败请联系微信:mofangbao
回复

使用道具 举报

新浪微博达人勋

发表于 2020-12-23 11:09:45 | 显示全部楼层
http://bbs.06climate.com/forum.p ... 2437&extra=page%3D2Python完美白化  专门增加了clabel的标识符,用以对clabel的返回值进行白化
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

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

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

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