基于栅格生成随机采样点_arcgis在有值的栅格数据上生成随机点-优快云博客

本文链接：https://blog.youkuaiyun.com/dwcai_mind/article/details/149066172

基于栅格生成随机采样点

1. 功能需求
2. 使用方法
3. 代码
4. 工具截图
5. 示例
6. 链接

1. 功能需求

在ArcGIS里，目前只有基于矢量数据的采样工具“创建随机点”，很多时候我们需要直接从栅格上面创建随机点，因此，我结合莫顿编码（Morton）规则，在ArcGIS Pro环境中，借助ArcPy站点包和numpy包，构建了“莫顿编码栅格随机采样”工具。以下代码可在ArcGIS Pro中作为脚本工具使用。。

2. 使用方法

将下面的代码保存为Python脚本工具箱.pyt 文件（例如 RasterSamplingTools.pyt）
在ArcGIS Pro中，右键点击工具箱目录，选择"Add Toolbox"
导航到保存的.pyt文件并添加该.pyt文件即可使用

3. 代码

# -*- coding: utf-8 -*-
# Copyright: @dwcai 
# 20250702

import arcpy
import numpy as np
from itertools import zip_longest

class Toolbox(object):
    def __init__(self):
        """定义工具箱（工具箱的名称是. pyt文件的名称）。"""
        self.label = "栅格随机采样工具"
        self.alias = "RasterSamplingTools"

        # 列出此工具箱中包含的所有工具类
        self.tools = [MortonRandomSampling]

class MortonRandomSampling(object):
    def __init__(self):
        """定义工具(工具名称是类的名称)"""
        self.label = "莫顿编码栅格随机抽样"
        self.description = "使用莫顿编码对栅格数据进行随机抽样，并输出采样点"
        self.canRunInBackground = True

    def getParameterInfo(self):
        """定义工具参数"""
        # 输入栅格参数
        param0 = arcpy.Parameter(
            displayName="输入栅格",
            name="in_raster",
            datatype="GPRasterLayer",
            parameterType="Required",
            direction="Input")
        
        # 抽样方式参数
        param1 = arcpy.Parameter(
            displayName="抽样方式",
            name="sampling_method",
            datatype="GPString",
            parameterType="Required",
            direction="Input")
        param1.filter.type = "ValueList"
        param1.filter.list = ["按百分比抽样", "按固定数量抽样"]
        param1.value = "按百分比抽样"
        
        # 抽样值参数
        param2 = arcpy.Parameter(
            displayName="抽样值",
            name="sample_value",
            datatype="GPDouble",
            parameterType="Required",
            direction="Input")
        param2.value = 10
        
        # 输出点要素参数
        param3 = arcpy.Parameter(
            displayName="输出点要素",
            name="out_points",
            datatype="DEFeatureClass",
            parameterType="Required",
            direction="Output")
        
        params = [param0, param1, param2, param3]
        return params

    def isLicensed(self):
        """设置工具是否可以执行"""
        return True

    def updateParameters(self, parameters):
        """根据用户输入修改参数属性"""
        # 当抽样方式改变时，调整抽样值的范围
        if parameters[1].valueAsText == "按百分比抽样":
            # parameters[2].value = 10  # 默认10%
            parameters[2].value = parameters[2].value  # 默认10%
            parameters[2].filter.type = "Range"
            parameters[2].filter.list = [0.01, 100]  # 0.01%到100%
        else:
            # parameters[2].value = 1000  # 默认1000个点
            parameters[2].value = parameters[2].value  # 默认1000个点
            parameters[2].filter.type = "Range"
            parameters[2].filter.list = [1, 1e8]  # 1到1亿个点
        
        return

    def updateMessages(self, parameters):
        """修改由内部验证为每个工具参数创建的消息"""
        return

    def mortonEncode(self, argy, argx):
        """莫顿编码"""
        biny, binx = np.binary_repr(argy, 32), np.binary_repr(argx, 32)
        binMD = ''.join(i for yx in zip_longest(biny, binx, fillvalue='0') for i in yx)
        decMD = int(binMD, 2)
        return decMD

    def mortonDecode(self, decMD):
        """莫顿解码"""
        binMD = np.binary_repr(decMD, 64)
        n = len(binMD)
        biny = ''.join(binMD[i] for i in np.arange(0, n, 2))
        binx = ''.join(binMD[i] for i in np.arange(1, n, 2))
        argy, argx = int(biny, 2), int(binx, 2)
        return argy, argx
    
    def pixelTypeNodata(self, pixelType):
        PIXEL_TYPE = {
            "U1": 0, "U2": 0, "U4": 0, "U8": 0, "S8": -128, "U16": 65535, "S16": -32768, 
            "U32": 4294967295, "S32": -2147483648, "F32": np.nan, "F64": np.nan,  
        }
        return PIXEL_TYPE[pixelType]

    def execute(self, parameters, messages):
        """执行工具"""
        arcpy.env.overwriteOutput = True
        
        # 获取参数
        in_raster = parameters[0].valueAsText
        sampling_method = parameters[1].valueAsText
        sample_value = parameters[2].value
        out_points = parameters[3].valueAsText
        
        # 确定抽样方式
        if sampling_method == "按百分比抽样":
            sample_pct = sample_value
            sample_size = None
        else:
            sample_pct = None
            sample_size = int(sample_value)
        
        # 开始处理
        arcpy.AddMessage("开始处理栅格随机抽样...")
        
        # 获取栅格信息
        ras = arcpy.sa.Raster(in_raster)
        ny, nx = ras.height, ras.width
        pixel_type = ras.pixelType
        # nodata = ras.noDataValue  # 有些栅格没有给定无效值
        nodata = self.pixelTypeNodata(pixel_type)
        
        # 读取栅格数据
        arcpy.AddMessage("正在读取栅格数据...")
        data = arcpy.RasterToNumPyArray(ras, nodata_to_value=nodata)
        
        # 生成莫顿矩阵
        arcpy.AddMessage("正在生成莫顿编码矩阵...")
        morton_mtx = np.zeros([ny, nx], dtype=np.uint64)
        for i in range(ny):
            for j in range(nx):
                if not (np.isnan(data[i, j]) | (data[i,j]==nodata)):
                    morton_mtx[i, j] = self.mortonEncode(i, j)
        
        # 提取有效莫顿编码
        not_valid_mask = (np.isnan(data) | (data==nodata))
        morton_scalar = morton_mtx[~not_valid_mask]
        
        # 确定样本大小
        if sample_size is None:
            sample_size = min(int(len(morton_scalar) * sample_pct / 100), len(morton_scalar))
            arcpy.AddMessage(f"将抽取 {sample_size} 个样本点 (占总有效像素的 {sample_pct}%)")
        else:
            sample_size = min(sample_size, len(morton_scalar))
            arcpy.AddMessage(f"将抽取 {sample_size} 个样本点")
        
        # 随机抽样
        arcpy.AddMessage("正在进行随机抽样...")
        morton_samples = np.random.choice(morton_scalar, size=sample_size, replace=False)
        
        # 创建结果矩阵
        arcpy.AddMessage("正在创建结果矩阵...")
        result_mtx = np.full(data.shape, nodata, dtype=data.dtype)
        
        # 解码并填充结果矩阵
        arcpy.AddMessage("正在解码莫顿编码并生成结果...")
        for code in morton_samples:
            row, col = self.mortonDecode(code)
            result_mtx[row, col] = data[row, col]
        
        # 将结果写入临时栅格
        arcpy.AddMessage("正在创建临时栅格...")
        temp_raster_path = arcpy.env.scratchGDB + "\\temp_sample_raster"
        arcpy.NumPyArrayToRaster(
            result_mtx, 
            arcpy.Point(ras.extent.XMin, ras.extent.YMin), 
            ras.meanCellWidth, 
            ras.meanCellHeight, 
            nodata
        ).save(temp_raster_path)
        # 后处理部分
        arcpy.AddMessage("定义临时栅格的投影信息...")
        arcpy.management.DefineProjection(temp_raster_path, ras.spatialReference)
        # 栅格转点
        arcpy.AddMessage("正在生成输出点要素...")
        arcpy.conversion.RasterToPoint(temp_raster_path, out_points, "VALUE")
        # 清理临时数据
        arcpy.management.Delete(temp_raster_path)
        arcpy.AddMessage("处理完成！")
        return