Python 计算土地分类的基本步骤

土地分类是将遥感影像中的像素划分为不同土地覆盖类型（如森林、草地、水体等）的过程。Python 中常用的土地分类工具包括 rasterio（用于读取栅格数据）和 numpy（用于数组操作）。以下是土地分类的基本步骤：

土地分类的基本步骤

1. 读取分类数据：使用 rasterio 读取土地分类栅格文件。
2. 统计各类面积：计算每个土地覆盖类型的像元数量，并根据像元大小转换为实际面积。
3. 计算面积占比：计算每个土地覆盖类型的面积占总面积的比例。
4. 可视化结果：使用 matplotlib 绘制饼图和柱状图，展示土地覆盖类型的面积分布。
5. 保存结果：将统计结果保存为 CSV 文件，将图表保存为 PNG 文件。

示例代码注释说明

import rasterio  # 导入 rasterio 库，用于读取栅格数据
import numpy as np  # 导入 numpy 库，用于数组操作
import pandas as pd  # 导入 pandas 库，用于数据处理
import matplotlib.pyplot as plt  # 导入 matplotlib 库，用于绘图
import os  # 导入 os 库，用于文件路径操作
import sys  # 导入 sys 库，用于系统相关操作
import io  # 导入 io 库，用于编码设置

# 设置标准输出的编码为 UTF-8
sys.stdout = io.TextIOWrapper(sys.stdout.buffer, encoding='utf-8')

def calculate_area_statistics(classification_file):
    """
    统计土地覆盖分类面积
    :param classification_file: 土地分类栅格文件路径
    :return: 统计结果的 DataFrame
    """
    print("开始统计土地覆盖分类面积...")
    
    # 定义土地覆盖类型的名称
    land_cover_names = {
        10: '树木覆盖',
        20: '灌木覆盖',
        30: '草地',
        40: '耕地',
        50: '建筑用地',
        60: '裸地/稀疏植被',
        70: '雪冰覆盖',
        80: '永久水体',
        90: '草本湿地',
        95: '红树林',
        100: '苔原'
    }
    
    # 打开土地分类栅格文件
    with rasterio.open(classification_file) as src:
        classification = src.read(1)  # 读取第一个波段的数据
        
        pixel_area = 100  # 每个像元的面积（单位：平方米）
        
        # 统计每个土地覆盖类型的像元数量
        unique_values, counts = np.unique(classification[classification != 0], return_counts=True)
        
        # 计算每个土地覆盖类型的面积（单位：平方公里）
        areas = (counts * pixel_area) / 1000000
        
        # 创建统计结果的 DataFrame
        results = pd.DataFrame({
            '类别编码': unique_values,
            '类别名称': [land_cover_names.get(x, f'未知类别{x}') for x in unique_values],
            '像元数量': counts,
            '面积(km²)': np.round(areas, 3),
            '占比(%)': np.round(areas / areas.sum() * 100, 2)
        })
        
        # 按面积从大到小排序
        results = results.sort_values('面积(km²)', ascending=False)
        
        # 设置 pandas 的显示选项，确保中文字符对齐
        pd.set_option('display.unicode.ambiguous_as_wide', True)
        pd.set_option('display.unicode.east_asian_width', True)
        
        # 打印统计结果
        print("\n土地覆盖分类面积统计结果：")
        print(results.to_string(index=False))
        
        # 创建绘图窗口
        plt.figure(figsize=(15, 6))
        
        # 绘制饼图
        plt.subplot(121)
        plt.pie(results['面积(km²)'], 
                labels=results['类别名称'],
                autopct='%1.1f%%', 
                startangle=90)
        plt.title('ESA WorldCover 土地覆盖类型面积比例')
        plt.axis('equal')
        
        # 绘制柱状图
        plt.subplot(122)
        bars = plt.bar(results['类别名称'], results['面积(km²)'])
        plt.title('ESA WorldCover 土地覆盖类型面积分布')
        plt.xlabel('土地覆盖类型')
        plt.ylabel('面积(km²)')
        plt.xticks(rotation=45, ha='right')
        
        # 在柱状图上添加数值标签
        for bar in bars:
            height = bar.get_height()
            plt.text(bar.get_x() + bar.get_width()/2., height,
                    f'{height:.3f}',
                    ha='center', va='bottom')
        
        # 调整布局
        plt.tight_layout()
        
        # 保存图表为 PNG 文件
        output_png = os.path.join(os.path.dirname(output_csv), 'ESA_WorldCover_面积统计.png')
        plt.savefig(output_png, dpi=300, bbox_inches='tight')
        print(f"\n图表已保存至：{output_png}")
        
        # 显示图表
        plt.show()
        
        # 返回统计结果
        return results

if __name__ == "__main__":
    # 设置 matplotlib 的中文字体和符号显示
    plt.rcParams['font.sans-serif'] = ['SimHei']
    plt.rcParams['axes.unicode_minus'] = False
    
    # 定义输入文件和输出路径
    classification_file = r'<土地覆盖分类数据路径>'
    output_csv = r'<统计结果输出路径(CSV格式)>'
    
    try:
        # 创建输出目录（如果不存在）
        os.makedirs(os.path.dirname(output_csv), exist_ok=True)
        
        # 调用函数统计土地覆盖分类面积
        results = calculate_area_statistics(classification_file)
        
        # 将统计结果保存为 CSV 文件
        results.to_csv(output_csv, index=False, encoding='utf-8-sig')
        print(f"\n统计结果已导出到：{output_csv}")
        
    except Exception as e:
        # 捕获并打印错误信息
        print(f"处理过程中出现错误: {str(e)}")