好久没有更新实用教程啦,也刚好赶上确实有这方面的需要,于是“火”不停蹄跟大伙分享这个实现过程。
故事背景主要发生在昨天需要处理一些指标数据,需要利用自然断点法进行重分类;但是考虑数据量过多,亟需利用批量来进行处理。(至于为啥是自然断点法是因为要运用到地理探测器,目前网上有关于处理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)
。