【ArcGIS自定义脚本工具】批量执行栅格计算器


系列文章目录: ArcGIS自定义脚本工具


一、功能介绍

在这里插入图片描述

  • 功能
    将多个栅格文件按照某以代数表达式执行栅格计算器(Raster Calculator) 工具,并将输出栅格文件保存至指定的文件夹中。该工具接受代数表达式的动态输入。通过修改代数表达式可以实现批量对栅格执行四则运算(加、减、乘、除)、计算掩膜、填充空值(固定值、插值)、设为空等操作。
  • 关于代数表达式
    提供的代数表达式需要满足栅格计算器工具代数表达式的基本语法,并且包含字符串"{A}“,字符串”{A}"代表输入栅格在代数表达式中的位置。在用例章节提供了几个案例。
    在这里插入图片描述

二、脚本代码

#!/usr/bin/python
# -*- coding: UTF-8 -*-
import arcpy
import os
import time

rasters = arcpy.GetParameterAsText(0)
expression = arcpy.GetParameterAsText(1) # expression including {A}
out_path = arcpy.GetParameterAsText(2)
prefix = arcpy.GetParameterAsText(3)

rasters = rasters.split(";")
nums = len(rasters)
num = 1
for raster in rasters:
    raster_path = raster.replace("'","")
    raster_dir,raster_name = os.path.split(raster_path)
    arcpy.env.workspace = raster_dir
    out_raster = os.path.join(out_path, prefix + raster_name)
    exp = expression.replace('{A}', r'"{0}"'.format(raster_name))
    arcpy.AddMessage(exp)
    if not os.path.exists(out_raster):
        try:
            arcpy.gp.RasterCalculator_sa(exp, out_raster)
            arcpy.AddMessage("%d/%d | %s completed" % (num, nums, out_raster))
        except Exception as err:
            arcpy.AddMessage("%d/%d | %s errored, %s" % (num, nums, out_raster, err))
    else:
        arcpy.AddMessage("%d/%d | %s already exists" % (num, nums, out_raster))
    num = num + 1

三、工具参数

在这里插入图片描述
在这里插入图片描述

ParameterExplanationData type
rasters待执行栅格计算器的多个栅格文件Raster Dataset
expression栅格计算器工具要执行的代数表达式。该表达式需要满足栅格计算器工具中代数表达式的基本语法,并且包含字符串"{A}“,字符串”{A}"代表输入栅格在代数表达式中的位置。String
out_path输出栅格文件的文件夹Folder
prefix输出栅格文件的文件名前缀,该栅格的文件名为"prefix_<旧文件名>.tif",该参数默认为"cal_"String

关于参数expression,提供了一些常用的示例

1.加

{A} + 10

2.减

{A} - 273.15

3.乘

{A} * 0.0001

4.除

{A} / 10000

5.四则运算

({A} * 3 + 4)/5

6.条件函数-计算掩膜

Con({A} > 0,1,0})

7.条件函数-计算植被覆盖度

Con({A}<0.1,0,Con({A}>=0.8,1,({A}-0.1)/0.7))

8.条件函数-固定值填充空值

Con(IsNull({A}), 0, {A})

9.条件函数-插值填充空值

Con(IsNull({A}), FocalStatistics({A}, NbrRectangle(5,5, "CELL"), "MEAN"), {A})

10.设为空-将满足条件的栅格赋值为空值

SetNull({A}<0,{A})


四、用例

通过修改参数expression,对栅格批量执行某一代数表达式。

4.1 批量单位换算

  • 输入文件:某区域的陆面温度(LST)栅格数据,单位为开氏温度
  • 处理目标:批量将输入文件的单位换算为摄氏温度
{A} - 273.15

在这里插入图片描述
在这里插入图片描述

4.2 批量计算植被覆盖度

假设我们根据NDVI的累计频率分布,已经确定了像元二分模型中的两个阈值分别为0.1、0.85。则可以通过将该工具中的expression设置为:

Con({A}<0.1,0,Con({A}>=0.85,1,({A}-0.1)/0.75))

在这里插入图片描述
在这里插入图片描述

4.3 批量填充空值

  • 输入文件:某区域的净初级生产力NPP栅格数据
  • 处理目标:将NPP的缺失值区域填充为0
Con(IsNull({A}), 0, {A})

在这里插入图片描述
在这里插入图片描述
还可以选择根据空白区域周围的有效像元插值进行填充,表达式可以参考

Con(IsNull({A}), FocalStatistics({A}, NbrRectangle(5,5, "CELL"), "MEAN"), {A})

关于上述表达式的说明,可以参考

缺失像元处的值就用以其为中心的5*5的矩形窗口内所有的有效像元的值平均后得到

FOCALMEAN中的长方形的尺寸可以被编辑以适应数据的需求。编辑"5,5"的值的大小可以改变长方形的尺寸。可能需要改成更大或者更小,这取决于数据。
原理是利用邻域均值代替空缺值。
可用的不同类型的邻域:
NbrAnnulus,NbrCircle,NbrRectangle,NbrWedge, NbrIrregular,和NbrWeight。默认邻域是正方形NbrRectangle,其宽度和高度为3个单元格。
要计算的统计类型:
MEAN(平均值)—计算附近单元格的平均值(平均值)。
MAJORITY —计算附近的单元格的多数(最常出现的值)。
MAXIMUM —计算附近的像元的最大值(最大值)。
MEDIAN —计算附近单元格的中位数。
MINIMUM —计算附近的像元的最小值(最小值)。
MINORITY-计算邻域中的单元格的少数(最少出现的值)。
RANGE(范围)—计算附近单元格的范围(最大值与最小值之间的差)。
STD —计算附近单元格的标准差。
SUM —计算附近单元格的总和(所有值的总和)。
VARIETY(变异度)—计算附近单元格的变异度(唯一值的数量)。
默认的统计类型是MEAN。
如果输入栅格为浮点,则只有MEAN,MAXIMUM,MINIMUM,RANGE,STD和SUM统计类型可用。

举个例子,对于上文使用的存在缺失区域的NPP栅格,使用相邻有效像元的值进行填充
在这里插入图片描述
在这里插入图片描述


  • 16
    点赞
  • 80
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 5
    评论
### 回答1: 批量裁剪栅格数据是指对多个栅格数据进行统一的裁剪操作,可以使用Python编写代码来实现该功能。 首先,需要导入必要的库,如`osgeo`库用于读取和处理栅格数据,`os`库用于文件操作。 ```python from osgeo import gdal, gdalnumeric import os ``` 接下来,可以定义一个函数来实现裁剪操作。 ```python def crop_raster(input_path, output_path, extent): # 打开输入栅格数据 ds = gdal.Open(input_path) if ds is None: print("打开栅格数据失败!") return # 读取输入栅格数据的信息 geotransform = ds.GetGeoTransform() projection = ds.GetProjection() x_size = ds.RasterXSize y_size = ds.RasterYSize # 计算裁剪后的输出栅格数据大小 ulx, uly = extent[0], extent[3] lrx, lry = extent[2], extent[1] px_start = int((ulx - geotransform[0]) / geotransform[1]) px_end = int((lrx - geotransform[0]) / geotransform[1]) py_start = int((uly - geotransform[3]) / geotransform[5]) py_end = int((lry - geotransform[3]) / geotransform[5]) px_count = px_end - px_start py_count = py_end - py_start # 创建输出栅格数据 driver = gdal.GetDriverByName("GTiff") output_ds = driver.Create(output_path, px_count, py_count, 1, gdal.GDT_Float32) # 设置输出栅格数据的信息 output_ds.SetProjection(projection) output_ds.SetGeoTransform((ulx, geotransform[1], 0, uly, 0, geotransform[5])) # 读取输入栅格数据的像素值,并写入输出栅格数据 input_data = ds.GetRasterBand(1).ReadAsArray(px_start, py_start, px_count, py_count) output_band = output_ds.GetRasterBand(1) output_band.WriteArray(input_data) # 关闭栅格数据集 ds = None output_ds = None ``` 在主程序中,可以通过遍历需要裁剪的栅格数据路径列表,调用`crop_raster`函数进行裁剪。 ```python if __name__ == "__main__": # 定义输入栅格数据路径列表 input_paths = ["input_1.tif", "input_2.tif", "input_3.tif"] # 定义输出栅格数据路径列表 output_paths = ["output_1.tif", "output_2.tif", "output_3.tif"] # 定义裁剪范围,格式为 [xmin, ymin, xmax, ymax] extent = [100, 50, 200, 150] # 遍历所有输入栅格数据 for i in range(len(input_paths)): input_path = input_paths[i] output_path = output_paths[i] # 调用裁剪函数进行裁剪 crop_raster(input_path, output_path, extent) print("裁剪完成!") ``` 以上代码片段可以实现对多个栅格数据进行批量裁剪,并将裁剪结果保存到指定的输出路径中。裁剪范围可以根据实际需求进行调整。 ### 回答2: 批量裁剪栅格数据是指对多个栅格数据进行裁剪操作,并将结果保存为新的栅格数据。以下是一个用 Python 编写的简单批量裁剪栅格数据的代码示例: ```python import arcpy # 设置工作空间 arcpy.env.workspace = "C:/RasterData" # 设置裁剪区域 clip_features = "C:/ClipData/clip.shp" # 获取所有栅格数据的文件名 raster_files = arcpy.ListRasters() # 遍历每个栅格数据进行裁剪 for raster_file in raster_files: # 构建裁剪后的栅格数据的输出路径 out_raster = "C:/OutputData/" + raster_file.split('.')[0] + "_clipped.tif" # 裁剪栅格数据 arcpy.Clip_management(raster_file, "#", out_raster, clip_features, "#", "ClippingGeometry") print("成功裁剪栅格数据:" + raster_file) print("批量裁剪栅格数据完成!") ``` 以上代码中,首先通过`arcpy.env.workspace`设置工作空间为包含待裁剪栅格数据的文件夹路径。然后使用`arcpy.ListRasters()`获取所有栅格数据的文件名。接下来通过遍历每个栅格数据文件,使用`arcpy.Clip_management()`函数进行裁剪操作,并将结果保存到指定文件夹。在裁剪完成后,程序会打印出成功裁剪的栅格数据文件名,并最终显示“批量裁剪栅格数据完成!”的提示。 ### 回答3: 批量裁剪栅格数据是一种对多个栅格数据进行裁剪操作的需求,可以使用Python实现。下面是一个大致的代码示例: ```python import os from osgeo import gdal from osgeo import ogr from osgeo import gdalconst # 设置要裁剪的区域范围 extent = [xmin, xmax, ymin, ymax] # 根据实际需要设置裁剪范围坐标 # 设置输入栅格数据的文件夹路径 input_folder = '/path/to/input/folder/' # 设置输出栅格数据的文件夹路径 output_folder = '/path/to/output/folder/' # 获取输入文件夹中的所有栅格数据文件 file_list = os.listdir(input_folder) tif_list = [file for file in file_list if file.endswith('.tif')] # 遍历每个栅格数据文件 for tif_file in tif_list: # 打开栅格数据 input_path = os.path.join(input_folder, tif_file) input_ds = gdal.Open(input_path, gdalconst.GA_ReadOnly) # 获取栅格数据的投影信息、地理变换等 prj = input_ds.GetProjection() geotransform = input_ds.GetGeoTransform() # 根据裁剪范围创建输出栅格数据的文件名和路径 output_path = os.path.join(output_folder, tif_file) # 创建输出栅格数据 output_ds = gdal.Warp(output_path, input_ds, format='GTiff', cutlineDSName=extent, cropToCutline=True) # 关闭栅格数据文件 input_ds = None output_ds = None ``` 以上代码使用了`gdal`库来进行栅格数据的读取和裁剪操作。首先,设置要裁剪的区域范围(extent)和输入、输出文件夹的路径。然后,遍历输入文件夹中的每个栅格数据文件,在每次循环中,打开栅格数据文件,获取其投影信息和地理变换参数,并根据裁剪范围创建输出栅格数据的文件名和路径。最后,使用`gdal.Warp`函数对栅格数据进行裁剪,并关闭输入和输出栅格数据文件。 注意:以上代码仅为大致示例,并未完全测试,实际使用中可能需要根据具体数据格式和需求进行适当调整和优化。

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 5
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

打赏作者

Salierib

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

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

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

打赏作者

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

抵扣说明:

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

余额充值