python对tiff格式的反演结果进行分级，并绘制不同年份不同等级的面积占比热力图_gdal 热力图

在进行遥感反演之后，我们常常遇到需要对反演结果进行分级、统计分析，例如NDVI、RESI等等，常常利用自然断点分级等方法进行分级，然后进行统计，例如统计面积，今天的方法，主要是利用python读取多年的反演结果数据，然后按你需要的方式进行分级，最后统计每一级的面积占总面积的比值，以热力图的形式绘制结果。

1、读取tiff

import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from osgeo import gdal
from matplotlib.ticker import PercentFormatter
from matplotlib.font_manager import FontProperties

import matplotlib
matplotlib.rc("font",family='SimSun',size=7)

reclass_path = r'D:\work\ZR_reclass.tif'

resi_paths = [
r'D:\work\esi_2019.tif',
r'D:\work\esi_2020.tif',
r'D:\work\esi_2021.tif',
r'D:\work\esi_2022.tif'
]

2、分级

#按0.2的间隔进行分级，数值范围为0-1

thresholds = [0,0.2, 0.4, 0.6, 0.8, 1]
num_classes = len(thresholds)-1

3、计算每个年份的各级别面积占比

area_counts = np.zeros((len(resi_paths), num_classes))

for i, resi_path in enumerate(resi_paths):
resi_dataset = gdal.Open(resi_path)
resi_data = resi_dataset.ReadAsArray()

for j in range(num_classes):
mask = np.logical_and(resi_data >= thresholds[j], resi_data < thresholds[j + 1])
area_counts[i, j] = np.sum(mask)

area_counts = area_counts * 900

total_area_per_image = np.sum(area_counts, axis=1)
total_area_percent_per_image = total_area_per_image / np.sum(total_area_per_image)

area_percent_per_range = np.sum(area_counts, axis=0) / np.sum(area_counts)
area_percent_per_class_per_image = area_counts / np.sum(area_counts, axis=1).reshape(len(resi_paths), 1)

4、绘图并保存图片为300dpi的tiff图片

resi_ranges = [0,0.2, 0.4, 0.6, 0.8, 1]

#fig, ax = plt.subplots(figsize=(3.543, 2.165), dpi=300)
fig, ax = plt.subplots(figsize=(6.5, 2.165), dpi=300)
cmap = 'RdBu_r'
#cmap ='coolwarm'
# 绘制热力图，将其置于最底层
im = ax.imshow(np.flipud(area_percent_per_class_per_image.T), cmap=cmap, aspect='auto')

for i in range(len(resi_paths)):
for j in range(num_classes-1, -1, -1):
text = ax.text(i, num_classes-j-1, format(area_percent_per_class_per_image[i, j], '.3%'),
ha="center", va="center", color="k", weight='bold', fontsize=7)


font_properties = FontProperties(fname=r'C:\Windows\Fonts\simsun.ttc', size=7)

#plt.xlabel('年份', fontproperties=font_properties)
plt.ylabel('生态指数数值范围', fontproperties=font_properties)
plt.colorbar(im, label='面积占比', orientation='vertical', format=PercentFormatter(1.0))

years = [2019, 2020, 2021, 2022]
years_with_suffix = [str(year) + '年' for year in years]

plt.tick_params(axis='x', length=3) # 调整 x 轴刻度线的长度为 3
plt.tick_params(axis='y', length=3) # 调整 y 轴刻度线的长度为 3

plt.xticks(range(len(resi_paths)), years_with_suffix, fontproperties=font_properties)
ax.set_yticks(np.arange(num_classes))
labels = ['0.2', '0.4', '0.6', '0.8','1']
plt.yticks(np.arange(num_classes)-0.5, labels[::-1], fontproperties=font_properties, rotation=0)

for i in range(len(total_area_per_image)):
print(f"景影像 {i + 1}:")
print(f"总面积: {total_area_per_image[i]} 平方公里")
total_percentage = np.sum(area_percent_per_class_per_image[i])
for j in range(num_classes):
print(f"分级范围 {j}: 面积占比 {area_percent_per_class_per_image[i, j]*100:.4f}%")
print(f"所有分级的面积占比总和: {total_percentage*100:.2f}%")
print()

plt.savefig(r'D:\work\热力图.tif', dpi=300, format='tiff', bbox_inches='tight')
plt.show()

这样就得到了一张显示不同年份不同等级的面积占比热力图