Python 空间绘图 - 等值线绘制
Python 空间绘图 - 等值线绘制
DataCharm
发布于 2021-02-22 12:03:13
发布于 2021-02-22 12:03:13
等值线是气象上比较常用的一种图形，特别是分析天气形势时，常用的地面气压、位势高度、气温等以等值线展示效果最好；在某些时候，我们还需要对等值线填色图进行进一步的美化。兹分别介绍之。
一、等值线基础的设定
从matplotlib的底层中，我们可以知道，等值线是基于绘图功能中的线条属性的，所以对于等值线来说，plt.plot命令的很多参数可以直接使用。不过需要注意的是——等值线每一根线条的值是不一样的，所以控制其颜色的参数为colors而不是color，当然，你如果在等值线里使用color='r'，程序并不会报错，但是颜色也不会改变；还有linewidth、linestyle也同理替换为linewidths、linestyles（所以lw、ls缩写也不能在contour中使用）。
二、等值线标签的问题与解决方法
与等值线填色图不一样，等值线的标签需要另外的clabel命令以绘制出来。其简要步骤如下：
ac=ax.contour(...)#这一步绘制等值线，并名为 ac
ax.clabel(cc)#这一步说明是在cc上绘制等值线标签然后问题来了，很多朋友在用再分析资料绘制时会发现该区域没有等值线标签，就几根光秃秃的线，不清楚其值究竟为多大。
在鸽了半个月之后，终于来解决这个问题。我推测这种情况出现的原因是在extent时，将标签的一部分或全部截去了，导致目标区域的等值线标签不全。目前有两种方法解决标签的放置问题：第一种，在绘图时仅取要绘制的部分，就仅取东经90-130，北纬20-50这部分的数据，这样标签应该是能够充分显示在本区域；第二种，使用clabel命令中的manual参数，强制每个标签的摆放位置，使全部标签都显示出来。可通过下面这个小程序理解：
import numpy as np
import matplotlib.pyplot as plt
import xarray as xr
import cartopy.crs as ccrs
import cartopy.feature as cf
import cartopy.io.shapereader as shpreader
from cartopy.mpl.ticker import LongitudeFormatter,LatitudeFormatter
import matplotlib.ticker as mticker
plt.rcParams['font.sans-serif']=['SimHei']
#########以下为数据准备工作################
proj=ccrs.PlateCarree()#简写投影
mapdata=shpreader.Reader(r'E:\hubei\hubei.shp')#读取shp文件
f=xr.open_dataset(r'E:\aaaa\datanc\he\fnl_20200715_12_00.nc')#读取再分析资料
lon=f['lon_0'][:]#读取经度
lat=f['lat_0'][:]#读取纬度
lons,lats=np.meshgrid(lon,lat)#网格化
HGT_P0_L100_GLL0=f['HGT_P0_L100_GLL0'][21][:][:]*0.1#我们常用的是位势十米，所以要乘0.1
RH_P0_L100_GLL0=f['RH_P0_L100_GLL0'][18][:][:]#相对湿度
UGRD_P0_L100_GLL0=f['UGRD_P0_L100_GLL0'][18][:][:]#纬向风
VGRD_P0_L100_GLL0=f['VGRD_P0_L100_GLL0'][18][:][:]#经向风
TMP_P0_L100_GLL0=f['TMP_P0_L100_GLL0'][21][:][:]-273.15#原来是开氏度，转化为摄氏度
##########封装地图函数########################
def create_map(ax):
ax.add_geometries(mapdata.geometries(),crs=proj,lw=0.45,facecolor='none',edgecolor='r',zorder=5)
ax.add_feature(cf.COASTLINE.with_scale('50m'),lw=0.6)
ax.add_feature(cf.RIVERS.with_scale('50m'),lw=0.5)
ax.set_xticks([90,95,100,105,110,115,120,125,130])#需要显示的经度，一般可用np.arange
ax.set_yticks([20,25,30,35,40,45,50])#需要显示的纬度
ax.xaxis.set_major_formatter(LongitudeFormatter())#将横坐标转换为经度格式
ax.yaxis.set_major_formatter(LatitudeFormatter())#将纵坐标转换为纬度格式
ax.tick_params(axis='both',labelsize=3,direction='in',length=2.75,width=0.55,right=True,top=True)#修改刻度样式
ax.grid(linewidth=0.4, color='k', alpha=0.45, linestyle='--')#开启网格线
ax.set_extent([90,130,20,50],crs=proj)
return ax
###########准备画布,绘制子图############################
fig=plt.figure(figsize=(2,5),dpi=500)
ax1=fig.add_axes([0,0,1,0.3],projection=proj)#添加三个子图
ax2=fig.add_axes([0,0.33,1,0.3],projection=proj)
ax3=fig.add_axes([0,0.66,1,0.3],projection=proj)
create_map(ax1)#让每个子图有地图与经纬度
create_map(ax2)
create_map(ax3)
###########首先是将数据全部绘制出来，不做取舍###################
a=ax3.contour(lons[:,:],lats[:,:],HGT_P0_L100_GLL0[:,:],linewidths=0.5,colors='k',levels=np.arange(500,600,4))
b=ax3.contour(lons[:,:],lats[:,:],TMP_P0_L100_GLL0[:,:],linewidths=0.5,linestyles='--',colors='r',levels=np.arange(-30,30,4))
ax3.clabel(a,inline=True,fmt='%.f',fontsize=3.5)
ax3.clabel(b,inline=True,fmt='%.f',fontsize=3.5)
ax3.barbs(lons[::2,::2],lats[::2,::2],UGRD_P0_L100_GLL0[::2,::2],VGRD_P0_L100_GLL0[::2,::2],linewidth=0.25,
barb_increments={'half':2,'full':4,'flag':20},sizes=dict(spacing=0.1,width=0.1,emptybarb=0.02,height=0.35),length=3.5,zorder=5)
ax3.set_ylabel('不做取舍\n(标签因为extent截取掉了)',fontsize=5)
##########然后是仅仅截取要绘制的部分数据来绘图##################
c=ax2.contour(lons[40:81,90:140],lats[40:81,90:140],HGT_P0_L100_GLL0[40:81,90:140],linewidths=0.5,colors='k',levels=np.arange(500,600,4))
d=ax2.contour(lons[40:81,90:140],lats[40:81,90:140],TMP_P0_L100_GLL0[40:81,90:140],linewidths=0.5,linestyles='--',colors='r',levels=np.arange(-30,30,4))
ax2.clabel(c,inline=True,fmt='%.f',fontsize=3.5)
ax2.clabel(d,inline=True,fmt='%.f',fontsize=3.5)
ax2.barbs(lons[40:81:2,90:140:2],lats[40:81:2,90:140:2],UGRD_P0_L100_GLL0[40:81:2,90:140:2],VGRD_P0_L100_GLL0[40:81:2,90:140:2],linewidth=0.25,
barb_increments={'half':2,'full':4,'flag':20},sizes=dict(spacing=0.1,width=0.1,emptybarb=0.02,height=0.35),length=3.5,zorder=5)
ax2.set_ylabel('仅取中国部分地区数据\n(可以正常显示了)',fontsize=5)
############第三种方法强制指定标签位置##################################
m=ax1.contour(lons[40:81,90:140],lats[40:81,90:140],HGT_P0_L100_GLL0[40:81,90:140],linewidths=0.5,colors='k',levels=[572,576,580,584,588])
n=ax1.contour(lons[40:81,90:140],lats[40:81,90:140],TMP_P0_L100_GLL0[40:81,90:140],linewidths=0.5,linestyles='--',colors='r',levels=np.arange(-30,30,4))
ax1.clabel(m,inline=True,fmt='%.f',fontsize=3.5,manual=[(100,47),(100,45),(100,43),(100,37),(125,27)])
ax1.clabel(n,inline=True,fmt='%.f',fontsize=3.5)
ax1.barbs(lons[40:81:2,90:140:2],lats[40:81:2,90:140:2],UGRD_P0_L100_GLL0[40:81:2,90:140:2],VGRD_P0_L100_GLL0[40:81:2,90:140:2],linewidth=0.25,
barb_increments={'half':2,'full':4,'flag':20},sizes=dict(spacing=0.1,width=0.1,emptybarb=0.02,height=0.35),length=3.5,zorder=5)
ax1.set_ylabel('指定标签位置',fontsize=5)
plt.title('各种方法下的等值线绘制',fontsize=8)
plt.savefig('a',bbox_inches='tight')
plt.show()
在上面最后一幅子图中，我们使用了manual参数，传入了一个存储了坐标的列表，列表中的坐标与等值线值一一对应。可以对比图2与图3，虽然都显示了足够的等值线标签，但图2的标签比较分散，图3的左侧标签统一集中在100°E这根线上。不过最后一种方式最好坐标与等值线一一对应，否则会出现不可预知的错误。
另外，我还查出了常用的一些等值线参数命令，希望能帮助到大家：
levels（contour） | 划分绘制的等值线间隔 |
|---|---|
fontsize（clabel） | 等值线标签的字号大小 |
colors（contour） | 等值线的颜色 |
inline（clabel） | True时，等值线在标签位置会断开 |
inline_spacing（clabel） | 等值线断裂的长度 |
fmt（clabel） | 标签的格式，可以调小数点位数等等 |
manual（clabel） | 强制指定标签的位置 |
rightside_up（clabel） | True时，标签将始终旋转90° |
use_clabeltext（clabel） | 修改标签的文本样式 |
下面是一个事例：
a=ax.contour(lons[40:81,90:140],lats[40:81,90:140],HGT_P0_L100_GLL0[40:81,90:140],colors='k',levels=np.arange(500,600,4))
ml=ax.clabel(a,inline=True,fmt='%.f',use_clabeltext=True)
for m in ml:
m.set_bbox({'fc': 'w'})#给每根线加上框
还有一个，在做白化时，平流层的萝卜的程序是无法把外面的标签白化的，晋陵小生做了优化，语皆在气象家园中。
三、等值线填色图的阴影区操作
在前面某些章节提到了等值线填色图的一些操作，下面是一个关于等值线填色图阴影绘图的方法。阴影区的绘制主要依靠contourf中的hatches参数进行设置，在前面我们应该已经大概了解了一些hatches的用法。
import numpy as np
import matplotlib as mpl
import pandas as pd
import cartopy.crs as ccrs
from cartopy.mpl.ticker import LongitudeFormatter,LatitudeFormatter
from cartopy.io.shapereader import Reader
from scipy.interpolate import Rbf
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
import sys
sys.path.append(r"C:\Users\lenovo\Desktop")
import maskout5
plt.rcParams['font.sans-serif']=['SimHei']
extent=[108.3,109.35,29.7,30.7]
proj= ccrs.PlateCarree()
fig = plt.figure(figsize=(4, 4),dpi=500)
ax = fig.subplots(1, 1, subplot_kw={'projection': proj})
reader = Reader(r'E:\家园\区划-省界\县.shp')
geo=list(reader.geometries())
ax.add_geometries(geo[1581:1582],crs=proj,edgecolor='k',facecolor='none',lw=1)
ax.set_extent(extent, crs=proj)
##############################读取文件打包数据###########################################
filename=r'C:\Users\lenovo\Desktop\利川累计2.xlsx'
df=pd.read_excel(filename)
lon=df['经度']#经度
lat=df['纬度']#纬度
name=df['站名']
lonlat=zip(lon,lat)
mapname=dict(zip(name,lonlat))
for key,value in mapname.items():
ax.scatter(value[0] , value[1] , marker='.' , s=10 , color = "k" , zorder = 3)
ax.text(value[0]-0.04 , value[1]+0.01 , key , fontsize = 4 , color = "k")
ax.set_xticks([108.5,108.73,108.94,109.15,109.349])
ax.set_yticks(np.arange(extent[2], extent[3]+0.2, 0.2))
ax.xaxis.set_major_formatter(LongitudeFormatter())
ax.yaxis.set_major_formatter(LatitudeFormatter())
ax.tick_params(axis='both',labelsize=5,direction='in',right=True,top=True)
ax.grid(linewidth=0.6, color='k', alpha=0.45, linestyle='--')
#####################################################################################
olon=np.linspace(108,111,110)
olat=np.linspace(29,32,110)
olon,olat=np.meshgrid(olon,olat)
rain=df['rain']
func=Rbf(lon,lat,rain,function='linear')
rain_new=func(olon,olat)
cs= ax.contourf(olon,olat,rain_new,colors='none',extend='both',hatches=['ooo','xxxx','*****','---','*','///','++++'])
clip=maskout5.shp2clip(cs,ax,r'E:\家园\区划-省界\县.shp',422802)
plt.title('累计降雨量',size=12)
position=fig.add_axes([0.78,0.17,0.02,0.3])
cb=fig.colorbar(cs,cax=position,extend='both',shrink=0.3,pad=0.01)
cb.set_label('累计降水量 $mm$',fontdict={'size':6})
cb.ax.tick_params(which='major',direction='in',length=6,labelsize=7)
plt.savefig("利川市累计降雨量.png",dpi=500, bbox_inches='tight')
plt.show()
上面这幅图，我们是关闭了填色（colors=None），这样就只能显示阴影显示填色了，当然你也可以打开颜色：
需要注意的是，hatches列表中的阴影样式数量必须与levels划分出来的间隔一致，否则会出现无法预知的错误。还请注意，使用了extend命令，使色条有尖尖时，尖尖也算一个等级。如上图，60-360划为了五份，但是两个小尖尖也算，所以有七个等级，hatches里面应该传入7个阴影样式。
你也可以在列表里放入None来使阴影仅展示在某个数值上，比如下图仅有180-240显示了阴影 ：
cs= ax.contourf(olon,olat,rain_new,levels=np.arange(60,361,60),cmap='Blues',extend='both',hatches=[None,None,None,'+++',None,None,None])
四、实验数据
链接：https://pan.baidu.com/s/1ZK47zL2XJjKn0e2ubDYNuw
提取码：oci0
谢谢阅读
本文参与 腾讯云自媒体分享计划，分享自微信公众号。
原始发表：2020-11-07，如有侵权请联系 cloudcommunity@tencent.com 删除
评论
登录后参与评论
推荐阅读