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

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))
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163

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