【实用教程】基于Python的栅格数据批量重分类(自然断点法),几乎无需配置环境!

发布于:2025-09-03 ⋅ 阅读:(13) ⋅ 点赞:(0)

好久没有更新实用教程啦,也刚好赶上确实有这方面的需要,于是“火”不停蹄跟大伙分享这个实现过程。

故事背景主要发生在昨天需要处理一些指标数据,需要利用自然断点法进行重分类;但是考虑数据量过多,亟需利用批量来进行处理。(至于为啥是自然断点法是因为要运用到地理探测器,目前网上有关于处理0-1这种固定区间,关于运用模型迭代器进行处理的教程,但尝试了一下效果不太行),于是考虑转战运用Python进行处理。

第二波折,AI的代码教程推荐运用Arcpy、GDAL的思路,毫无疑问这思路是最快高效解决需求,因为能直接调用相关库,但是问题在于需要进行一下环境配置,如虚拟环境等;于是再次发问,不完全依赖这两种。

采用rasterio、jenkspy,刚好赶上数据的量不是很大(这两个库需要自己安装一下哈,非常方便,用pip install 就行),所有数据加起来大概50MB,其中也遇到一些报错,但好在能用AI及时修改,最后解决了这个问题,并且在ArcGIS中进行了研究发现结果与跑出来的情况一样的。

综上,当面临大量数据时做一遍没问题,其实如果立马去重分类可能当写完这篇推文的时候已经做完了,但是这就不能水这篇推文了呀doge;主要还是考虑到可能有小伙伴也会遇到类似的情况,更何况,地理探测器的分类情况还需要不断的尝试,这样运行的时候,小火也能摸鱼干其他事情啦,比如更一篇推文呢是吧。假如报错全文复制给AI让它帮你改,结束后也进行随机进行一下验证哈。觉得有需要可以转发给身边需要的小伙伴,祝大家开学快乐

下面附一下相关代码,改一下输入输出路径名就行运行的,点赞+收藏=提升成功率哈

import os
import numpy as np
import rasterio
import jenkspy
from glob import glob
def read_raster(raster_path):
    """读取栅格数据及元信息"""
    with rasterio.open(raster_path) as src:
        data = src.read(1)
        profile = src.profile
        nodata = src.nodata
    return data, profile, nodata
def calculate_jenks_breaks(data, num_classes, nodata):
    """计算自然断点"""
    if nodata is not None:
        valid_data = data[data != nodata]
    else:
        valid_data = data[~np.isnan(data)]
    
    if len(valid_data) < num_classes:
        raise ValueError(f"有效数据点({len(valid_data)})少于分类数({num_classes})")
    
    return jenkspy.jenks_breaks(valid_data, n_classes=num_classes)
def reclassify_raster(data, breaks, nodata):
    """重分类数据"""
    reclassified = np.zeros_like(data, dtype=np.int32)
    
    for i in range(1, len(breaks)):
        if i == 1:
            mask = (data >= breaks[i-1]) & (data <= breaks[i])
        else:
            mask = (data > breaks[i-1]) & (data <= breaks[i])
        
        if nodata is not None:
            mask = mask & (data != nodata)
        else:
            mask = mask & (~np.isnan(data))
            
        reclassified[mask] = i
    
    if nodata is not None:
        reclassified[data == nodata] = 0
    else:
        reclassified[np.isnan(data)] = 0
        
    return reclassified
def write_reclassified_raster(output_path, data, profile):
    """写入重分类结果"""
    profile.update(
        dtype=rasterio.int32,
        count=1,
        nodata=0
    )
    
    with rasterio.open(output_path, 'w', **profile) as dst:
        dst.write(data, 1)
def batch_process(input_dir, output_dir, num_classes=5):
    """批量处理栅格文件"""
    os.makedirs(output_dir, exist_ok=True)
    
    # 修复路径拼接错误:os.path.join.join → os.path.join
    raster_files = glob(os.path.join(input_dir, "*.tif")) + glob(os.path.join(input_dir, "*.tiff"))
    
    if not raster_files:
        print(f"在{input_dir}中未找到TIFF格式栅格文件")
        return
    
    for file in raster_files:
        try:
            print(f"处理文件: {os.path.basename(file)}")
            data, profile, nodata = read_raster(file)
            breaks = calculate_jenks_breaks(data, num_classes, nodata)
            print(f"自然断点: {breaks}")
            
            reclassified = reclassify_raster(data, breaks, nodata)
            base_name = os.path.splitext(os.path.basename(file))[0]
            output_path = os.path.join(output_dir, f"{base_name}_reclass.tif")
            
            write_reclassified_raster(output_path, reclassified, profile)
            print(f"已保存至: {output_path}\n")
            
        except Exception as e:
            print(f"处理{file}时出错: {str(e)}\n")
if __name__ == "__main__":
    # 注意:此处也修复了拼写错误 outout_FOLDER → OUTPUT_FOLDER
    INPUT_FOLDER = "输入路径名"
    OUTPUT_FOLDER = "输出路径名"
    NUM_CLASSES = 6
    
    batch_process(INPUT_FOLDER, OUTPUT_FOLDER, NUM_CLASSES)


网站公告

今日签到

点亮在社区的每一天
去签到