Python气象处理绘图第八弹–中国区域极端气温指标

努力努力再努力搬砖

已于 2023-12-14 18:31:20 修改

阅读量1.4k

点赞数 5

分类专栏：大气科学文章标签： python 数学建模开发语言 matplotlib

于 2023-06-08 01:09:59 首次发布

本文链接：https://blog.csdn.net/weixin_43347581/article/details/131083929

版权

大气科学专栏收录该内容

2 篇文章 1 订阅

订阅专栏

Python气象处理绘图第八弹–中国区域极端气温指标

pcolormesh绘制

Python气象处理绘图第八弹–中国区域极端气温指标
前言
一、极端指标的计算
- 1.引入库
- 2.读入数据并处理
二、绘制图形
总结

前言

根据ETCCDI中定义的极端气温指标，可以计算20年内tmax，tmin的极端气温指标：TXx，TXn，TNx，TNn，ID，FD，SU，TR等等。并可以绘制极端气温指标的多年平均的空间分布，主要分为简单的绘图，添加网格线和添加主、次要刻度三种子图。

一、极端指标的计算

1.引入库

代码如下（示例）：

import warnings
warnings.filterwarnings("ignore")
import xarray as xr
import numpy as np
import pandas as pd

2.读入数据并处理

代码如下（示例）：

#open file
#years:1986-2005
path='E:/dataset/CN05/'
data1=xr.open_dataset(path+'CN05.1_Tmax_1961_2021_daily_025x025.nc')
tmax=data1.tmax.loc['19860101':'20051231']
data2=xr.open_dataset(path+'CN05.1_Tmin_1961_2021_daily_025x025.nc')
tmin=data2.tmin.loc['19860101':'20051231']
#units
tmax = tmax.assign_attrs(units="degC")
tmin = tmin.assign_attrs(units="degC")
#delete the Jan 2nd
tmax=tmax.sel(time=~((tmax.time.dt.month == 2) & (tmax.time.dt.day == 29)))
tmin=tmin.sel(time=~((tmin.time.dt.month == 2) & (tmin.time.dt.day == 29)))
#Absolute indices
tmax_y=tmax.groupby(tmax.time.dt.year).mean()
txx = tmax.groupby(tmax.time.dt.year).max()
txn = tmax.groupby(tmax.time.dt.year).min()
tnx = tmin.groupby(tmin.time.dt.year).max()
tnn = tmin.groupby(tmin.time.dt.year).min()
#Threshold indices
T = np.array(tmax).reshape((20,365,163,283))
tr = np.nansum(T>25,axis=1)
T = np.array(tmin).reshape((20,365,163,283))
fd = np.nansum(T<0,axis=1)
su = np.nansum(T>25,axis=1)
#insert nan
fd1=np.where(tmax_y.isnull(),np.nan,fd)
su1=np.where(tmax_y.isnull(),np.nan,su)
tr1=np.where(tmax_y.isnull(),np.nan,tr)
#climatology
txx_c=np.nanmean(txx,axis=0)
txn_c=np.nanmean(txn,axis=0)
tnn_c=np.nanmean(tnn,axis=0)
tnx_c=np.nanmean(tnx,axis=0)
fd_c=np.nanmean(fd1,axis=0)
tr_c=np.nanmean(tr1,axis=0)
su_c=np.nanmean(su1,axis=0)
tasin=np.stack((txx_c,txn_c,tnx_c,tnn_c,fd_c,tr_c,su_c))

二、绘制图形

1.简要图

代码如下（示例）：

import cmaps
c1=cmaps.temp_19lev
#不带投影的格式
import matplotlib.ticker as ticker
fig, axs = plt.subplots(nrows=2, ncols=2, figsize=(12, 8))
lon = np.arange(69.75, 140.5, 0.25)
lat = np.arange(14.75, 55.5, 0.25)
x, y = np.meshgrid(lon, lat)
for i in range(4):
    r, c = divmod(i, 2)
    im=axs[r, c].pcolormesh(x, y, tasin[i,:,:],vmin=-30,vmax=30,cmap=c1)
# 创建格式化器，将轴标签转换为度数符号
xticks_fmt = ticker.StrMethodFormatter('{x:.0f}°E')
yticks_fmt = ticker.StrMethodFormatter('{x:.0f}°N')
for ax in axs.flat:
    ax.xaxis.set_major_formatter(xticks_fmt)
    ax.yaxis.set_major_formatter(yticks_fmt)
# 设置 x 和 y 轴的间距，并将 Locator 添加到每个轴上
xlocator = ticker.MultipleLocator(base=20.0)
ylocator = ticker.MultipleLocator(base=10.0)
for ax in axs.flat:
    ax.xaxis.set_major_locator(xlocator)
    ax.yaxis.set_major_locator(ylocator)
# 调整颜色条大小并添加至子图末尾
cbar = fig.colorbar(im, ax=axs.ravel().tolist(), shrink=0.8, extend='both',location='right')
cbar.ax.tick_params(labelsize=12)
plt.savefig('Figures/etccdi1.eps')
plt.savefig('Figures/etccdi1.png',dpi=500)

在这里插入图片描述

2.带网格线

代码如下（示例）：

import matplotlib.ticker as ticker
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.io.shapereader import Reader
from cartopy.mpl.gridliner import LongitudeFormatter, LatitudeFormatter
import cmaps
c1=cmaps.temp_19lev
# 设置字体类型
plt.rcParams['font.family'] = 'Times New Roman'
# 将全部图形元素的字体大小设置为18
plt.rc('font', size=18)
fig, axs = plt.subplots(nrows=2, ncols=2, figsize=(12, 8), subplot_kw={'projection': ccrs.PlateCarree()})
lon = np.arange(69.75, 140.5, 0.25)
lat = np.arange(14.75, 55.5, 0.25)
x, y = np.meshgrid(lon, lat)
shpp = 'E:/shapefile/国界省界十段线南海诸岛/'
titles = ['TXx', 'TXn', 'TNx', 'TNn']
for i in range(4):
    r, c = divmod(i, 2)
    im = axs[r, c].pcolormesh(x, y, tasin[i,:,:], cmap=c1)
    cbar = fig.colorbar(im, ax=axs[r, c], shrink=0.6, extend='both')
    # 将colorbar分成8个刻度显示
    tick_locator = ticker.MaxNLocator(nbins=8)
    cbar.locator = tick_locator
    cbar.update_ticks()
    cbar.ax.tick_params(labelsize=18)
    # 添加图例
    axs[r, c].text(0.1, 0.9, f'{chr(97+i)})', transform=axs[r, c].transAxes, fontsize=18)
    # 增加子图副标题
    axs[r, c].set_title(titles[i], loc='left', fontsize=18)
for ax in axs.flat:
    ax.add_geometries(Reader(shpp+'china10.shp').geometries(),
                    crs=ccrs.PlateCarree(), facecolor='none',edgecolor='k',linewidth=0.5)
    gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True, linewidth=0.5, color='gray', alpha=0.5, linestyle='--',
                      x_inline=False, y_inline=False, clip_path=None)
    gl.xlabels_top = False
    gl.ylabels_right = False
    gl.xlabel_style = {'size': 12}
    gl.ylabel_style = {'size': 12}
    # 设置刻度线样式和间距
    gl.ylocator = ticker.FixedLocator(np.arange(22, 56, 8))
    gl.yformatter = LatitudeFormatter(degree_symbol='°')
    gl.xlocator = ticker.FixedLocator(np.arange(70, 142, 16))
    gl.xformatter = LongitudeFormatter(degree_symbol='°')

plt.tight_layout()#使得子图清晰展示
plt.savefig('Figures/etccdi2.eps')
plt.savefig('Figures/etccdi2.png',dpi=500)
plt.show()

输出eps的时候会遇到warning：

The PostScript backend does not support transparency; partially transparent artists will be rendered opaque.
eps打开之后发现后面三个子图缺少了刻度、图例和子标题星信息
png格式如下：

3.带主次刻度

代码如下（示例）：

#外刻度线
import matplotlib.ticker as ticker
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.io.shapereader import Reader
from cartopy.mpl.gridliner import LongitudeFormatter, LatitudeFormatter
import cmaps
c1=cmaps.temp_19lev
# 设置字体类型
plt.rcParams['font.family'] = 'Times New Roman'
# 将全部图形元素的字体大小设置为18
plt.rc('font', size=18)
proj=ccrs.PlateCarree()
fig, axs = plt.subplots(nrows=2, ncols=2, figsize=(12, 8), subplot_kw={'projection': ccrs.PlateCarree()})
lon = np.arange(69.75, 140.5, 0.25)
lat = np.arange(14.75, 55.5, 0.25)
x, y = np.meshgrid(lon, lat)
shpp='E:/shapefile/国界省界十段线南海诸岛/'
titles = ['TXx', 'TXn', 'TNx', 'TNn']
for i in range(4):
    r, c = divmod(i, 2)
    im = axs[r, c].pcolormesh(x, y, tasin[i,:,:], cmap=c1)
    cbar = fig.colorbar(im, ax=axs[r, c], shrink=0.6, extend='both')
    # 将colorbar分成8个刻度显示
    tick_locator = ticker.MaxNLocator(nbins=8)
    cbar.locator = tick_locator
    cbar.update_ticks()
    cbar.ax.tick_params(labelsize=18)
    # 添加图例
    axs[r, c].text(0.1, 0.9, f'{chr(97+i)})', transform=axs[r, c].transAxes, fontsize=20)
   # 增加子图副标题
    axs[r, c].set_title(titles[i], loc='left', fontsize=18)
for ax in axs.flat:
    ax.add_geometries(Reader(shpp+'china10.shp').geometries(),
                    crs=ccrs.PlateCarree(), facecolor='none',edgecolor='k',linewidth=0.5)
    # 设置大刻度和小刻度
    ax.set_xticks(np.arange(70, 145 , 20), crs=proj)
    ax.set_xticks(np.arange(70, 145 , 10), minor=True, crs=proj)
    ax.set_yticks(np.arange(15, 60 , 10), crs=proj)
    ax.set_yticks(np.arange(15, 55 , 5), minor=True, crs=proj)
    # 利用Formatter格式化刻度标签
    ax.xaxis.set_major_formatter(LongitudeFormatter())
    ax.yaxis.set_major_formatter(LatitudeFormatter())
plt.tight_layout()#使得子图清晰展示
plt.savefig('Figures/etccdi3.eps')
plt.savefig('Figures/etccdi3.png',dpi=500)

在这里插入图片描述

代码中使用 divmod() 函数将索引 i 除以 2，其中 r 为商，c 为余数。这个函数返回一个元组，包含商和余数。
在这个代码中，r 和 c 分别代表行号和列号，用于确定当前子图在子图数组中的位置。如果 i 是偶数，divmod() 函数返回 (i/2, 0)，表示子图在当前行第一列的位置；如果 i 是奇数，divmod() 函数返回 (i//2, 1)，表示子图在当前行第二列的位置。!](https://img-blog.csdnimg.cn/c9d04b544ea042f49d23d2abf2f5ccf5.png)
divmod

4.绘制一张子图

#外刻度线
import matplotlib.ticker as ticker
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.io.shapereader import Reader
from cartopy.mpl.gridliner import LongitudeFormatter, LatitudeFormatter
import cmaps
c1=cmaps.temp_19lev
# 设置字体类型
plt.rcParams['font.family'] = 'Times New Roman'
# 将全部图形元素的字体大小设置为18
plt.rc('font', size=18)
proj=ccrs.PlateCarree()
fig, axs = plt.subplots(nrows=1, ncols=1, figsize=(12, 8), subplot_kw={'projection': ccrs.PlateCarree()})
lon = np.arange(69.75, 140.5, 0.25)
lat = np.arange(14.75, 55.5, 0.25)
#x, y = np.meshgrid(lon, lat)
shpp='F:/shapefile/国界省界十段线南海诸岛/'
titles = ['TXx', 'TXn', 'TNx', 'TNn']
#im = axs.pcolormesh(x, y, tasin[0,:,:], cmap=c1)
im = axs.contourf(lon, lat, tasin[0,:,:], cmap=c1)
cbar = fig.colorbar(im, ax=axs, shrink=0.6, extend='both')
# 将colorbar分成8个刻度显示
tick_locator = ticker.MaxNLocator(nbins=8)
cbar.locator = tick_locator
cbar.update_ticks()
cbar.ax.tick_params(labelsize=18)
# 添加图例
axs.text(0.1, 0.9, f'{chr(97+0)})', transform=axs.transAxes, fontsize=20)
# 增加子图副标题
axs.set_title(titles[0], loc='left', fontsize=18)
#添加文件的shp数据
axs.add_geometries(Reader(shpp+'china10.shp').geometries(),
                crs=ccrs.PlateCarree(), facecolor='none',edgecolor='k',linewidth=0.5)
# 设置大刻度和小刻度
axs.set_xticks(np.arange(70, 145 , 20), crs=proj)
axs.set_xticks(np.arange(70, 145 , 10), minor=True, crs=proj)
axs.set_yticks(np.arange(15, 60 , 10), crs=proj)
axs.set_yticks(np.arange(15, 55 , 5), minor=True, crs=proj)
# 利用Formatter格式化刻度标签
axs.xaxis.set_major_formatter(LongitudeFormatter())
axs.yaxis.set_major_formatter(LatitudeFormatter())
plt.tight_layout()#使得子图清晰展示
plt.savefig('Figures/etccdi4.eps')
plt.savefig('Figures/etccdi4.png',dpi=500)

在这里插入图片描述
利用sel函数可以索引感兴趣的区域，例如下面可以剔除闰年2月29日的数据

tmax=tmax.sel(time=~((tmax.time.dt.month == 2) & (tmax.time.dt.day == 29)))
t1=tmax.sel(time=(tmax.time.dt.season == 'MAM') )

总结

本章主要简单介绍了绝对阈值的计算，此外还可以应用xclim进行其它极端气候指标的计算，详情之后有时间再更新。值得注意的是，这里剔除2月29日的方法比较简单高效；再进行网格线eps图片的绘制中，会丢失很多信息（暂时还没有解决）。强烈推荐第三种画法，但是还仍有南海子图的绘制未完善。

努力努力再努力搬砖

关注

5
点赞
踩
39

收藏

觉得还不错? 一键收藏
9
评论
Python气象处理绘图第八弹–中国区域极端气温指标

根据ETCCDI中定义的极端气温指标，可以计算20年内tmax，tmin的极端气温指标：TXx，TXn，TNx，TNn，ID，FD，SU，TR等等。并可以绘制极端气温指标的多年平均的空间分布，主要分为简单的绘图，添加网格线和添加主、次要刻度三种子图。本章主要简单介绍了绝对阈值的计算，此外还可以应用xclim进行其它极端气候指标的计算，详情之后有时间再更新。值得注意的是，这里剔除2月29日的方法比较简单高效；再进行网格线eps图片的绘制中，会丢失很多信息（暂时还没有解决）。
复制链接

扫一扫