基于Arcpy的栅格数据批量处理＋简单说明

最近写毕业论文，把气象数据和卫星影像插值的插值，拼接的拼接，临到要计算了，问题成堆的出现，把人累的半死。
本文的目的是给出一些常用的批量处理栅格数据的arcoy代码，做出简单说明，并说一些踩过的坑。

1、代码

批量裁剪(clip)代码：

import arcpy

arcpy.env.workspace="工作路径"
rasters=arcpy.ListRasters("*","tif")
mask="相应文件"
#此处只能输入相应shp文件
#否侧可能有：The clip feature is outside the raster extent.

for raster in rasters:
   out="路径"+"test_clip"+raster
   print out+"begin project "
   #文件名称不可超过13个字符，且不可以数字为开头
   arcpy.Clip_management(raster, "482801.137013 4440341.783397 1262801.137013 5196351.760010", out, mask, "", "ClippingGeometry","MAINTAIN_EXTENT")

   #raster为输入数据，即被裁剪栅格图层，四个数字代表裁剪图层的最小外接矩形，
   #如果该参数为默认值，即
   #，表示不特殊制定边界，与裁剪的矢量范围保持一致，out是输出文件路径及名称，
   #mask是裁剪图层，“”代表该参数使用默认值，ClippingGeometry为按照几何边界裁剪。
   

批量掩膜代码：

# Import system modules

import arcpy

from arcpy import env

from arcpy.sa import *

 

# Set environment settings

env.workspace = "D:/data/"

rasterList = arcpy.ListRasters("*","tif")

 

#输出路径

output_path = "D:/8/"

# mask shp

inMaskData = "E:/haolu/XLGLM/XLGLM.shp"

 

for raster in rasterList:

    print raster

    # Set local variables

    inRaster = raster

    rd = arcpy.sa.Raster(inRaster)

    # Set the extent environment as the raster, very important for clip with different vector

    arcpy.env.extent = rd.extent

    # Execute ExtractByMask

    outExtractByMask = ExtractByMask(inRaster, inMaskData)

    # Save the output

    out = output_path + inRaster #对生成文件进行命名

    outExtractByMask.save(out)

批量重投影代码：

import arcpy,os,os.path
 

def projectRaster(rootPath):

    try:

        

        ##arcpy工作目录

        root_path = r"D:\data\LAIresult\clip"

        arcpy.env.workspace ="D:\data\LAIresult\clip"

 

        ##待处理文件所在目录(相对于根目录)

        input_path = "2005"

        output_path = "result"

 

        ##源坐标系 "CGCS2000_3_Degree_GK_CM_123E" 

        #sourceSR = arcpy.SpatialReference("WGS 1984 Albers")

        ##目标坐标系(WGS 1984 Web Mercator Auxiliary Sphere)

        targetSR = arcpy.SpatialReference("Asia North Albers Equal Area Conic")

 

        ##遍历目录，查找栅格数据

        files = os.listdir(root_path+os.sep+input_path)

        for f in files:

            if os.path.splitext(f)[1].upper() == ".TIF":

                fileName = os.path.splitext(f)[0] + ".tif"

                in_dataset = input_path + os.sep + fileName

                out_dataset = output_path + os.sep + fileName

 

                #print "begin project "+in_dataset+" from: " +sourceSR.name+" to: "+targetSR.name
                #arcpy.ProjectRaster_management(in_dataset, out_dataset,
                print "begin project "+in_dataset

                #arcpy.BatchProject_management(in_dataset, out_dataset,targetSR )

                arcpy.ProjectRaster_management(in_dataset, out_dataset, targetSR, "NEAREST",\

                                       "#", "#", "#")
    

 

        print "project success!"

        

    except arcpy.ExecuteError:

        print "Project Raster example failed."

        print arcpy.GetMessages()

 

################################################


if __name__ == '__main__':

 

    #指定处理文件根目录
    root_path = r"D:\data\ET\clip"

    projectRaster(root_path)

注意:批量重投影时，可以选择原有投影与目标投影，若不知道原有投影，直接选择目标投影也可。
投影的选择参考ArcGIS下自带投影文档。
批量重采样代码：

import arcpy
 
arcpy.env.workspace = r"D:/data/LAIresult/clip/2005/"
rasterList = arcpy.ListRasters("*","tif")

for raster in rasterList:
    print raster
    inRaster="D:/data/LAIresult/clip/2005/"+raster
    # Set local variables
    # Execute 
    out = "D:/data/LAIresult/clip/"+"resample"+raster
    arcpy.Resample_management(inRaster, out, "1000", "CUBIC")

以上就是我使用全部Arcpy批量处理代码，如果你并非地理信息专业，这些操作可以帮助你完成大部分的运算。
看到这里，相信大家也都能看出来，使用Arcpy进行栅格批量处理步骤是十分固定的，如果还想完成其他操作，自行去看相关官方文档即可。

坑

我在论文这次要做到的是，对两个栅格影像进行因果分析。因为因果分析代码是用matlab写的，而要把栅格导入matlab中，两幅影像完全对齐进行运算，是万分麻烦的事情，在过程中，一定要记住，如果计算完成后发现边界出错，请主要检查以下两点：1、投影 2、分辨率。
请牢记。
如果有空还会添加例子说明。