Python 空间绘图 - 等值线绘制

ac=ax.contour(...)#这一步绘制等值线，并名为 ac
ax.clabel(cc)#这一步说明是在cc上绘制等值线标签

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()

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'})#给每根线加上框

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()

cs= ax.contourf(olon,olat,rain_new,levels=np.arange(60,361,60),cmap='Blues',extend='both',hatches=[None,None,None,'+++',None,None,None])