代码收藏家 Python地理数据处理 27：基于Arcpy批量处理已矫正的worldclim2.1未来气候数据——投影、重采样、多波段拆分以及裁剪
Arcpy批量处理已矫正的worldclim2.1未来气候数据
1. 写在前面
前面我写了一篇关于如何使用ArcGIS自带的Python工具处理worldclim数据的多波段数据的文章，而这只是处理该数据的其中一步。要想得到满足要求的数据，还需要其他操作，依次为投影为指定投影坐标系（Albers）、重采样为1000m空间分辨率、多波段拆分以及裁剪。
今天我把以上所有的相关操作基于一套代码进行演示，可以实现一键操作。并且我使用的是已矫正的worldclim2.1未来气候数据，该数据排除了异常值。
2.实现代码
使用以上代码之前，你需要建立一个文件夹，例如：C:\Users\Jacksonzhao\Desktop\sspdata\ssp126。再将worldclim数据放入其中，这里我选择了worldclim的ssp126_2030s和ssp126_2050s两幅影像数据进行处理，之后便可以运行代码，喝一杯茶的功夫就好了：
# -*- coding: cp936 -*-
import arcpy
import shutil
from arcpy.sa import *
import os
import time
# 检查并获取 Spatial Analyst 扩展许可
arcpy.CheckOutExtension('Spatial')
start_time = time.time() # 开始时间
# *****************************************************************************************
# 基础路径和子文件夹
subfolder = 'ssp126' # 需要修改
base_path = os.path.join(r'C:\Users\Jacksonzhao\Desktop\sspdata', subfolder)
# 设置工作空间
arcpy.env.workspace = base_path
# 获取工作空间中的所有tif格式文件
raster_lists = arcpy.ListRasters("*.tif")
print(raster_lists)
for raster_list in raster_lists:
#print(raster_list)
raster_list = [raster_list]
# 检查是否找到了tif文件
if raster_list:
# 遍历文件列表
for raster in raster_list:
print("*****************************************************************************")
print("正在处理文件：{}************************************************".format(raster))
# 移除.tif扩展名
raster_name = os.path.splitext(raster)[0]
#print(raster_name) # 打印不含扩展名的文件名
# 子文件夹模板列表
subfolders = [
'{}/Albers998m'.format(raster_name),
'{}/Albers1000m'.format(raster_name),
'{}/Albers1000m_19bio'.format(raster_name),
'{}_Finish'.format(raster_name)
]
print("=========================文件夹创建=========================")
# 遍历子文件夹模板列表，为每个raster_name创建子文件夹
for subfolder_template in subfolders:
subfolder_path = subfolder_template.format(raster_name=raster_name)
folder_path = os.path.join(base_path, subfolder_path)
# 检查文件夹是否已经存在
if not os.path.exists(folder_path):
os.makedirs(folder_path)
print("文件夹 '{}' 已创建".format(os.path.basename(folder_path)))
else:
print("文件夹 '{}' 已存在".format(os.path.basename(folder_path)))
print("所有文件夹创建完毕，耗时：{:.2f}分====================================".format((time.time() - start_time)/60))
# 1、投影为albers
print("=========================投影为Albers=========================")
arcpy.env.workspace = base_path
#rasters = arcpy.ListRasters("*","tif")
for raster in raster_list:
out = os.path.join(base_path, raster_name, 'Albers998m', raster) # 输出路径
print(os.path.join(base_path, raster_name, 'Albers998m', raster))
arcpy.ProjectRaster_management(raster, out, "PROJCS['Asia_North_Albers_WGS84_LCR',GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Albers'],PARAMETER['False_Easting',0.0],PARAMETER['False_Northing',0.0],PARAMETER['Central_Meridian',105.0],PARAMETER['Standard_Parallel_1',25.0],PARAMETER['Standard_Parallel_2',47.0],PARAMETER['Latitude_Of_Origin',0.0],UNIT['Meter',1.0]]")
print("操作完毕，耗时：{:.2f}分====================================".format((time.time() - start_time)/60))
# 2、重采样为1000m
print("=========================重采样为1000m=========================")
arcpy.env.workspace = os.path.join(base_path, raster_name, 'Albers998m')
input_folder = os.path.join(base_path, raster_name, 'Albers998m')
output_folder = os.path.join(base_path, raster_name, 'Albers1000m')
resampling_resolution = "1000" # 设置重采样分辨率
tif_list = arcpy.ListFiles("*.tif")
# 循环处理每个tif文件
for tif_file in tif_list:
#构建输入和输出路径
input_path = os.path.join(input_folder, tif_file)
output_path = os.path.join(output_folder, tif_file)
# 执行重采样
arcpy.Resample_management(input_path, output_path, resampling_resolution, "NEAREST")
print(output_path)
print("操作完毕，耗时：{:.2f}分====================================".format((time.time() - start_time)/60))
# 3、拆分为单波段
print("=========================多波段拆分为单波段=========================")
arcpy.env.workspace = os.path.join(base_path, raster_name, 'Albers1000m')
output_directory = os.path.join(base_path, raster_name, 'Albers1000m_19bio')
print(output_directory)
# 获取栅格列表
raster_list = arcpy.ListRasters("*.tif")
# 遍历栅格集合
for raster in raster_list:
# 获取栅格数据集的全路径
raster_full_path = os.path.join(arcpy.env.workspace, raster)
# 创建栅格对象
raster_obj = arcpy.Raster(raster_full_path)
# 获取波段数量
number_of_bands = raster_obj.bandCount
# 构建基本的文件名
base_filename = "wc2.1BCC_" + raster_name
# 遍历所有波段
for i in range(1, number_of_bands + 1):
# 提取单波段并创建栅格图层
band_name = "band_" + str(i)
arcpy.MakeRasterLayer_management(raster_full_path, band_name, "", "", i)
# 保存单波段影像
output_path = os.path.join(output_directory, base_filename + "_Bio{0}.tif".format(i))
arcpy.CopyRaster_management(band_name, output_path)
print("Saved:", output_path)
print("操作完毕，耗时：{:.2f}分====================================".format(time.time() - start_time))
# 4、裁剪
print("=========================影像裁剪=========================")
finish_folder = os.path.join(base_path, '{}_Finish'.format(raster_name))
arcpy.env.workspace = os.path.join(base_path, raster_name, 'Albers1000m_19bio')
rasters = arcpy.ListRasters("*", "tif")
mask ="C:/Users/Jacksonzhao/Desktop/CSCD/Data/Bio_Topo/envi1km/wc2.1_30s_bio_1.tif"
for raster in rasters:
print(raster)
out = os.path.join(finish_folder, raster)
arcpy.gp.ExtractByMask_sa(raster, mask, out)
print("Clip_" + raster + "has done!")
print("操作完毕，耗时：{:.2f}分====================================".format((time.time() - start_time)/60))
# 5、删除文件夹
folder_paths = [
os.path.join(base_path, raster_name, 'Albers998m'),
os.path.join(base_path, raster_name, 'Albers1000m_19bio')
]
# 检查并删除文件夹
for folder_path in folder_paths:
if os.path.isdir(folder_path):
shutil.rmtree(folder_path)
print("文件夹 : {} 及其内容已删除".format(folder_path))
else:
print("文件夹 : {} 不存在".format(folder_path))
end_time = time.time() # 开始时间
print("所有操作完毕，耗时：{:.2f}分*****************************************".format((end_time - start_time)/60))
else:
print("没有找到tif文件。")
print("所有操作完毕!!!，耗时：{:.2f}分*****************************************".format((end_time - start_time)/60))
处理结果：
BCC_ssp126_2030s_Bio1结果如下：
作者：Jackson的生态模型
