Python地理数据处理 21：基于arcpy批量操作（三）

1. 实现一（批量裁剪的高级用法）

实现将给定的 .shp 文件中的所有省份作为裁剪范围，对给定的 .tif 文件进行裁剪，输出所有省份的单独 .tif 文件：

# -*- coding: cp936 -*-
import arcpy

arcpy.CheckOutExtension('Spatial')

# 定义输入文件的路径
shp_path = r"C:\Users\map.shp"
tif_path = r"C:\Users\dem.tif"

# 读取文件
shp = arcpy.mapping.Layer(shp_path)
tif = arcpy.Raster(tif_path)

# 获取.tif 文件的空间参考信息
spatial_ref = tif.spatialReference

# 创建一个用于存储输出的文件夹
output_folder = r"C:\Users\output"
arcpy.CreateFolder_management(output_folder, "clipped_tifs")

# 循环遍历 .shp 文件中的每个省份
for row in arcpy.da.SearchCursor(shp, ["SHAPE@", "name"]):
    # 构建输出.tif 文件的路径
    output_tif_path = output_folder + "\\clipped_tifs\\" + row[1] + ".tif"
    
    # 使用.shp 文件对.tif 文件进行裁剪
    arcpy.Clip_management(tif_path, "#", output_tif_path, row[0], "#", "ClippingGeometry")

    # 打印输出信息
    print row

# 删除中间变量
del shp, tif

2. 实现二（栅格计算器求和）

实现对文件名前14个字符一样的tif图像进行栅格运算求和：
如：XXXX_XXX_2003.M01_Mean、XXXX_XXX_2003.M02_Mean、XXXX_XXX_2003.M03_Mean；XXXX_XXX_2004.M01_Mean、XXXX_XXX_2004.M02_Mean、XXXX_XXX_2004.M03_Mean；XXXX_XXX_2005.M01_Mean、XXXX_XXX_2005.M02_Mean、XXXX_XXX_2005.M03_Mean。。。。。。

# -*- coding: cp936 -*-
import arcpy
import os

arcpy.CheckOutExtension('Spatial')

# 设置工作空间
arcpy.env.workspace = r"D:\Datasets"

# 获取所有.tif文件
tif_list = arcpy.ListRasters("*", "TIF")

# 创建字典，用于存储同一年份的栅格图像
tif_dict = {}

# 遍历所有.tif文件，将同一年份的栅格图像添加到字典中
for tif in tif_list:
    year = tif.split("_")[2][:4]
    if year in tif_dict:
        tif_dict[year].append(tif)
    else:
        tif_dict[year] = [tif]

# 对每个年份的栅格图像进行栅格运算求和，并保存结果
for year in tif_dict:
    tif_sum = None
    for tif in tif_dict[year]:
        tif_path = os.path.join(arcpy.env.workspace, tif)
        tif_raster = arcpy.Raster(tif_path)
        if tif_sum is None:
            tif_sum = tif_raster
        else:
            tif_sum += tif_raster
    output_name = "D:\\Datasets\\TIFClip\\" + year
    output_path = os.path.join(arcpy.env.workspace, output_name)
    tif_sum.save(output_path)
    print year

3. 实现三（使用栅格计算器删除指定值）

1、如果我们想删除栅格图像的属性表中的某个空值或者指定值，如9999或52668，如果在arcgis中操作的话，可以使用栅格计算器，输入以下代码即可：

Con((input.tif != 9999) & (input.tif != 52668), input.tif)

2、也可以使用arcpy进行批量处理：

# -*- coding: cp936 -*-
import arcpy
import os

arcpy.CheckOutExtension("spatial")
# 定义输入和输出文件夹
input_folder = r"D:\Datasets"
output_folder = r"D:\Datasets\tifclipCon"

# 循环遍历文件夹下的所有tif文件
for filename in os.listdir(input_folder):
    if filename.endswith(".tif"):
        # 构建输入和输出路径
        input_path = os.path.join(input_folder, filename)
        output_path = os.path.join(output_folder, filename)

        # 使用Map Algebra表达式删除指定的属性表值
        expression = "Con((\"" + input_path + "\" != 9999) & (\"" + input_path + "\" != 52668), \"" + input_path + "\")"
        arcpy.gp.RasterCalculator_sa(expression, output_path)
        
        print filename

4. 实现四（读取两个文件夹中的tif文件进行相乘并保存）

读取两个文件夹中的tif文件，并对两个相对应的文件进行相乘：

# -*- coding: cp936 -*-
import arcpy
import os
import warnings
warnings.filterwarnings('ignore')

arcpy.CheckOutExtension('Spatial')

# 设置工作空间
arcpy.env.workspace = r"D:\RSE\dataset\1\output"

# 设置输入文件夹路径
t_folder = r"D:\RSE\dataset\t\output"
p_folder = r"D:\RSE\dataset\p\output"

# 设置输出文件夹路径
out_folder = r"D:\RSE\dataset\tbyp"

# 获取输入文件夹中以TEM开头的tif文件列表
t_list = arcpy.ListRasters("TEM*", "TIF")

# 遍历输入文件夹中的tif文件
for t_file in t_list:
    
    # 获取tif文件名的后三个字符
    suffix = t_file[3:6]
    
    # 构造PRE文件路径
    p_file = os.path.join(p_folder, "PRE" + suffix + ".tif")
    
    # 判断PRE文件是否存在
    if not arcpy.Exists(p_file):
        print("PRE file does not exist for " + t_file)
        continue

    print "*******************************************************************"
    # 构造输出文件路径
    out_file = os.path.join(out_folder, "output" + suffix + ".tif")
    
    # 使用栅格计算器对两个tif文件进行相乘
    expression = "Times('{}', '{}')".format(t_file, p_file)
    arcpy.gp.RasterCalculator_sa(expression, out_file)
    
    print "Processed " + out_file
   

本代码使用Times()函数对两个栅格进行相乘

5. 实现五（文件夹中所有tif文件相加）

实现文件夹中所有tif文件相加：

# -*- coding: cp936 -*-
import arcpy
import os

arcpy.CheckOutExtension('Spatial')

# 设置工作空间
arcpy.env.workspace = r"D:\RSE\dataset\TbyP"

# 设置输出文件夹路径
out_folder = r"D:\RSE\dataset\TbyP"

# 获取输入文件夹中的所有tif文件
tif_list = arcpy.ListRasters("*", "TIF")

# 如果没有tif文件，则输出提示信息
if not tif_list:
    print("There are no TIFF files in the input folder.")
else:
    # 构造栅格计算器表达式
    #expression = " + ".join(tif_list)
    expression = ' + '.join(['"' + tif + '"' for tif in tif_list])

    # 构造输出文件路径
    out_file = os.path.join(out_folder, "outputAll.tif")

    # 使用栅格计算器对所有tif文件进行相加
    arcpy.gp.RasterCalculator_sa(expression, out_file)

    print"The output file is located at " + out_file