【跟着SCI学作图】Cartopy和Matplotlib fill_betweenx可视化tif数据

Cartopy和Matplotlib fill_betweenx可视化tif数据

在这里插入图片描述

01 引言:

接着昨天的推文在参考文献【African Lightning and its Relation to Rainfall and Climate Change in a Convection‐Permitting Mode】中的图表的排版布局的基础上,加入自己的思考。在空间数据的的上方和右侧加入折线图,右上角加入频率直方图,最终结果如下图所示。

02 结果如下:

在这里插入图片描述

03 代码如下:

# -*- encoding: utf-8 -*-
'''
@File    :   wwlln_plot.py
@Time    :   2022/03/09 14:16:28
@Author  :   HMX 
@Version :   1.0
@Contact :   kzdhb8023@163.com
'''
# here put the import lib
import os
import time
import cartopy.crs as ccrs
import cmaps
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from cartopy.io.shapereader import Reader
from cartopy.mpl.ticker import LatitudeFormatter, LongitudeFormatter
from matplotlib.ticker import MultipleLocator
import rioxarray as rxr

def cm2inch(value): 
    return value/2.54

t1 = time.time()

size1 = 10.5
fontdict = {'weight': 'bold','size':size1,'color':'k','family':'SimHei'}
mpl.rcParams.update(
    {
    'text.usetex': False,
    'font.family': 'stixgeneral',
    'mathtext.fontset': 'stix',
    "font.family":'serif',
    "font.size": size1,
    "mathtext.fontset":'stix',
    "font.serif": ['Times New Roman'],
    }
    )


tif_path=r'E:\Project\XINAN\WWLLN_NEW\CG\ymean_CG.tif'
shp_path = r'E:\Project\XINAN\SHP\xinan_Dissolve.shp'
proj=ccrs.PlateCarree()
extent = [96,112,20,36]
fig = plt.figure(figsize=(cm2inch(16),cm2inch(16)),dpi=100)

ax = fig.add_axes([0.12,0.15,0.5,0.5],projection = proj)#c
ax1 = fig.add_axes([0.12,0.73,0.5,0.25])#a
ax2 = fig.add_axes([0.70,0.15,0.25,0.5])#d
ax3 = fig.add_axes([0.12,0.08,0.5,0.02])#colorbar
ax4 = fig.add_axes([0.70,0.73,0.25,0.25])#b

ax.set_xticks(np.arange(extent[0], extent[1] + 1, 4), crs = proj)
ax.set_yticks(np.arange(extent[-2], extent[-1] + 1, 4), crs = proj)
ax.xaxis.set_major_formatter(LongitudeFormatter(zero_direction_label=False))
ax.yaxis.set_major_formatter(LatitudeFormatter())
ax.set_extent(extent, crs=ccrs.PlateCarree())
ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=False,
linewidth=.8, linestyle=':', color='k', alpha=0.5)

xminorticks = MultipleLocator(1)
yminorticks = MultipleLocator(1)
ax.xaxis.set_minor_locator(xminorticks)
ax.yaxis.set_minor_locator(yminorticks)

# 绘制矢量边界
ax.add_geometries(Reader(shp_path).geometries(), proj, edgecolor='k', facecolor='none',alpha=1, linewidth=0.65)

# 绘制tif
tif = rxr.open_rasterio(tif_path)
lon,lat = np.meshgrid(tif['x'],tif['y'])
lev=[1,5,8,12,17,24,34,49]
cf=ax.contourf(lon,lat,tif[0],levels=lev,extend='both',transform=ccrs.PlateCarree(),cmap=cmaps.NCV_jaisnd)
b=plt.colorbar(cf,shrink=0.93,orientation='horizontal',extend='both',pad=0.035,aspect=30,ticks=lev,cax = ax3)
b.ax.set_xlabel(r'云地闪频次/$\mathrm{(fl·a^{-1})}$',fontdict = fontdict)

# 绘制经纬度曲线图
data = tif[0].values/(11.028*11.028)
x,y = tif['x'].values,tif['y'].values

def betweenx():
    # 计算25th 75th
    q1s, q4s, mus = [], [], []
    m,n = data.shape
    for i in range(m):
        line = data[i,:]
        res = pd.DataFrame({'cg':line})
        res = res[res.cg>0]        
        mu = res.describe().values[1][0]
        mus.append(mu)
        q1, q4 =res.describe().values[4][0], res.describe().values[6][0]
        q1s.append(q1)
        q4s.append(q4)
    return(q1s, q4s, mus)

def betweeny():
    # 计算25th 75th
    q1s, q4s, mus = [], [], []
    m,n = data.shape
    for i in range(n):
        line = data[:,i]
        res = pd.DataFrame({'cg':line})
        res = res[res.cg>0]
        mu = res.describe().values[1][0]
        mus.append(mu)
        q1, q4 =res.describe().values[4][0], res.describe().values[6][0]
        q1s.append(q1)
        q4s.append(q4)
    return(q1s, q4s, mus)

q1s_x, q4s_x, mus_x = betweenx()
q1s_y, q4s_y, mus_y = betweeny()

ax1.plot(x,mus_y)
ax1.fill_between(x,q1s_y,q4s_y, alpha=0.2,label = ' 25th to 75th',color = 'tab:blue')
ax2.plot(mus_x,y)
ax2.fill_betweenx(y,q1s_x,q4s_x, alpha=0.2,label = ' 25th to 75th',color = 'tab:blue')

ax1.set_xlim(96,112)
ax1.set_ylim(0,0.15) 
ax2.set_ylim(20,36)
ax2.set_xlim(0,0.15)

ax1.set_xticks(np.arange(96,113,4))
ax2.set_yticks(np.arange(20,37,4))
ax1.set_xticklabels([])
ax2.set_yticklabels([])
ax1.set_ylabel('云地闪密度/$\mathrm{(fl·km^{-2}·a^{-1})}$',fontdict = fontdict)
ax2.set_xlabel('云地闪密度/$\mathrm{(fl·km^{-2}·a^{-1})}$',fontdict = fontdict)
xminorticks1= MultipleLocator(1)
yminorticks1= MultipleLocator(1)
ax1.xaxis.set_minor_locator(xminorticks1)
ax2.yaxis.set_minor_locator(yminorticks1)

# 绘制频率直方图
fre = (tif[0].values[~np.isnan(tif[0].values)]).flatten()
dfre = pd.DataFrame({'box':fre,'cg':fre})
for i in range(0,50):
    dfre.box[(dfre.cg>=i)&((dfre.cg<(i+1)))] = i
dsfinal = dfre.groupby('box').count()
vcount = dsfinal.cg.sum()
dsfinal.cg = dsfinal.cg/vcount*100

ax4.bar(dsfinal.index,dsfinal.cg,width = 1,edgecolor = 'white')
ax4.set_xlabel(r'云地闪频次/$\mathrm{(fl·a^{-1})}$',fontdict = fontdict)
ax4.set_ylabel('频数/%',fontdict = fontdict,labelpad = -.5)
ax4.set_yticks(np.arange(0,13,4))
ax4.set_xticks(np.arange(0,51,10))
ax4.set_xlim(0,50)

# 设置序号
ax.set_title('c    ',loc = 'right',y = 0.9)
ax1.set_title('a    ',loc = 'right',y = 0.8)
ax2.set_title('d    ',loc = 'right',y = 0.9)
ax4.set_title('b    ',loc = 'right',y = 0.8)

plt.tight_layout()
plt.savefig(r'D:\公众号\N0.26\cg.png',dpi=600)
plt.show()    
t2 = time.time()
print('共计用时{:.2f}秒'.format(t2-t1))

​如果对你有帮助的话,请‘点赞’、‘收藏’,‘关注’,你们的支持是我更新的动力。
欢迎关注公众号【森气笔记】。

  • 4
    点赞
  • 15
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值