[python]matplotlib绘制遥感两波段散点图和密度散点图

用python绘制遥感两波段的散点图和密度散点图，包括1：1线、拟合直线、R2、RMSE等，自留用，密度散点图的colorbar图例似乎还有点问题，欢迎大家讨论。

#二维散点图
from scipy import stats
import numpy as np
import matplotlib.pyplot as plt
from osgeo import gdal

# 读取数据
dataset=gdal.Open('读取第一个数据')
band=dataset.GetRasterBand(1)
p1=band.ReadAsArray()
p1=p1.flatten()

dataset1=gdal.Open('读取第二个数据')
band1=dataset1.GetRasterBand(1)
p2=band1.ReadAsArray()
p2=p2.flatten()

#计算数据最大值
max1=np.max(p1)
max2=np.max(p2)
max_final=max(max1,max2)

fig, ax=plt.subplots()
plt.scatter(p1, p2, alpha=0.5,s=3)#画散点图
#设置坐标轴范围
plt.ylim([0,max_final])#设置坐标轴范围为0到最大值
plt.xlim([0,max_final])
plt.plot([0,max_final],[0,max_final],'k--',lw=1.0) # 绘制1:1线

#对散点数据进行线性拟合 获取斜率 截距 R2
slope, intercept, r_value, p_value, std_err = stats.linregress(p1, p2) #斜率 截距 R2

#画拟合线
X1 = np.arange (0,int(max_final),100)
Y1 = np.array([ intercept+ slope * p1 for p1 in X1])
plt.plot(X1,Y1,'y--',lw=1.0)

plt.rcParams['font.sans-serif'] = ['SimSun']#设置中文字符

plt.xlabel('x轴标题',font={'size':10})  # 横坐标轴标题
plt.ylabel('y轴标题',font={'size':10})  # 纵坐标轴标题

#写入公式以及R2
plt.text(max_final*0.05,max_final*0.95,'R$^2$=%s'%(np.around(r_value**2,2)),fontsize=12)
plt.text(max_final*0.05,max_final*0.9,'y=%s*x+%s'%(np.around(slope,2),np.around(intercept,2)),fontsize=12)

plt.title('标题',font={'size':17})
#设置没有上边框和右边框
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
plt.savefig('存储位置',dpi=250)
plt.show()

二维散点图

#二维密度散点图
from statistics import mean
import matplotlib.pyplot as plt

from osgeo import gdal
import numpy as np
from sklearn.metrics import explained_variance_score,r2_score,median_absolute_error,mean_squared_error,mean_absolute_error
from scipy import stats

# 读取数据
dataset=gdal.Open('数据1')
band=dataset.GetRasterBand(1)
p1=band.ReadAsArray()
p1=p1.flatten()

dataset1=gdal.Open('数据2')
band1=dataset1.GetRasterBand(1)
p2=band1.ReadAsArray()
p2=p2.flatten()

#计算数据最大值
max1=np.max(p1)
max2=np.max(p2)
max_final=max(max1,max2)

# ==========计算评价指标==========
BIAS = mean(p1 - p2)
MSE = mean_squared_error(p1, p2)
RMSE = np.power(MSE, 0.5)
R2 = r2_score(p1, p2)
MAE = mean_absolute_error(p1, p2)
EV = explained_variance_score(p1, p2)

print('==========算法评价指标==========')
print('BIAS:', '%.3f' % (BIAS))
print('Explained Variance(EV):', '%.3f' % (EV))
print('Mean Absolute Error(MAE):', '%.3f' % (MAE))
print('Mean squared error(MSE):', '%.3f' % (MSE))
print('Root Mean Squard Error(RMSE):', '%.3f' % (RMSE))
print('R_squared:', '%.3f' % (R2))

fig, ax=plt.subplots()

#对散点数据进行线性拟合 获取斜率 截距 R2
slope, intercept, r_value, p_value, std_err = stats.linregress(p1, p2) #斜率 截距 R2
# ===========Calculate the point density==========
xy = np.vstack([p1, p2])
z = stats.gaussian_kde(xy)(xy)
# ===========Sort the points by density, so that the densest points are plotted last===========
idx = z.argsort()
p1, p2, z = p1[idx], p2[idx], z[idx]


# edgecolor: 设置轮廓颜色的,none是无色
# 将c设置成一个y值列表并使用参数 cmap.cm.XX来使用某个颜色渐变(映射)
# https://matplotlib.org/examples/color/colormaps_reference.html 从这个网页可以获得colorbard的名称
# c=是设置数据点颜色的,在2.0的matplotlib中,edgecolor默认为'none'
# c=还可以使用RGB颜色,例如 c = (0,0,0.8) 分别表示 红绿蓝,值越接近0,颜色越深.越接近1,颜色越浅
plt.scatter(p1, p2, c=z, s=7, edgecolor='none', cmap='jet')
plt.plot([0, max_final], [0, max_final], 'k--', lw=1.0)  # 画的1:1线，线的颜色为black，线宽为1.0

#画拟合线
X1 = np.arange (0,int(max_final),100)
Y1 = np.array([ intercept+ slope * p1 for p1 in X1])
plt.plot(X1,Y1,'y--',lw=1.0)

plt.axis([0, max_final, 0, max_final])  # 设置线的范围

plt.rcParams['font.sans-serif'] = ['SimSun']#设置中文字符
plt.xlabel('x轴标题',font={'size':10})  # 横坐标轴标题
plt.ylabel('y轴标题',font={'size':10})  # 纵坐标轴标题
plt.title('标题',font={'size':17})

cb=plt.colorbar()

plt.text(max_final*0.05,max_final*0.95, '$N=%.f$' % len(p2), family = 'Times New Roman',fontsize=9) # text的位置需要根据x,y的大小范围进行调整。
plt.text(max_final*0.05,max_final*0.9, '$R^2=%.3f$' % R2, family = 'Times New Roman',fontsize=9)
plt.text(max_final*0.05,max_final*0.85, '$BIAS=%.4f$' % BIAS, family = 'Times New Roman',fontsize=9)
plt.text(max_final*0.05,max_final*0.8, '$RMSE=%.3f$' % RMSE, family = 'Times New Roman',fontsize=9)
plt.xlim(0, max_final)                                  # 设置x坐标轴的显示范围
plt.ylim(0, max_final)                                  # 设置y坐标轴的显示范围

#设置没有上边框和右边框
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
plt.savefig('存储位置',dpi=250)

plt.show()
plt.pause(20)            # 暂停3秒钟