Python地理数据处理 27：基于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结果如下：