【python掩膜及多子图colorbar】

2024-07-16 13:39| 来源: 网络整理| 查看: 265

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

mask和colorbar设定 Python气象处理绘图（三）--青藏高原空间分布白化前言一、maskout子程序思路？二、绘制白化图形1.引入库2.读入数据并处理指标3.绘制图形总结

前言

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

绘制的图形：在这里插入图片描述

shp文件的读取，特征值选取，maskout的改写由于cordex数据是旋转坐标，因此需要进行插值成非旋转坐标的值才能够进行常规坐标分析，在后文会对IDW插值的运用进行更新多子图colorbar的设定核心：返回ac值。如果不返回ac的话，绘制的colorbar则无法适应子图。气象绘图cmap、cbar超详细版(附示例)

【本文地址】

公司简介

联系我们