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

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

起因

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

脚本实现

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

# -*- coding: utf-8 -*-import arcpyimport osimport numpy as npimport gdalimport sysfrom 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 filelistdef 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, projectiondef 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 = Nonedef 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_resultdef 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/