使用Arcgis对landsat7和landsat8进行配准、裁剪、导出

一、安装ArcGIS10.8并在pycharm中配置编译器

1、安装ArcGIS，网上有很多教程，大家自行搜索
2、在pycharm中配置编译环境
（1）首先确保在安装ArcGIS时已经安装了python2.7，找到该文件夹，我的安装路径是默认的：C:\Python27\ArcGIS10.8
（2）打开pycharm，新建项目，选择“文件”->“设置”->“项目：XXX”->“项目解释器”，点击右面设置，选择add

（3）选择“现有环境”，找到刚刚安装的python27路径，选择ArcGIS10.8里面的python.exe，点击确定，最后点击应用即可使用。

二、使用arcpy库

修改完编译器之后在pycharm中“import arcpy”就没有问题啦~
下面的代码是使用arcpy将多个单波段合称为多波段：

# -*- coding: utf-8 -*-
import os, arcpy, datetime
# 读取和存放文件夹的位置
layerspath = r'D:\cloud_dataset\Peizhun\boduanhecheng_input\landsat8_2'  # 波段存放文件夹
path = r'D:\cloud_dataset\Peizhun\boduanhecheng_output'  # 输出结果存放文件夹

# start calculate time
startTime = datetime.datetime.now()
print startTime

# set the intermediate data folder
intermediateDataPath = path + "\\" + "IntermediateData"
# set result data folder
resultDataPath = path + "\\" + "Result"

# determine if the folder exists
if os.path.exists(intermediateDataPath):
    print "IntermediateData floder exists"
else:
    # create a intermediate data floder
    arcpy.CreateFolder_management(path, "IntermediateData")

if os.path.exists(resultDataPath):
    print("Result floder exists")
else:
    # create a result floder
    arcpy.CreateFolder_management(path, "Result")
print("-----------------------------------------------------------")
print("Under calculation......")
print("Please do not close the window.")

# 这段为数据处理代码
# function:按顺序拼接字符串,作为波段组合函数的的参数
def getInputParam(files):
    ret = ""
    rets = []
    layers_num = len(files)

    i = 1
    for file in files:  # 拼接其他的
        print file
        ret = ret + file + ";"
        if i % 3 == 0:  # 这里的9代表读取9个数据进行波段相加，如果你要合成3个波段，修改即可
            ret = ret[:-1]  # 去除最后的;
            rets.append(ret)
            ret = ""
        i += 1
    return rets

arcpy.env.workspace = layerspath  # 栅格影像所在文件夹
print layerspath + "路径下包含以下数据:"

arcpy.env.overwriteOutput = True
files = os.listdir(layerspath)
print "-----------------------------------------"
input_param = getInputParam(files)
print "-----------------------------------------"
print input_param
for param in input_param:
    print "------------param------------------------"
    print param
    print "开始波段组合运算！"
    print "-----------------------------------------"
    # landsat8
    resultfileName = param[:37]
    # landsat7
    # resultfileName = param[:21]
    arcpy.CompositeBands_management(param, resultDataPath + "\\" + resultfileName + "_StackLayers.TIF")
    print resultfileName + "计算完成！"

print "Finish!"
endTime = datetime.datetime.now()
print "Time use: " + str((endTime - startTime).seconds) + " (second)"

三、使用ArcGIS进行配准、裁剪并导出

1、配准
（1）首先，将刚刚合成的tif作为图层导入到ArcGIS中


此时，两张遥感影像已经在坐标上自动对齐了
2、创建裁剪矢量图–shp
（1）首先，点击目录，然后在保存矢量图的文件夹上右键选择新建shapefile

对其进行命名，要素类型选择折线：

此时，会弹出警告窗口，返回上一步选择合适的坐标系，如何选择可看文末

（2）编辑矢量图
选择“编辑器”，然后点击“开始编辑”

选中“裁剪2”，点击继续

右边菜单栏点击刚刚新建的裁剪2，然后选择下面的“构造工具”，这里选择矩形

绘制的裁剪区域如下图所示：

3、裁剪
点击“地理处理”，选择“ArcToolbox”，选择“数据管理工具”->“栅格”->“栅格处理”->“裁剪”

输入栅格表示待裁剪影像，输出范围即刚刚创建的shp矢量矩形，勾选“使用输入要素裁剪几何”，设置输出文件名称即可


4、加载刚刚裁剪后的栅格数据
双击RGB，调整好合适的rgb波段

在刚刚的图层上点击右键，选择“数据”->“导出数据”

选择“使用渲染器”，“强制为RGB颜色”，然后设置路径、名称以及格式，点击保存即可

如果想要保存为png或者jpg格式的话，在“格式”处进行修改即可

至此，大功告成！

ps：小问题
进行裁剪时，显示“ArcGIS出现：ERROR 001566：创建栅格数据集失败。原因是：裁剪要素超出了栅格范围。”，原因是两个图层坐标系不是一样的，解决办法：
（1）首先右键待裁剪图层，选择属性

查看XY坐标系，复制

在创建shp文件时，选择相同坐标系即可