关于等值线叠加地图，裁剪后地图区域外遗留等值线标签和重新设置标签位置的方法

heidicat

一直以来，从论坛收获了很多，非常感谢大神们不吝惜的分享！！！！！
潜水多年，感觉是时候回报下社会了！因此，想把自己最近解决的一个问题跟大家分享一下。
问题来源是： 画等值线叠加地图时，对地图进行裁剪后发现，保留的地图区域外的等值线确实是裁掉了，但是保留地图区域外的等值线标签还在。
论坛有个类似的贴子也提出了这个问题：
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

topmad

如果底图是确定的那几个  可以预设一个字典 把label的location都定义好..

heidicat

topmad 发表于 2020-12-22 13:30
如果底图是确定的那几个  可以预设一个字典 把label的location都定义好..

那也要有等值线从那里经过吧

来自苹果APP客户端

little卓忘密码

{:5_235:}

edwardli

http://bbs.06climate.com/forum.p ... 2437&extra=page%3D2Python完美白化  专门增加了clabel的标识符，用以对clabel的返回值进行白化