1. 批量裁剪進(jìn)階
代碼描述：遍歷a文件夾下的所有tif影像，并使用每個(gè)a文件夾中的tif影像對(duì)b文件夾下的所有tif影像進(jìn)行裁剪。裁剪后的柵格將以兩個(gè)tif文件進(jìn)行組合命名，并保存到另一個(gè)文件夾中。
# -*- 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"
# 設(shè)置工作環(huán)境
arcpy.env.workspace = a_folder
# 獲取a文件夾下的所有tif影像
a_rasters = arcpy.ListRasters("*", "TIF")
# 創(chuàng)建輸出文件夾
if not os.path.exists(output_folder):
os.makedirs(output_folder)
# 遍歷a文件夾下的每個(gè)tif影像
for indexa, a_raster in enumerate(a_rasters):
# 構(gòu)建輸出文件名前綴
a_name = os.path.splitext(a_raster)[0]
a_raster = a_folder + "\\" + a_name + ".tif"
# 設(shè)置工作環(huán)境為b文件夾
arcpy.env.workspace = b_folder
# 獲取b文件夾下的所有tif影像
b_rasters = arcpy.ListRasters("*", "TIF")
# 遍歷b文件夾下的每個(gè)tif影像
for indexb, b_raster in enumerate(b_rasters):
# 構(gòu)建輸出文件名
output_name = a_name + "_" + os.path.splitext(b_raster)[0] + ".tif"
# 構(gòu)建輸出路徑
output_path = os.path.join(output_folder, output_name)
# 設(shè)置工作環(huán)境為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. 統(tǒng)計(jì)運(yùn)算
獲取柵格數(shù)據(jù)的平均值，并輸出程序運(yùn)行進(jìn)度：
# -*- 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!")
程序運(yùn)行進(jìn)度：
3. 柵格批量縮小n倍
某文件夾中包含多個(gè)子文件夾，如：“2003clip”, “2004clip”, “2005clip”, “2006clip”, “2007clip”, “2008clip”, “2009clip”， 每個(gè)子文件夾中包含多個(gè)tif圖像，現(xiàn)在需要對(duì)這些圖像乘以縮放因子，如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")]
# 遍歷每個(gè)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對(duì)象
raster = arcpy.Raster(input_tif)
# 過濾像元值大于10000和小于0的像元
filtered_raster = arcpy.sa.Con((raster >= 0) & (raster <= 10000), raster)
# 對(duì)過濾后的像元乘以0.0001
scaled_raster = filtered_raster * 0.0001
# 保存輸出tif文件
scaled_raster.save(output_tif)
print folder_name + "OK!!!"
print "Finish!!!"
4. 建立屬性表（簡(jiǎn)化、普適）
之前已經(jīng)寫過如何建立柵格屬性表：Python地理數(shù)據(jù)處理 十九：arcpy批量處理數(shù)據(jù)之為柵格數(shù)據(jù)建立屬性表并導(dǎo)出
現(xiàn)在，結(jié)合之前的，寫一個(gè)更加簡(jiǎn)單普適的代碼，方便今后使用：文章來源：http://www.zghlxwxcb.cn/news/detail-445371.html
# -*- coding: utf-8 -*-
import arcpy
import os
arcpy.CheckOutExtension("Spatial")
# 設(shè)置工作空間
arcpy.env.workspace = r"D:\dataset"
# 獲取文件夾中的所有柵格數(shù)據(jù)
raster_list = arcpy.ListRasters()
# 檢查是否有至少一個(gè)柵格數(shù)據(jù)
if len(raster_list) < 1:
print("文件夾中至少需要有一個(gè)柵格數(shù)據(jù)才能建立屬性表。")
exit()
# 循環(huán)處理每個(gè)柵格數(shù)據(jù)
for raster in raster_list:
raster_path = os.path.join(arcpy.env.workspace, raster)
# 建立屬性表
arcpy.BuildRasterAttributeTable_management(raster_path)
print "屬性表已建立完成: {}。".format(raster)
print("屬性表已建立完成。")
5. 計(jì)算土地利用未變化區(qū)域（LUCC）
土地利用數(shù)據(jù)每年都在變化，如何判斷20年的土地利用變化數(shù)據(jù)中，哪些地方是未變化的？這是需要解決的一個(gè)問題，該代碼實(shí)現(xiàn)以上內(nèi)容，將變化的數(shù)據(jù)設(shè)置為100。
數(shù)據(jù)來源：GEE11：2個(gè)土地覆蓋數(shù)據(jù)（LUCC）分享和下載文章來源地址http://www.zghlxwxcb.cn/news/detail-445371.html
# -*- coding: utf-8 -*-
import arcpy
import os
arcpy.CheckOutExtension("Spatial")
# 設(shè)置工作空間
arcpy.env.workspace = r"D:\dataset"
# 獲取文件夾中的所有tif影像
tif_list = arcpy.ListRasters("*.tif")
# 檢查是否有至少2個(gè)tif影像
if len(tif_list) < 2:
print("文件夾中至少需要有2個(gè)tif影像才能進(jìn)行計(jì)算。")
exit()
# 選擇第一個(gè)tif影像作為基準(zhǔn)影像
base_tif = tif_list[0]
base_tif_path = os.path.join(arcpy.env.workspace, base_tif)
# 讀取基準(zhǔn)影像的像素?cái)?shù)組
base_raster = arcpy.Raster(base_tif_path)
base_array = arcpy.RasterToNumPyArray(base_raster)
# 循環(huán)處理每個(gè)tif影像
for tif in tif_list[1:]:
tif_path = os.path.join(arcpy.env.workspace, tif)
# 讀取當(dāng)前影像的像素?cái)?shù)組
current_raster = arcpy.Raster(tif_path)
current_array = arcpy.RasterToNumPyArray(current_raster)
# 對(duì)比基準(zhǔn)數(shù)組和當(dāng)前數(shù)組的像素值
comparison_array = (base_array == current_array)
# 將不相等的像素值置為100
result_array = current_array.copy()
result_array[~comparison_array] = 100
# 更新基準(zhǔn)數(shù)組為結(jié)果數(shù)組
base_array = result_array
print ("正常處理{}".format(tif))
# 創(chuàng)建輸出柵格
output_raster = arcpy.NumPyArrayToRaster(result_array, arcpy.Point(base_raster.extent.XMin, base_raster.extent.YMin),
base_raster.meanCellWidth, base_raster.meanCellHeight)
# 設(shè)置輸出路徑和名稱
output_path = os.path.join(arcpy.env.workspace, "unchanged_landuse.tif")
# 保存輸出柵格
output_raster.save(output_path)
print("未發(fā)生變化的土地利用柵格數(shù)據(jù)圖像已生成")
到了這里，關(guān)于Python地理數(shù)據(jù)處理 22：基于arcpy批量操作（四）的文章就介紹完了。如果您還想了解更多內(nèi)容，請(qǐng)?jiān)谟疑辖撬阉鱐OY模板網(wǎng)以前的文章或繼續(xù)瀏覽下面的相關(guān)文章，希望大家以后多多支持TOY模板網(wǎng)！
