Arcpy批量处理已矫正的worldclim2.1未来气候数据
- 1. 写在前面
- 2.实现代码
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、投影为albersprint("=========================投影为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、重采样为1000mprint("=========================重采样为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结果如下：
