Python绘制可以表示密度的散点图

  本文介绍基于Python语言的matplotlib模块，对Excel表格文件中的指定数据，加以密度散点图绘制的方法。

  首先，明确一下本文的需求。

  现有一个.csv格式的表格文件，其各列数据的开头部分如下图所示。其中，对于名称为26的这1列（左侧紫色框内数据），我们希望提取其数值等于1的所有行，并对这些行中的NIR_predict列与NIR_true列（右侧紫色框内数据）的数值加以密度散点图的绘制。

  明确了需求，即可开始代码的撰写。本文所用代码如下。

# -*- coding: utf-8 -*-
"""
Created on Mon Apr  1 12:14:38 2024

@author: fkxxgis
"""

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.stats import gaussian_kde

csv_file_path = r"E:\04_Reconstruction\99_MODIS\Train_Model_0715_Main_Combine.csv"
picture_file_path = r"E:\04_Reconstruction\99_MODIS\Scatter_result.png"

data = pd.read_csv(csv_file_path)

x = data[data["26"] == 1]['NIR_true']
y = data[data["26"] == 1]['NIR_predict']

xy = np.vstack([x,y])
z = gaussian_kde(xy)(xy)
idx = z.argsort()
plt.scatter(x, y, c = z, s = 10, cmap = "Spectral")
plt.colorbar()

plt.rc("font", family = "Times New Roman")
x_line = np.linspace(min(min(x), min(y)), max(max(x), max(y)), 100)
plt.plot(x_line, x_line, color='black', linestyle='--')
plt.xlabel('NIR_true')
plt.ylabel('NIR_predict')
plt.savefig(picture_file_path, dpi = 400)
plt.show()

  首先，我们通过import语句导入所需模块。其中，numpy用于数值计算，pandas用于数据处理，matplotlib.pyplot用于绘图，scipy.stats.gaussian_kde用于计算核密度估计。

  其次，通过csv_file_path定义待绘图的.csv格式文件的路径，通过picture_file_path定义存储所得图片结果的路径。

  随后，使用pd.read_csv()从.csv格式文件中读取数据，并存储在名为data的DataFrame中。通过筛选条件data["26"] == 1从DataFrame中获取符合指定条件的数据，并分别存储在x和y中。

  接下来，使用np.vstack()将x和y垂直堆叠为一个二维数组xy，并使用scipy.stats.gaussian_kde()计算二维数据的核密度估计值，并将其存储在z中；使用z.argsort()对z进行排序，返回索引值，并将其存储在idx中。使用plt.scatter()绘制散点图，其中x和y是散点的横纵坐标，c是颜色值，s是散点的大小，cmap是颜色映射，并使用plt.colorbar()添加颜色条。

  紧接着，使用plt.rc()设置字体为Times New Roman；随后，生成一条直线的横坐标范围，使用np.linspace()生成一系列横坐标值，并存储在x_line中；这些点将组成后续所得散点图中的x = y线；使用plt.plot()绘制直线，颜色为黑色，线型为虚线。此外，使用plt.xlabel()和plt.ylabel()添加x轴和y轴的标签，使用plt.savefig()将图形保存为指定路径的图片文件，设置dpi值为400。最后，使用plt.show()显示图形。

  执行上述代码，即可在结果文件夹中看到所得图片；如下图所示。

  可以看到，我们已经绘制得到了指定数据之间的密度散点图。当然，我这里所选色带，将密度较低的区域标记为红色系，密度较高的区域标记为了蓝色系，可能和一般情况下大家常用的色系相反——我是一开始选错了，后面也没有修改，这里大家理解即可；如果需要修改这个色系，大家修改上述代码中的cmap = "Spectral"部分即可。

  至此，大功告成。