批量裁剪
- 1. 批量裁剪进阶
- 2. 统计运算
- 3. 栅格批量缩小n倍
- 4. 建立属性表（简化、普适）
- 5. 计算土地利用未变化区域（LUCC）
1. 批量裁剪进阶
代码描述：遍历a文件夹下的所有tif影像，并使用每个a文件夹中的tif影像对b文件夹下的所有tif影像进行裁剪。裁剪后的栅格将以两个tif文件进行组合命名，并保存到另一个文件夹中。
# -*- coding: cp936 -*-
import arcpy
import os
import time
start = time.clock()
arcpy.CheckOutExtension("spatial")
arcpy.CheckOutExtension("spatial")
a_folder = r"D:\dataset\a"
b_folder = r"D:\dataset\b"
output_folder = r"D:\dataset\output_tif"
# 设置工作环境
arcpy.env.workspace = a_folder
# 获取a文件夹下的所有tif影像
a_rasters = arcpy.ListRasters("*", "TIF")
# 创建输出文件夹
if not os.path.exists(output_folder):
os.makedirs(output_folder)
# 遍历a文件夹下的每个tif影像
for indexa, a_raster in enumerate(a_rasters):
# 构建输出文件名前缀
a_name = os.path.splitext(a_raster)[0]
a_raster = a_folder + "\\" + a_name + ".tif"
# 设置工作环境为b文件夹
arcpy.env.workspace = b_folder
# 获取b文件夹下的所有tif影像
b_rasters = arcpy.ListRasters("*", "TIF")
# 遍历b文件夹下的每个tif影像
for indexb, b_raster in enumerate(b_rasters):
# 构建输出文件名
output_name = a_name + "_" + os.path.splitext(b_raster)[0] + ".tif"
# 构建输出路径
output_path = os.path.join(output_folder, output_name)
# 设置工作环境为a文件夹
arcpy.env.workspace = a_folder
b_raster = b_folder + "\\" + b_raster
# 裁剪b影像到输出路径
arcpy.gp.ExtractByMask_sa(b_raster, a_raster, output_path)
index = indexa + indexb + 1
print "{} have been completed and the current file is ".format(index) + output_name
end = time.clock()
print "All finish!!!"
2. 统计运算
获取栅格数据的平均值，并输出程序运行进度：
# -*- coding: utf-8 -*-
import arcpy
import os
import glob
from arcpy.sa import *
arcpy.CheckOutExtension("ImageAnalyst")
arcpy.CheckOutExtension("spatial")
input_folder = r"D:\dataset"
output_file = r'D:\dataset\Meandata.csv'
rasters = glob.glob(os.path.join(input_folder, "*.tif"))
where_clause = "VALUE = -32556"
total_rasters = len(rasters)
processed_rasters = 0
with open(output_file, 'w') as output:
for raster in rasters:
outSetNull = SetNull(raster, raster, where_clause)
meanValueInfo = arcpy.GetRasterProperties_management(outSetNull, 'MEAN')
meanValue = meanValueInfo.getOutput(0)
output.write(os.path.basename(raster).split('.')[0] + ',' + str(meanValue) + '\n')
processed_rasters += 1
progress = float(processed_rasters) / total_rasters * 100
print("Processed {} out of {} rasters. Progress: {:.2f}%".format(processed_rasters, total_rasters, progress))
print("All processing is done!")
程序运行进度：
3. 栅格批量缩小n倍
某文件夹中包含多个子文件夹，如：“2003clip”, “2004clip”, “2005clip”, “2006clip”, “2007clip”, “2008clip”, “2009clip”， 每个子文件夹中包含多个tif图像，现在需要对这些图像乘以缩放因子，如0.0001
# -*- coding: cp936 -*-
import arcpy
import os
arcpy.CheckOutExtension('Spatial')
# 输入文件夹和输出文件夹路径
input_folder = r"D:\Datasets"
output_folder = r"D:\Datasets\缩小10000倍"
# 遍历输入文件夹中的所有子文件夹
for folder_name in ["2003clip", "2004clip", "2005clip", "2006clip", "2007clip", "2008clip", "2009clip"]:
print folder_name
# 子文件夹路径
folder_path = os.path.join(input_folder, folder_name)
# 获取子文件夹中的所有tif文件
tif_files = [file for file in os.listdir(folder_path) if file.endswith(".tif")]
# 遍历每个tif文件
for tif_file in tif_files:
# 输入tif文件路径
input_tif = os.path.join(folder_path, tif_file)
# 输出tif文件路径
output_tif = os.path.join(output_folder, tif_file)
# 打开Raster对象
raster = arcpy.Raster(input_tif)
# 过滤像元值大于10000和小于0的像元
filtered_raster = arcpy.sa.Con((raster >= 0) & (raster <= 10000), raster)
# 对过滤后的像元乘以0.0001
scaled_raster = filtered_raster * 0.0001
# 保存输出tif文件
scaled_raster.save(output_tif)
print folder_name + "OK!!!"
print "Finish!!!"
4. 建立属性表（简化、普适）
之前已经写过如何建立栅格属性表：Python地理数据处理 十九：arcpy批量处理数据之为栅格数据建立属性表并导出
现在，结合之前的，写一个更加简单普适的代码，方便今后使用：
# -*- coding: utf-8 -*-
import arcpy
import os
arcpy.CheckOutExtension("Spatial")
# 设置工作空间
arcpy.env.workspace = r"D:\dataset"
# 获取文件夹中的所有栅格数据
raster_list = arcpy.ListRasters()
# 检查是否有至少一个栅格数据
if len(raster_list) < 1:
print("文件夹中至少需要有一个栅格数据才能建立属性表。")
exit()
# 循环处理每个栅格数据
for raster in raster_list:
raster_path = os.path.join(arcpy.env.workspace, raster)
# 建立属性表
arcpy.BuildRasterAttributeTable_management(raster_path)
print "属性表已建立完成: {}。".format(raster)
print("属性表已建立完成。")
5. 计算土地利用未变化区域（LUCC）
土地利用数据每年都在变化，如何判断20年的土地利用变化数据中，哪些地方是未变化的？这是需要解决的一个问题，该代码实现以上内容，将变化的数据设置为100。
数据来源：GEE11：2个土地覆盖数据（LUCC）分享和下载
# -*- coding: utf-8 -*-
import arcpy
import os
arcpy.CheckOutExtension("Spatial")
# 设置工作空间
arcpy.env.workspace = r"D:\dataset"
# 获取文件夹中的所有tif影像
tif_list = arcpy.ListRasters("*.tif")
# 检查是否有至少2个tif影像
if len(tif_list) < 2:
print("文件夹中至少需要有2个tif影像才能进行计算。")
exit()
# 选择第一个tif影像作为基准影像
base_tif = tif_list[0]
base_tif_path = os.path.join(arcpy.env.workspace, base_tif)
# 读取基准影像的像素数组
base_raster = arcpy.Raster(base_tif_path)
base_array = arcpy.RasterToNumPyArray(base_raster)
# 循环处理每个tif影像
for tif in tif_list[1:]:
tif_path = os.path.join(arcpy.env.workspace, tif)
# 读取当前影像的像素数组
current_raster = arcpy.Raster(tif_path)
current_array = arcpy.RasterToNumPyArray(current_raster)
# 对比基准数组和当前数组的像素值
comparison_array = (base_array == current_array)
# 将不相等的像素值置为100
result_array = current_array.copy()
result_array[~comparison_array] = 100
# 更新基准数组为结果数组
base_array = result_array
print ("正常处理{}".format(tif))
# 创建输出栅格
output_raster = arcpy.NumPyArrayToRaster(result_array, arcpy.Point(base_raster.extent.XMin, base_raster.extent.YMin),
base_raster.meanCellWidth, base_raster.meanCellHeight)
# 设置输出路径和名称
output_path = os.path.join(arcpy.env.workspace, "unchanged_landuse.tif")
# 保存输出栅格
output_raster.save(output_path)
print("未发生变化的土地利用栅格数据图像已生成")
