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

基于Arcpy和MATLAB批量提取指定经纬度点的栅格数据并转换为矩阵格式

article 2024/12/24 11:39:00

前言

Python代码（ArcPy）

这段Python代码使用ArcPy模块进行GIS操作，具体来说是从指定目录下的多个栅格文件中提取点处的栅格值。

MATLAB代码

这段MATLAB代码从指定文件夹中读取shapefile文件，提取每个文件中名为RASTERVALU的列，并将数据保存为矩阵。
我们也可以直接逐一提取对应点的值，但是那样会比较慢，上面的方法会很快，而且省事。

简述

代码各部分的功能

导入模块：从ArcPy中导入必要的模块，用于GIS处理和空间分析。
设置环境工作空间：指定文件夹路径，其中包含所有的.tif栅格文件。
列出栅格文件：arcpy.ListRasters("*", "tif")获取指定文件夹中的所有.tif栅格文件。
定义点要素文件：inPointFeatures指定输入点要素文件的路径，该文件包含提取栅格值的位置。
遍历栅格文件：对rasterlist中的每个栅格文件进行处理：
    通过栅格文件名定义输出的点要素文件路径。
    检出Spatial Analyst扩展许可证以启用空间分析功能。
    使用ExtractValuesToPoints函数提取每个栅格在指定点位置的值。
    打印消息表示该栅格文件的提取已经完成。

整体目的

这段代码遍历指定目录下的所有.tif栅格文件，从Newpoints.shp文件中的点位置提取栅格值，并将每个输出保存为一个单独的shapefile。

二、代码实操

1.python

# 导入系统模块
import arcpy
from arcpy import env
from arcpy.sa import *

# 设置环境工作空间
arcpy.env.workspace = "E:/LHASADATA/"
# 列出工作空间中的所有.tif栅格文件
rasterlist = arcpy.ListRasters("*", "tif")

# 设置输入点要素
inPointFeatures = "E:/LHASADATA//Newpoints.shp"

# 遍历每一个栅格文件
for raster in rasterlist:
    # 定义每个栅格的输出点要素路径，使用栅格文件名生成对应的输出文件名
    outPointFeatures = "E:/LHASADATA/Dryandwet/" + raster + ".shp"  

    # 检出ArcGIS的Spatial Analyst扩展模块许可证
    arcpy.CheckOutExtension("Spatial")

    # 执行ExtractValuesToPoints函数，将栅格值提取到指定点
    ExtractValuesToPoints(inPointFeatures, raster, outPointFeatures, "NONE", "VALUE_ONLY")
    
    # 打印消息，表示该栅格的提取完成
    print(raster + " has done")

2.matlab

% 清空命令窗口和工作空间
clc
clear all

% 指定包含shapefile文件的文件夹路径
folderPath = 'E:\LHA\';

% 获取文件夹中所有shapefile文件的列表
shapefiles = dir(fullfile(folderPath, '*.shp'));

% 初始化一个单元格数组，用于存储提取的'RASTERVALU'数据
rasterValuesCell = cell(length(shapefiles), 1);

% 遍历文件夹中的每个shapefile文件
for k = 1:length(shapefiles)
    % 获取当前shapefile文件的完整路径
    shapefilePath = fullfile(folderPath, shapefiles(k).name);
    
    % 读取shapefile数据
    shapefileData = shaperead(shapefilePath);
    
    % 提取'RASTERVALU'字段（如果存在）
    if isfield(shapefileData, 'RASTERVALU')
        % 将'RASTERVALU'列数据提取为数值数组
        rasterValues = [shapefileData.RASTERVALU];
    else
        % 如果'RASTERVALU'字段不存在，则设置为空数组
        rasterValues = [];
        warning('RASTERVALU field not found in %s', shapefiles(k).name);
    end
    
    % 将提取的值存储在单元格数组中
    rasterValuesCell{k} = rasterValues;
end

% 确定所有shapefile文件中'RASTERVALU'数据的最大长度
maxLength = max(cellfun(@length, rasterValuesCell));

% 初始化矩阵并填充NaN值，较短的行将用NaN填充
rasterValuesMatrix = NaN(length(shapefiles), maxLength);

% 将每个shapefile文件的'RASTERVALU'数据填充到矩阵中
for k = 1:length(shapefiles)
    numValues = length(rasterValuesCell{k});
    rasterValuesMatrix(k, 1:numValues) = rasterValuesCell{k};
end
% 显示完成消息
disp('done')