Python地理数据处理 28:基于Arcpy批量操作实现——按属性提取和分区统计

1. 批量按属性提取

# -*- coding: cp936 -*-
"""
PROJECT_NAME: ArcPy 
FILE_NAME: batch_attribute_extract 
AUTHOR: JacksonZhao
DATE: 2025/04/05
"""
import arcpy
import os

arcpy.CheckOutExtension('Spatial')

# Set the workspace environment
input_folder = r"C:\Users\YiPei\Desktop\2Reclassify"
output_folder = r"C:\Users\YiPei\Desktop\2Reclassify_reclassify4"

# Check if the output folder exists, and create it if it doesn't
if not os.path.exists(output_folder):
    os.makedirs(output_folder)

# Get the list of all raster data
arcpy.env.workspace = input_folder
raster_list = arcpy.ListRasters("*", "TIF")

# Define the list of conditions
conditions = ["Value = 0", "Value = 1", "Value = 2", "Value = 3"]

# Iterate through each raster and perform the extraction for each condition
for raster in raster_list:
    print "Processing raster: %s" % raster
    for condition in conditions:
        # Define the output raster name with the condition appended
        base_name = os.path.splitext(raster)[0]
        output_raster_name = "%s_%s.tif" % (base_name, condition.replace(' ', ''))
        out_raster = os.path.join(output_folder, output_raster_name)
        
        # Execute the raster extraction by attributes
        arcpy.gp.ExtractByAttributes_sa(raster, condition, out_raster)
        print "Extracted raster saved to: %s" % out_raster

2. 批量分区统计(最大值、最小值和像元个数等)

# -*- coding: cp936 -*-
# -*- coding: utf-8 -*-
"""
PROJECT_NAME: ArcPy 
FILE_NAME: batch_UHI 
AUTHOR: JacksonZhao
DATE: 2025/02/23
"""
import arcpy
import os

# Check out the Spatial Analyst extension
arcpy.CheckOutExtension('Spatial')

# Set the input and output folders
input_folder = r"C:\Users\YiPei\Desktop\2Reclassify_reclassify4"
output_folder = r"C:\Users\YiPei\Desktop\2Reclassify_reclassify4_table"

# Set the path to the shapefile
shapefile = r"C:\Users\YiPei\Desktop\GEE_data\ChinaAlbers\ChinaAlbers.shp"

# Check if the output folder exists, and create it if it doesn't
if not os.path.exists(output_folder):
    os.makedirs(output_folder)

# Get the list of all raster files in the input folder
arcpy.env.workspace = input_folder
raster_list = arcpy.ListRasters("*", "TIF")

# Iterate through each raster file
for raster in raster_list:
    print "Processing raster: %s" % raster
    
    # Define the output CSV file name
    base_name = os.path.splitext(raster)[0]
    output_csv = os.path.join(output_folder, base_name + "_ZonalStats.csv")
    
    # Perform Zonal Statistics as Table
    arcpy.gp.ZonalStatisticsAsTable(
        shapefile,
        "name",
        raster,
        output_csv,
        "DATA",
        "ALL"
    )
    
    print "Zonal Statistics results saved to: %s" % output_csv

### ArcGIS 中均值方差的计算方法 在 ArcGIS 中,可以通过多种方式来计算栅格数据(如 DEM)中的均值方差。以下是几种常见的方法: #### 使用【计算统计数据】工具 对于单个栅格文件,可以直接使用 ArcToolbox 中的数据管理工具下的【计算统计数据】功能[^1]。该工具会自动计算并更新栅格属性表中的统计信息,包括最小值、最大值、平均值以及标准偏差。 ```python import arcpy arcpy.CalculateStatistics_management(in_raster="path_to_your_dem.tif", ignore_values="", x_skip_factor=0, y_skip_factor=0, skip_nulls="SKIP_NULLS", area_of_interest="#") ``` 此命令适用于处理单一影像或少量影像的情况;当面对大量影像时,则需考虑编程批量执行上述操作。 #### 利用 R 语言实现批量化计算 针对多个 TIFF 文件组成的 DEM 数据集,利用 R 软件包 `raster` 及其相关函数可以高效完成多景影像上指定区域内高程均值与方差的求解工作[^2]。下面给出一段简单的脚本用于示范如何读取目录下所有的 .tif 文件,并对其逐个进行分析: ```r library(raster) # 设置输入输出路径 input_dir <- "D:/dem_data" output_csv <- "D:/results.csv" files_list <- list.files(path=input_dir,pattern='*.tif$',full.names=T) result_df <- data.frame() for(file in files_list){ r <- raster(file) # 提取感兴趣区(ROI),此处假设整个研究范围即为 ROI mean_val <- cellStats(r,'mean') var_val <- cellStats(r,'var') temp_df<-data.frame(filename=file,basename(basename(file)),mean_value=mean_val,variance_value=var_val) result_df=rbind(result_df,temp_df) } write.table(result_df,file=output_csv,row.names=F,sep=",") ``` 这段代码实现了遍历给定文件夹内的所有 tif 文件,依次获取各文件内数值型像元的整体平均数(mean) 变异数(variance),并将结果保存至 CSV 表格中以便后续查阅。 #### 结合空间分区要素类进行局部统计 如果希望按照特定的空间单元(比如行政边界、流域等)分别汇总内部各地形因子指标,则可通过叠加分析的方式先将目标矢量图层转成掩膜(mask),再调用条件运算符 (Con) 或者 Zonal Statistics as Table 工具来进行更精细层面的信息挖掘[^4]。 ```python from arcpy.sa import * mask_layer = "county_boundaries.shp" # 替换成实际使用的面状要素名称 elevation_ras = Raster("dem.tif") outZSaT = ZonalStatisticsAsTable(mask_layer,"FID", elevation_ras , "in_memory/zstat","DATA","MEAN;STD") ``` 以上 Python 语句展示了怎样基于县级行政区划轮廓作为参照系,对镶嵌于其中的不同海拔高度实施分组聚合,最终获得每县范围内地面起伏程度的具体描述——即均值(MEAN) 标准差(STD) 的组合表现形式。
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

Jackson的生态模型

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值