本文共 3330 字，大约阅读时间需要 11 分钟。

ArcGIS中创建Python脚本工具：计算多幅栅格影像的众数

起因

ArcMap自带了丰富的工具库，但在某些特殊需求下，自定义工具可能更适合项目需求。最近项目中需要计算多幅栅格影像的众数，现有工具中缺少这一功能，因此决定编写一个定制化的Python脚本。

脚本实现

为了实现计算多幅栅格影像的众数，主要使用了以下库：numpy、gdal和ArcGIS自带的arcpy。以下是代码实现：

# -*- coding: utf-8 -*-
import arcpy
import os
import numpy as np
import gdal
import sys
from gdalconst import *
def Search_File(base_directory, file_type, filelist=None):
    if filelist is None:
        filelist = []
    for file in os.listdir(base_directory):
        path = os.path.join(base_directory, file)
        if os.path.isdir(path):
            Search_File(path, file_type, filelist)
        else:
            if os.path.splitext(file)[1] == file_type:
                filelist.append(path)
    return filelist
def Readinfo(raster_file):
    ds = gdal.Open(raster_file, GA_ReadOnly)
    if ds is None:
        print('无法打开文件：{}'.format(raster_file))
        sys.exit(1)
    rows = ds.RasterYSize
    cols = ds.RasterXSize
    band = ds.GetRasterBand(1)
    no_data_value = band.GetNoDataValue()
    projection = ds.GetProjection()
    geotransform = ds.GetGeoTransform()
    return rows, cols, no_data_value, geotransform, projection
def WriteGTiffFile(filename, n_rows, n_cols, data, geotrans, proj, no_data_value, gdal_type):
    format = "GTiff"
    driver = gdal.GetDriverByName(format)
    ds = driver.Create(filename, n_cols, n_rows, 1, gdal_type)
    ds.SetGeoTransform(geotrans)
    ds.SetProjection(proj)
    ds.GetRasterBand(1).SetNoDataValue(no_data_value)
    ds.GetRasterBand(1).WriteArray(data)
    ds = None
def ras_to_array(ras_list):
    arcpy.AddMessage('将栅格转换为2D数组...')
    array_list = []
    for ras in ras_list:
        nodata = -1
        array = arcpy.RasterToNumPyArray(ras, '', '', '', nodata)
        array_list.append(array)
    num_array = np.array(array_list, dtype='int8')
    shape = num_array.shape
    res_array = np.full((shape[1], shape[2]), -1, dtype='int8')
    arcpy.AddMessage('计算大矩阵众数...')
    for i in range(shape[1]):
        for j in range(shape[2]):
            if np.min(num_array[:, i, j]) < 0:
                continue
            counts = np.bincount(num_array[:, i, j])
            res_array[i, j] = np.argmax(counts)
            print('计算结果：{}'.format(res_array[i, j]))
    return res_array, ras_list[0]
def array_to_ras(res_array, ras_info):
    arcpy.AddMessage('将数组转换为栅格...')
    my_raster = arcpy.Raster(ras_info)
    x_min = my_raster.extent.XMin
    y_min = my_raster.extent.YMin
    my_raster_result = arcpy.NumPyArrayToRaster(res_array, arcpy.Point(x_min, y_min),
                                                  my_raster.meanCellWidth, my_raster.meanCellHeight, -1)
    arcpy.AddMessage('转换完成!')
    return my_raster_result
def Mode_cal(ras_array, ras_info, out_name):
    result = array_to_ras(ras_array, ras_info)
    result.save(out_name)
    arcpy.AddMessage('输出文件：{}'.format(out_name))
arcpy.env.overwriteOutput = True
# 工具参数
raster_workspace = arcpy.GetParameterAsText(0)
out_name = arcpy.GetParameterAsText(1)
raslist = Search_File(raster_workspace, '.tif')
res = ras_to_array(raslist)
gdal_array_to_ras(res[0], res[1], out_name)

工具脚本创建

代码实现完成后，可以按照以下步骤创建ArcGIS脚本工具：

通过以上步骤，可以自定义ArcGIS工具，满足项目需求。如有疑问或优化建议，欢迎在评论区留言！

转载地址：http://vxka.baihongyu.com/