Python气象处理绘图（三）–青藏高原空间分布白化

mask和colorbar设定

Python气象处理绘图（三）--青藏高原空间分布白化
前言
一、maskout子程序思路？
二、绘制白化图形
总结

前言

前期进行了黄河流域降水分布的白化之后，最近尝试了对于青藏高原气象要素空间分布进行了白化。尝试绘制多子图，并学习了如何添加多子图的colorbar。

一、maskout子程序思路？

代码如下（参考微信公众号：气象水文科研猫发布的mask内容）：
在进行特征量选取的时候，要选择有效的shape_rec,可以提前用arcmap等绘图软件检查shapeRecords的相关信息。本次使用的就是青藏高原的边界shp，因此不用选择特别的region，直接应用即可。

import shapefile
from matplotlib.path import Path
from matplotlib.patches import PathPatch
from shapely.geometry import Point as ShapelyPoint
from shapely.geometry import Polygon as ShapelyPolygon
from collections import Iterable
def shp2clip(originfig,ax,region_shpfile,clabel=None,vcplot=None):
    sf = shapefile.Reader(region_shpfile)
    for shape_rec in sf.shapeRecords():
    	#if shape_rec.record[0] == 1:#Hunan Sheng
        vertices = []
        codes = []
        pts = shape_rec.shape.points
        prt = list(shape_rec.shape.parts) + [len(pts)]
        for i in range(len(prt) - 1):
            for j in range(prt[i], prt[i+1]):
                vertices.append((pts[j][0], pts[j][1]))
            codes += [Path.MOVETO]
            codes += [Path.LINETO] * (prt[i+1] - prt[i] -2)
            codes += [Path.CLOSEPOLY]
        clip = Path(vertices, codes)
        clip = PathPatch(clip, transform=ax.transData)
    if vcplot:
        if isinstance(originfig,Iterable):
            for ivec in originfig:
                ivec.set_clip_path(clip)
        else:
            originfig.set_clip_path(clip)
    else:
        for contour in originfig.collections:
            contour.set_clip_path(clip)
    if  clabel:
        clip_map_shapely = ShapelyPolygon(vertices)
        for text_object in clabel:
            if not clip_map_shapely.contains(ShapelyPoint(text_object.get_position())):
                text_object.set_visible(False)    
    return clip

在这里插入图片描述
青藏高原边界shp利用arcmap打开如上所示

或者利用geopandas打开文件查询属性信息

import fiona
import geopandas
shp_path = 'E:/shapefile/青藏高原/青藏高原边界数据总集/TPBoundary_new(2021)/'
shp = geopandas.read_file(shp_path + 'TPBoundary_new(2021).shp')

在这里插入图片描述

二、绘制白化图形

1.引入库

代码如下（示例）：

import cartopy.crs as ccrs
from matplotlib.path import Path
from matplotlib.patches import PathPatch
import matplotlib.pyplot as plt
from osgeo import gdal
import numpy as np
import cartopy.crs as ccrs
import shapefile
import os
from cartopy.io.shapereader import Reader
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
import numpy
import pandas as pd
import xarray as xr
import matplotlib as mpl
from matplotlib.font_manager import FontProperties#自己设定字体的时候运用
import netCDF4 as nc
import cartopy.feature as cfeature
import cartopy.mpl.ticker as cticker
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER

2.读入数据并处理指标

代码如下（示例）：

data=xr.open_dataset('E:/CORDEX-TP/hfls_1981-2005.nc')
cclm_hf=data.cclm_hf
wrf_hf=data.wrf_hf
had_hf=data.had_hf
remo_hf=data.remo_hf
cclm_hf_clim=np.nanmean(cclm_hf,axis=0)
wrf_hf_clim=np.nanmean(wrf_hf,axis=0)
had_hf_clim=np.nanmean(had_hf,axis=0)
remo_hf_clim=np.nanmean(remo_hf,axis=0)

该处使用的url网络请求的数据。

3.绘制图形

代码如下（示例）：

leftlon=69
rightlon=105
lowerlat=24
upperlat=38
def contour_map(fig,img_extent,spec1,spec2,var,clevels):
    fig.set_extent(img_extent, crs=ccrs.PlateCarree())
    fig.set_xticks(np.arange(leftlon,rightlon+spec1,spec1), crs=ccrs.PlateCarree())
    fig.set_yticks(np.arange(lowerlat,upperlat+spec2,spec2), crs=ccrs.PlateCarree())
    lon_formatter = cticker.LongitudeFormatter()
    lat_formatter = cticker.LatitudeFormatter()
    fig.xaxis.set_major_formatter(lon_formatter)
    fig.yaxis.set_major_formatter(lat_formatter)
    fig.tick_params(labelsize=10)
    labels = fig.get_xticklabels() + fig.get_yticklabels()
    [label.set_fontname('Arial') for label in labels]
    ac=fig.contourf(lon,lat, var,clevels , cmap=cmap,norm = norm,
                        extend = 'both', transform=ccrs.PlateCarree() )
    fig.tick_params(labelsize=18)
    fig.add_geometries(Reader('E:/shapefile/青藏高原/青藏高原边界数据总集/TPBoundary_new(2021)/TPBoundary_new(2021).shp').geometries(),
                    crs=ccrs.PlateCarree(), facecolor='none',edgecolor='k',linewidth=1)
    return ac
    #%%绘图主程序
#公共设置
proj = ccrs.PlateCarree(central_longitude=85)
img_extent = [leftlon, rightlon, lowerlat, upperlat]
spec1=10
spec2=10
lon = np.arange(69.75, 140.5, 0.25)
lat = np.arange(14.75, 55.5, 0.25)
fig = plt.figure(figsize=[20,20],dpi=180)

# 采样简单的colormap线性映射关系https://www.cnblogs.com/cainiao-chuanqi/p/11301471.html
cdict = ['#191970','#6495ED','#024CEB','#87CEFA', '#E1FFFF', 
         '#FFFAFA', '#02BBA9', '#65FF00', '#FEFF00', '#FF8800', '#D40608','#BA55D3','#9400D3']
cmap = ListedColormap(cdict)
clevels1=[-150,-75,-25,-10,-5,0.0,5,10,25,50,75,150,200]#
norm = BoundaryNorm(clevels1, cmap.N)
#原文链接：https://blog.csdn.net/qq_35189715/article/details/109288053
#https://www.heywhale.com/mw/project/626965dcaf519200174517d7
#原文链接：https://blog.csdn.net/qq_35189715/article/details/109288053
#https://www.heywhale.com/mw/project/626965dcaf519200174517d7
#绘制子图
shp1='E:/shapefile/青藏高原/青藏高原边界数据总集/TPBoundary_new(2021)/TPBoundary_new(2021).shp'
#子图1
ax1=plt.subplot(5,4,1,projection=ccrs.PlateCarree())
ax1.set_title('(a) ERA5 DJF',loc='left',fontsize=12)#分别修改不同的季节并修改绘图区间
cf1=contour_map(ax1,img_extent,spec1,spec2,era5_hfls_DJF,clevels1)
ax1.set_xticks([])
ax1.set_yticks([])
clip1=shp2clip(ax1, ax1,shp1)
#...

#利用返回的ac值绘制共同的colorbar公用colorbar
核心点在contour_map中返回ac才能顺利地进行colorbar绘制

ac=contour_map(ax20,img_extent,spec1,spec2,era5_hfls_DJF,clevels1)
clip=shp2clip(ac, ax20,shp1)
l,b,w,h = 0.25, 0.08, 0.5, 0.02
rect = [l,b,w,h]#位置[左,下,右,上]
cbar_ax = fig.add_axes(rect)
cb = fig.colorbar(ac, cax = cbar_ax,orientation='horizontal')#spacing='proportional'
cb.set_label('Units:W/(m**2)',fontsize=18)
plt.suptitle(f'Surface latent heat flux',fontsize=32,x=0.5,y=0.92)
plt.savefig("hfls1.png", dpi=180)