起因

ArcMap中本身就已集成了非常丰富的工具库，但有时候由于一些特殊的使用需求，我们需要构建更适用于当前项目的工具。最近由于项目中需要计算多幅栅格影像的众数，查了查已有的工具发现栅格计算中发现只有较为常见的最大值、最小值、平均值等函数，于是决定自己写个脚本来做个定制化工具。

脚本实现

求取多幅栅格影像的众数的代码如下，主要用到了numpy、gdal和arcgis自带的arcpy这三个库：

# -*- coding: utf-8 -*-import arcpyimport osimport numpy as npimport gdalimport sysfrom gdalconst import *#Create a tool box for calculating mode number of rster datasets#查询当前文件夹下所有同类型后缀名的文件def Search_File(BaseDirectory,fileType,filelist=None):   if filelist is None:       filelist=[]   file=os.listdir(BaseDirectory)   for i in file:      x=BaseDirectory+'\\'+i      if os.path.isdir(x):         Search_File(x,fileType,filelist)      elif os.path.isfile(x):         if os.path.splitext(x)[1]==fileType:            filelist.append(x)   return filelist#gdal库读取栅格信息def Readinfo(RasterFile):    ds = gdal.Open(RasterFile,GA_ReadOnly)    if ds is None:        print 'Cannot open ',RasterFile        sys.exit(1)    cols = ds.RasterXSize    rows = ds.RasterYSize    band = ds.GetRasterBand(1)    noDataValue = band.GetNoDataValue()    projection=ds.GetProjection()    geotransform = ds.GetGeoTransform()    return rows,cols,noDataValue,geotransform,projection#gdal库保存栅格影像def WriteGTiffFile(filename, nRows, nCols, data,geotrans,proj, noDataValue, gdalType):    format = "GTiff"    driver = gdal.GetDriverByName(format)    ds = driver.Create(filename, nCols, nRows, 1, gdalType)    ds.SetGeoTransform(geotrans)    ds.SetProjection(proj)    ds.GetRasterBand(1).SetNoDataValue(noDataValue)    ds.GetRasterBand(1).WriteArray(data)    ds = None#Transform raster to 2d array，栅格转数组def ras_to_array(ras_list):    arcpy.AddMessage(u"Transforming raster to 2d array......")    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(u"Calculating ModeNumber of big matrix......")    #calculate every single array，计算众数    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 res_array[i, j]    # np.save('F:\Arcgis_process\Test\modenumber.npy',res_array)    # np.savetxt('F:\Arcgis_process\Test\modenumber.txt',res_array)    arcpy.AddMessage(u"Value Calculated!")    return res_array,ras_list[0]#numpy to raster，数组转栅格，arcpy内置方法，对行列数较多的栅格影像来说，容易报内存错误def array_to_ras(res_array,ras_info):    arcpy.AddMessage(u"Converting array to raster......")    myRaster = arcpy.Raster(ras_info)    #lower left coordinate    x=myRaster.extent.XMin    y=myRaster.extent.YMin    print myRaster.meanCellWidth,myRaster.meanCellHeight    myRaster_result=arcpy.NumPyArrayToRaster(res_array,arcpy.Point(x,y),myRaster.meanCellWidth,myRaster.meanCellHeight,-1)    arcpy.AddMessage(u"Converted over!")    return myRaster_result#用gdal库中的方法实现数组转栅格def gdal_array_to_ras(res_array,ras_info,out_name):    rows, cols, noDataValue, geotransform, projection=Readinfo(ras_info)    WriteGTiffFile(out_name, rows, cols, res_array, geotransform, projection, -1, GDT_Int16)#save result结果保存def Mode_cal(ras_array,ras_info,out_name):    result=array_to_ras(ras_array, ras_info)    result.save(out_name)    arcpy.AddMessage(u"Exported the image!")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官方帮助文档：，创建脚本工具中需要注意的主要有以下三点：

python脚本中尽量不要出现中文字符，不论是注释还是文件名，否则很容易报错，血泪的教训；

有可能的话不要使用自建的包，最好完全使用arcpy中自带的类和函数，不然会报模块不存在的错误（这点我不清楚为什么，希望有知道的朋友能解答一下）

工具界面中输入项的direction属性是input，输出项的属性是output，千万不能忘记选择，并且还要明确自己输入输出数据的类型，是文件还是文件夹。

以上就是今天实现利用python创建脚本工具的过程与心得，与大家共勉！

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

你可能感兴趣的文章

MySQL 证明为什么用limit时，offset很大会影响性能