Arcgis-Tools_02影像处理

前言

介绍影像处理中常用的一些批量处理方法。

创建金字塔

批量创建金字塔并无技巧可言，只需遍历到每一个影像文件，调用构建金字塔工具。

示例代码

import arcpy, os

in_folder = r"C:\Users\Admin\Desktop\Raster"
walk = arcpy.da.Walk(in_folder, datatype="RasterDataset")
for dirpath, dirnames, filenames in walk:
    for filename in filenames:
        in_raster = os.path.join(dirpath + filename)
        print(in_raster)
        arcpy.BuildPyramids_management(in_raster)

格式转换

经常获取到高度压缩的 ers 格式影像（只知为ER Mapper产出），用 ArcMap 无法为其定义投影，就需要先将其转换为 tif 格式。

可使用栅格转其他格式(批量)工具来批量转换。此工具位置：转换工具——转为栅格——栅格转其他格式(批量)。

定义投影

在遇到影像未进行投影，而数量很多，采用的投影又多样的情况下，可以根据影像 X 坐标形式来反推坐标。

注意：至少要知道影像的地理坐标，如不能把 XIAN80 定义成 CGCS2000 或者 WGS84。

import arcpy, os

in_folder = r"C:\Users\Admin\Desktop\Raster"
sr_folder = os.getcwd() + os.sep + "CGCS2000"
for root, dirs, files in os.walk(in_folder):
    for file in files:
        if file[-4:] == ".tif":
            print file
            in_raster = root + os.sep + file
            left = float(arcpy.GetRasterProperties_management(in_raster, "LEFT")[0])
            x_length = len(str(int(left)))
            sr_name = ""
            if x_length == 6:
                sr_name = "CGCS2000 GK CM " + str(6 * int(file[-10:-8]) - 3) + "E.prj"
            elif x_length == 8:
                sr_name = "CGCS2000 GK Zone " + str(int(file[-10:-8])) + ".prj"
            elif x_length == 2 or x_length == 3:
                sr_name = "China Geodetic Coordinate System 2000.prj"
            print sr_name
            sr = arcpy.SpatialReference(sr_folder + os.sep + sr_name)
            arcpy.DefineProjection_management(in_raster, sr)

注：

• 影像名称为标准名称 zy3_030430_015127_20170702_19_fus.tif ，可从中提取到带号。

• 影像幅面较大，这里判断带号后采用六度带。，

• 脚本路径下面有文件夹存储 CGCS2000 的 prj文件。

获取影像边界

在工作中经常会用到影像的矢量边界，其中存储时间、像元大小、波段等信息，时间有时候我们可以根据影像的文件名称来判断，像元大小、波段等影像信息可以通过 Describe 函数 获取栅格属性。对于矢量边界，除了人工手画之外，可通过多个Arcgis工具组合以达到提取的目的。

获取流程：

创建地理数据库——创建镶嵌数据集——添加栅格至镶嵌数据集——构建轮廓——导出镶嵌数据集几何

针对以上流程，可手动按步骤处理或者使用模型构建器创建工具达到重复使用的目的。

也可使用以下Python代码

# -*- coding: utf-8 -*-
import arcpy, os

# in_folder = arcpy.GetParameterAsText(1)
# spatialReference = arcpy.GetParameter(2)
# out_rasterBoundary = arcpy.GetParameterAsText(3)

in_folder = r"C:\Users\Admin\Desktop\tif"
spatialReference = arcpy.SpatialReference(4490)
out_rasterBoundary = os.getcwd() + os.sep + "RasterBoundary.shp"
#   设置并行处理为1
arcpy.env.parallelProcessingFactor = 1
try:
    arcpy.CreateFileGDB_management(os.getcwd(), "Temp.gdb")
    arcpy.CreateMosaicDataset_management(os.getcwd() + os.sep + "Temp.gdb", "Mosaic", spatialReference)
    in_mosaic = os.getcwd() + os.sep + "Temp.gdb" + os.sep + "Mosaic"
    # Process: 添加栅格至镶嵌数据集
    arcpy.AddRastersToMosaicDataset_management(in_mosaic, "Raster Dataset", in_folder, "UPDATE_CELL_SIZES")
    # Process: 构建轮廓
    arcpy.BuildFootprints_management(in_mosaic)
    # Process: 导出镶嵌数据集几何
    arcpy.ExportMosaicDatasetGeometry_management(in_mosaic, out_rasterBoundary, "", "FOOTPRINT")
    arcpy.Delete_management(os.getcwd() + os.sep + "Temp.gdb")
    print("Completed")
    arcpy.AddMessage("Completed")
except:
    for msg in range(0, arcpy.GetMessageCount()):
        if arcpy.GetSeverity(msg) == 2:
            print(msg)
            arcpy.AddReturnMessage(msg)

裁剪影像

裁剪影像首先要明确两个对象，用什么去裁剪，裁剪什么。一对一即单个矢量图形对应单个影像文件，多对一即多个矢量对应一个影像文件。

一对一裁剪

要明确裁剪的矢量图形记录，对应的影像文件路径及名称，输出的影像文件路径及名称。这些都可以存储在矢量文件的属性表中。

Python代码示例：

# -*- coding: cp936 -*-
# 脚本适合一对一裁剪，即一条矢量记录对应一块相同名称的影像，矢量可以为多部件或者孔洞多边形
# 当矢量与影像无重叠时，会导出对应名称的影像全部
# 影像和矢量的单位必须一致，同为经纬度或者米
# 当输入或者输出格式不为tif时，需要更改输入对应影像名

import arcpy, os

# 输入影像文件夹
in_folder = r"F:\MyRaster\陕西\ecw".decode('gb2312')
# 裁剪框
in_polygon = r"F:\MyRaster\陕西\shp\ClipShp.shp".decode('gb2312')
# 导出路径
out_folder = r"F:\陕西分县".decode('gb2312')
# 影像名称字段
field = "RasterName"
# 区域字段
zone_field = "XIAN_NAME"

#   设置并行处理为1
arcpy.env.parallelProcessingFactor = 1

try:
    # 创建一个游标
    cursor = arcpy.SearchCursor(in_polygon)
    # 获取裁剪框总数
    total_count = arcpy.GetCount_management(in_polygon)
    i = 1
    for row in cursor:
        raster_name = row.getValue(field)
        # 输出当前分割的影像名称
        print (raster_name)
        # 当前分割的影像的路径
        raster = in_folder + os.sep + raster_name + ".ers"
        print raster
        # 判断影像是否存在
        if arcpy.Exists(raster):
            # 获取裁剪框的范围
            mask = row.getValue("Shape")
            extent = str(mask.extent)
            # 输出路径
            out_path = out_folder + os.sep + row.getValue(zone_field) + "2019" + ".tif"
            print out_path
            # 调用裁剪工具进行裁剪
            arcpy.Clip_management(raster, extent, out_path, mask, "0", "ClippingGeometry")
            print "finished " + str(i) + "  of " + str(total_count) + " " + row.getValue(zone_field)
            i = i + 1
        else:
            print raster_name + " not exist"
except:
    for msg in range(0, arcpy.GetMessageCount()):
        if arcpy.GetSeverity(msg) == 2:
            print(msg)
            arcpy.AddReturnMessage(msg)

多对一裁剪

这种情况只是一对一裁剪的特殊情况，只不过把对应的影像文件路径及名称换成了单一值，产出结果也无需拼接。

拼接影像

拼接影像要明确的是要将哪些影像拼接、输出影像的名称及路径。

注意：拼接的影像要根据像元大小保留一定重叠比例（裁剪或者数据准备阶段需注意），防止中间出现缝隙或者重叠过大。

Python代码示例（代码为拼接各个文件夹内的所有影像文件，输出影像名称为文件夹名）：

import arcpy, os

in_folder = r"F:\Export"
out_folder = r"F:\MERGE"
for root, dirs, files in os.walk(in_folder):
    for dir in dirs:
        print dir
        arcpy.CreateRasterDataset_management(out_folder, dir + ".tif", "", "", "", 3)
        arcpy.WorkspaceToRasterDataset_management(root, out_folder + os.sep + dir + ".tif", nodata_value="0")