当前位置：首页 > news >正文

【Arcpy】入门学习笔记（四）-栅格数据

news 2026/5/12 7:28:53

前言
栅格数据处理
- 地形因子提取
- 影像数据处理
- 计算植被覆盖度
一些注意点
- 数据量的问题
- 栅格计算器
总结
参考

前言

栅格数据和矢量数据相关的知识点比较多，因此建议根据网络上的案例边练边学。本次列举了三个简单的例子，分别是地形因子的提取，影像数据处理，计算植被覆盖度。

官方说明：

栅格数据集通过将世界分割成在格网上布局的离散方形或矩形像元来表示地理要素。每个像元都具有一个值，用于表示该位置的某个特征，例如温度、高程或光谱值。

栅格数据集常用于表示和管理影像、数字高程模型及许多其他现象。通常，栅格是用于表示点、线和面要素的一种方法。在下面的示例中，可以看到如何将一系列面表示为栅格数据集。

栅格可用于表示所有地理信息（要素、影像和表面），并且具有丰富的分析地理处理运算符。除了在 GIS 中作为保存影像的通用数据类型外，栅格还经常用于表示要素，从而允许在基于栅格的建模和分析中使用所有地理对象。

栅格数据处理

栅格数据处理有很多方法，推荐找几个例子练习一下就容易记住，下面是我随便找的几个数据做的练习，因为涉及到数据处理，所以代码的很大部分也是需要熟悉Python内部的一些文件操作方法和函数

地形因子提取

这个代码提取了六种地形因子，其中山体阴影需要根据实际情况进行修改；主要涉及到表面分析工具的关于DEM的方法

原始数据：

输出数据：

代码：

# -*- coding:utf-8 -*-# -------------------------------------------------------------------------------
# Description  : 
# @Time        : 2024-07-24 9:32
# @Author      : zbh
# -------------------------------------------------------------------------------import arcpy
from arcpy.sa import *
import os# 设置环境
workspace = r"./data"
arcpy.CheckOutExtension("spatial")
arcpy.env.parallelProcessingFactor = '0'
arcpy.env.overwriteOutput = "True"
arcpy.env.workspace = workspacedef calculate_all_factors(input_raster, output_folder):# 获取dem数据的文件信息raster = arcpy.Raster(input_raster)raster_name = os.path.basename(input_raster)base_name, ext = os.path.splitext(raster_name)# 计算坡度output_slope = Slope(raster, "Degree")output_slope_path = os.path.join(output_folder, base_name + "_" + "slope" + ext)output_slope.save(output_slope_path)# 计算sosoutput_sos = Slope(output_slope, "Degree")output_sos_path = os.path.join(output_folder, base_name + "_" + "sos" + ext)output_sos.save(output_sos_path)del output_slopedel output_sos# 计算坡度output_aspect = Aspect(raster)output_aspect_path = os.path.join(output_folder, base_name + "_" + "aspect" + ext)output_aspect.save(output_aspect_path)# 计算soaoutput_soa = Slope(output_aspect, "Degree")output_soa_path = os.path.join(output_folder, base_name + "_" + "soa" + ext)output_soa.save(output_soa_path)del output_aspectdel output_soa# 计算曲率output_cur = Curvature(raster)output_cur_path = os.path.join(output_folder, base_name + "_" + "cur" + ext)output_cur.save(output_cur_path)del output_cur# 计算山体阴影output_shade = Hillshade(raster, 90, 45, "SHADOWS")output_shade_path = os.path.join(output_folder, base_name + "_" + "shade" + ext)output_shade.save(output_shade_path)del output_shadeprint("calculate_all_factors() success!")if __name__ == '__main__':calculate_all_factors(r'.\data\dem_sanqu.tif', r'.\data\output')print("ok")

影像数据处理

以一幅Landsat 8-9的 C2L2级别数据为例，批量裁剪，放缩与偏移，波段合成；主要涉及到os模块、正则表达式与Arcpy批量处理的结合

原始数据：

输出数据：

代码：

# -*- coding:utf-8 -*-# -------------------------------------------------------------------------------
# Description  : python 2.7, process Landsat8-9 C2L2 data batch
# @Time        : 2024-07-23 22:36
# @Author      : zbh
# -------------------------------------------------------------------------------import arcpy
from arcpy.sa import *
import os
import re
import fnmatchworkspace = r".\data"
arcpy.CheckOutExtension("spatial")
arcpy.env.parallelProcessingFactor = '0'
arcpy.env.overwriteOutput = "True"
arcpy.env.workspace = workspace# 批量裁剪
def clip_batch(input_folder, clip_feature, output_folder=None):# 输出文件路径不存在if not os.path.exists(output_folder):try:os.makedirs(output_folder)except OSError as e:raise OSError("Error creating output folder: {}".format(e))else:passraster_list = []# 从文件夹中获取输入数据, 并添加至列表files = os.listdir(input_folder)tif_files = [f for f in files if fnmatch.fnmatch(f, '*B*.TIF')]for tif_file in tif_files:full_path = os.path.join(input_folder, tif_file)if arcpy.Exists(full_path):raster_list.append(full_path)# 栅格数据批量裁剪for raster in raster_list:filename = os.path.basename(raster)base_name, ext = os.path.splitext(filename)last_underscore_index = base_name.rfind('_')if last_underscore_index != -1:extracted_name = base_name[last_underscore_index + 1:]else:extracted_name = base_nameout_clip_path = os.path.join(output_folder, extracted_name + "_clip" + ext)arcpy.management.Clip(raster, "#", out_clip_path, clip_feature, 0, "ClippingGeometry")print("clip_batch() sucess!")# 批量放缩和偏移
def process_batch(input_folder, is_composite=None, output_folder=None):# 输出文件路径不存在if not os.path.exists(output_folder):try:os.makedirs(output_folder)except OSError as e:raise OSError("Error creating output folder: {}".format(e))else:passraster_list = []# 从文件夹中获取输入数据files = os.listdir(input_folder)pattern = r'B(?:[1-9]|10)_clip+\.TIF$'tif_files = [f for f in files if re.match(pattern, f)]for tif_file in tif_files:full_path = os.path.join(input_folder, tif_file)if arcpy.Exists(full_path):raster_list.append(full_path)process_raster_str = ""for raster in raster_list:filename = os.path.basename(raster)base_name, ext = os.path.splitext(filename)out_process_path = os.path.join(output_folder, base_name + "_process" + ext)if "B10" in filename:raster = (arcpy.Raster(raster) * 0.00341802) + 149 - 273.15process_raster = Con(raster < -123, -123, Con(raster > 86, 86, raster))else:raster = (arcpy.Raster(raster) * 0.0000275) - 0.2process_raster = Con(raster < 0, 0, Con(raster > 1, 1, raster))process_raster_str += (out_process_path + ";")process_raster.save(out_process_path)if is_composite == "Composite":arcpy.CompositeBands_management(process_raster_str, os.path.join(output_folder, "L8" + ext))else:passprint("process_batch() success!")if __name__ == '__main__':# 要搜索的路径clip_input_folder = r"G:\WATER_TEST\LC08_L2SP_119042_20220712_20220722_02_T1"# 裁剪要素clip_feature = r".\data\fz_sub\fz_sub.shp"# 构建 clip 文件夹的路径clip_output_folder = os.path.join(r".\data\clip")# 裁剪clip_batch(clip_input_folder, clip_feature, clip_output_folder)# 要搜索的路径process_input_folder = r".\data\clip"# 构建 process 文件夹的路径process_output_folder = os.path.join(r".\data\process")# 像元值转换与波段合成process_batch(process_input_folder, "Composite", process_output_folder)print("OK")

计算植被覆盖度

使用处理好的Red和NIR波段，计算NDVI，之后再计算FVC。主要涉及到RasterToNumPyArray和RasterCalculator_sa的相关操作

原始数据：

输出数据：

代码：

# -*- coding:utf-8 -*-# -------------------------------------------------------------------------------
# Description  : python 2.7, calculate fvc
# @Time        : 2024-07-24 0:42
# @Author      : zbh
# -------------------------------------------------------------------------------import arcpy
from arcpy.sa import *
import os
import numpy as npworkspace = r".\data"
arcpy.CheckOutExtension("spatial")
arcpy.env.parallelProcessingFactor = '0'
arcpy.env.overwriteOutput = "True"
arcpy.env.workspace = workspace# 计算ndvi和植被覆盖度
def ndvi_fvc(red, nir, left_percent, right_percent, output_folder):# 输出文件夹不存在if not os.path.exists(output_folder):try:os.makedirs(output_folder)except OSError as e:raise OSError("Error creating output folder: {}".format(e))else:pass# 计算ndvindvi_raster_path = os.path.join(output_folder, "ndvi.tif")red_band = arcpy.Raster(red)nir_band = arcpy.Raster(nir)ndvi_raster = Float(arcpy.sa.Minus(nir_band, red_band)) / Float(arcpy.sa.Plus(nir_band, red_band))# 异常值处理ndvi_bounded = Con(ndvi_raster < -1, -1, Con(ndvi_raster > 1, 1, ndvi_raster))ndvi_bounded.save(ndvi_raster_path)# 计算置信区间值arr = arcpy.RasterToNumPyArray(ndvi_raster)nodata_value = ndvi_raster.noDataValuevalid_values = arr[arr != nodata_value]flat_arr = np.sort(valid_values)left_index = int(len(flat_arr) * float(left_percent))right_index = int(len(flat_arr) * float(right_percent))left_value = flat_arr[left_index]right_value = flat_arr[right_index]# 计算fvcvalues_range = float(right_value) - float(left_value)condition_fvc = r"""Con(("{0}"<={1}),0,Con(("{0}">={2}),1,("{0}"-{1})/{3}))""".format(ndvi_raster, left_value,right_value, values_range)fvc_raster_path = os.path.join(output_folder, "fvc.tif")arcpy.gp.RasterCalculator_sa(condition_fvc, fvc_raster_path)print(ndvi_raster_path, left_value, right_value, fvc_raster_path)print("ndvi_fvc() success!")if __name__ == "__main__":# red和nir波段red = r".\data\process\B4_clip_process.TIF"nir = r".\data\process\B5_clip_process.TIF"# 构建输出文件夹的路径ndvi_output_folder = os.path.join(r"./data/veg")ndvi_fvc(red, nir, 0.05, 0.95, ndvi_output_folder)print("OK")