【跟着SCI学作图】Python matplotlib绘制带有边界的散点图

前段时间看文献【Smartforests Canada: A Network of Monitoring Plots for Forest Management Under Environmental Change】看到了这么一幅图，可以展示不同森林区域的气温降水数据，并且还可以展示站点数据。由此想借鉴模仿一下。

由于无法获取文中数据，本例运用ERA5气象数据，和中科院的森林类型数据；并根据城市经纬度提取ERA5数据模拟站点数据。
最终结果如下：

代码如下，详见代码注释：

# -*- encoding: utf-8 -*-
'''
@File    :   prep_temp.py
@Time    :   2022/05/12 10:23:43
@Author  :   HMX 
@Version :   1.0
@Contact :   kzdhb8023@163.com
'''

# here put the import lib
import matplotlib as mpl
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from osgeo import gdal
from scipy.spatial import ConvexHull
import alphashape
from descartes import PolygonPatch
import matplotlib.colors as colors



def tif2arr(tif):
    ds = gdal.Open(tif)
    arr = ds.ReadAsArray().flatten()
    return arr


def pplot(data,ls='-',ec = 'k',label='label'):
    points = data.iloc[:,:-1].values
    alpha_shape = alphashape.alphashape(points,0)
    ax.add_patch(PolygonPatch(alpha_shape, alpha=1,fc = 'none',ec = ec,linestyle = ls,label = label))


def ExValues2Points(tif,df):
    arr = []
    dr = gdal.Open(tif)
    transform = dr.GetGeoTransform()
    x_origin = transform[0]
    y_origin = transform[3]
    pixel_width = transform[1]
    pixel_height = transform[5]
    for i in range(len(df)):
        x = df['经度'].values[i]
        y = df['纬度'].values[i]
        x_offset = int((x - x_origin) / pixel_width)
        y_offset = int((y - y_origin) / pixel_height)
        dt = dr.ReadAsArray()[x_offset, y_offset]
        arr.append(dt)
    return arr


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


# 数据读取
prep= r"E:ProjectCODEOtherZwenyidataexEXprep.tif"
temp= r'E:ProjectCODEOtherZwenyidataexEXtemp2020_07.tif'
lcpt= r'E:ProjectCODEOtherZwenyidataLCExmerge_Resample.tif'

df = pd.DataFrame({'tt':tif2arr(temp),'pp':tif2arr(prep),'lc':tif2arr(lcpt)})
nvlist = df.iloc[0,:]
for i in range(len(nvlist)):
    df = df[df.iloc[:,i] != nvlist[i]]
df.tt = df.tt-272.15
res = df[(df.lc>=1)&(df.lc<=5)]


fig,ax = plt.subplots(1,1,figsize = (8,6))
im = plt.hexbin(res.tt,res.pp,gridsize=20,mincnt=1,cmap = 'plasma_r')
ax.set_xlabel('TEMP/℃')
ax.set_ylabel('PREP/mm')
ax.set_title('西南地区不同林型气温降水图',fontdict=fontdict,fontsize = 18)
ax.set_xlim(5,30)
ax.set_ylim(500,4000)
cb = plt.colorbar()
cb.set_label('Count',fontsize = 10.5)

# 森林类型数据获取
res1 = df[df.lc==1]
res2 = df[df.lc==2]
res3 = df[df.lc==3]
res4 = df[df.lc==4]
res5 = df[df.lc==5]

# 森林类型外框绘制
pplot(res1,ls = ':',label = '灌木林')
pplot(res2,ls = '--',label = '季风雨林')
pplot(res3,ls = '-',label = '常绿阔叶林')
pplot(res4,ls = '-.',label = '常绿针叶林')
pplot(res5,ls = (0, (8, 10)),label = '混交林')

# 站点数据的获取
dfstd = pd.read_excel(r'E:ProjectCODEOtherZwenyidatastationxinan.xlsx',header=0)#读取站点经纬度
std = pd.DataFrame({'prep':ExValues2Points(prep,dfstd),'temp':ExValues2Points(temp,dfstd)})
std.temp = std.temp-272.15
std = std.round(2)
std = pd.concat([dfstd,std],axis=1)

# 绘制站点
clist = ['k','b','r','g']
for i in range(len(std)):
    plt.scatter(std.temp.values[i],std.prep.values[i],marker='o',facecolor = clist[i],edgecolors='k',s = 100,label = std['省会'].values[i])
ax.legend(prop = fontdict)
plt.tight_layout()
plt.savefig('pt.png',dpi = 600)
plt.show()

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