提示：文章末尾有强化学习代码资源 : )
前言
在本教程中，我们将学习使用 Python 和地理空间数据抽象库 GDAL 自动处理栅格数据的基本技术。
栅格文件通常用于存储地形模型和遥感数据及其衍生产品，例如植被指数和其他环境数据集。 栅格文件往往很大（例如，想象一个以 30m x 30m 分辨率覆盖全球的栅格数据集）并且通常以较小的部分（切片）交付和处理。 因此，对此类大型数据集的高效处理和分析需要自动化。 在教程中，您将学习如何读写常见的栅格格式，以及如何使用 Python 中的 GDAL/OGR API 和 GDAL 命令行实用程序对一批文件进行基本的栅格数据处理。
在本文总结部分将一些实用的python剪裁批处理，python影像镶嵌以及一些Arcpy相关代码进行链接下载，供有需求的同学进行学习与使用。
一、使用 GDAL 在 Python 中打开栅格文件
使用 GDAL，可以在 Python 中读取和写入多种不同的栅格格式。 导入 GDAL 模块时，Python 会自动注册所有已知的 GDAL 驱动程序以读取支持的格式。 最常见的文件格式包括 TIFF 和 GeoTIFF、ASCII Grid 和 Erdas Imagine .img 文件。
Landsat 8 波段作为单独的 GeoTIFF 文件存储在原始包中。 每个波段包含来自不同电磁波谱范围的表面反射率信息。
from osgeo import gdal
filepath = r"LandsatData/LC81910182016153LGN00_sr_band4.tif"
# Open the file:
raster = gdal.Open(filepath)
# Check type of the variable 'raster'
type(raster)(1) 读取栅格文件属性
卫星图像现在作为 GDAL 数据集对象存储在变量栅格中。 让我们仔细看看文件的属性：
# Projection
raster.GetProjection()
# Dimensions
raster.RasterXSize
raster.RasterYSize
# Number of bands
raster.RasterCount
# Metadata for the raster dataset
raster.GetMetadata()(2) 获得栅格影像波段
在我们的例子中，Landsat 8 的所有波段都存储为单独的文件。 rasterCount 为 1，因为我们只打开了一个包含 Landsat 8 波段 4 的 GeoTiff。但是，卫星图像的不同波段通常堆叠在一个栅格数据集中，在这种情况下 rasterCount 会大于 1。
# Read the raster band as separate variable
band = raster.GetRasterBand(1)
# Check type of the variable 'band'
type(band)
# Data type of the values
gdal.GetDataTypeName(band.DataType)现在我们有一个存储在变量 band 中的 GDAL Raster Band 对象。
波段的数据类型可以在像素数据类型的 GDAL 文档的帮助下进行解释。 无符号整数始终等于或大于零，有符号整数也可以存储负值。 例如，一个无符号的 16 位整数可以存储 2^16 (=65,536) 个值，范围从 0 到 65,535。
（3） 波段统计
接下来，让我们看一下存储在 band 中的值。 在能够打印出任何信息之前，您可能需要计算栅格的统计数据。
# Compute statistics if needed
if band.GetMinimum() is None or band.GetMaximum()is None:
band.ComputeStatistics(0)
print("Statistics computed.")
# Fetch metadata for the band
band.GetMetadata()
# Print only selected metadata:
print ("[ NO DATA VALUE ] = ", band.GetNoDataValue()) # none
print ("[ MIN ] = ", band.GetMinimum())
print ("[ MAX ] = ", band.GetMaximum())
二、将栅格文件作为数值数组打开
在访问地理空间栅格数据方面，GDAL 是一个功能强大的库，但它没有提供很多计算功能。 对于更高级的计算，我们将以数值数组的形式读入栅格数据，以便使用 NumPy 库中的功能。
如果您想将现有的 Gdal 数据集或 Band 转换为 numpy 数组，您可以使用 ReadAsArray 方法将其转换：
In [ ]:
# read raster data as numeric array from Gdal Dataset
rasterArray = raster.ReadAsArray()那有什么区别呢？
In [ ]:
#Check the datatype of variables
type(rasterArray)
type(raster)如您所见，我们现在已将相同的栅格数据存储到两种不同类型的变量中：
-raster 是一个 Gdal 数据集
-rasterArray 是一个 numpy 数组
GDAL 库还附带一个模块 gdal_array，用作 NumPy 和 GDAL 之间的接口。 Gdal_array 直接从文件读取和写入栅格文件到 NumPy 数组：
from osgeo import gdal_array
#from osgeo import gdalnumeric
# Read raster data as numeric array from file
rasterArray = gdal_array.LoadFile(filepath)（1）注意NODATA数据
# What is the minimum value of the array?
rasterArray.min() #-9999如您所见，numpy 数组仍然包含原始无数据值。 如果不注意这些，任何计算都将是错误的。
import numpy as np
# Get nodata value from the GDAL band object
nodata = band.GetNoDataValue()
#Create a masked array for making calculations without nodata values
rasterArray = np.ma.masked_equal(rasterArray, nodata)
type(rasterArray)
# Check again array statistics
rasterArray.min()现在您有一个二维数组，可以进行进一步的计算。 但是我们将使用一个非常简单的命令行工具 gdal_calc.py。
(2)关闭栅格数据
如果您想释放资源并从内存中删除不必要的变量，则在代码中间关闭现有的 GDAL 对象可能很有用。
raster = None
band = None三、GDAL 命令行实用程序
我们现在已经测试了 Python GDAL/OGR API 中用于读取和检查栅格文件的一些基本函数。 然而，GDAL 还包括其他强大的数据转换和处理功能，这些功能没有直接在库中实现。 我们将仔细研究几个这样的函数：
- gdalwarp 用于裁剪、镶嵌、重投影和其他过程
- gdal_merge.py 用于拼接/堆叠图像
- gdal_calc.py 用于栅格计算
这些工具需要从终端/命令提示符或作为 Python 中的子进程运行。 我们现在将在终端窗口中快速测试这些工具。
(1) 使用 gdalwarp 裁剪图像
在其他技巧中，gdalwarp 是一个非常方便的快速裁剪图像的工具。 我们现在将练习如何基于边界框裁剪卫星图像波段。使用选项 -te 指定输出文件的所需范围：
gdalwarp -te xmin ymin xmax ymax inputfile.tif outputfile.tif接下来，让我们同时对所有其余波段重复剪裁。 为此，我们将使用 Python 为场景中的每个光谱带生成命令。
import glob
import os
# List filepaths for all bands in the scence
FileList = glob.glob(os.path.join(r'/home/geo/LandsatData','*band*.tif'))
# Define clipping extent
xmin, ymin, xmax,ymax = (0, 0, 0, 0) # INSERT HERE THE CORRECT COORDINATES
# Generate gdalwarp command for each band
command = ""
for fp in FileList:
inputfile = fp
outputfile = inputfile[:-4] + "_clip.tif"
command += "gdalwarp -te %s %s %s %s %s %s \n" % (xmin, ymin, xmax, ymax, inputfile, outputfile)
# Write the commands to an .sh file
cmd_file = "ClipTurkufromLandsat.sh"
f = open(os.path.join(cmd_file), 'w')
f.write(command)
f.close()注意：如果您在 Windows 环境中工作，请将 .sh 扩展名更改为 .bat。
运行上述脚本后，您的工作目录中应该有一个文件“ClipTurkufromLandsat.sh”。 打开文件（使用文本编辑器）并检查命令是否已正确写入文件。
接下来，在终端窗口中运行该文件：
bash ClipTurkufromLandsat.sh(2) 使用 gdal_merge.py 生成layer stack
裁剪图像后，可以堆叠波段 3（绿色）、4（红色）和 5（近红外）以假彩色合成。 使用 gdal_merge.py 合并图层并使用 -separate 选项指示您希望将输入保存为输出文件中的单独波段。
让我们尝试在 python 中将命令作为子进程运行：
import os
# Define input and output files
inputfiles = "band3_clip.tif band4_clip.tif band5_clip.tif"
outputfile = "Landsat8_GreenRedNir.tif"
# Generate the command
command = "gdal_merge.py -separate %s -o %s" % (inputfiles, outputfile)
# Run the command. os.system() returns value zero if the command was executed succesfully
os.system(command)(3) Gdal_calc.py
Gdal_calc.py 是一个命令行栅格计算器，可用于栅格数据的简单重复计算。 打开终端窗口并输入：
Gdal_calc.py然后按回车。 您应该看到有关该工具的用法和选项的说明。gdal_calc.py 的基本语法如下：
gdal_calc.py -A input1.tif - B input2.tif [other_options] --outfile=outputfile.tif
从其他选项中，至少要注意用于指定计算语法的参数 --calc 和用于控制输出文件大小的 --creation-option（或 --co）参数：
- 在两个输入文件的情况下 --calc="A+B" 会将文件 A 和 B 添加到一起。
- 默认情况下，输出文件往往很大，这将很快导致磁盘大小和内存出现问题。 使用 gdal_calc.py 您可以添加参数 --co="COMPRESS=LZW" 以减小输出文件的大小。
总结
本博客介绍了Python GDAL对遥感数据的处理。下面附三个代码链接供有需要的同学进行强化学习。
（1）基于Geopandas的矢量文件裁剪矢量文件方法：矢量文件剪裁矢量文件（Python）-Python文档类资源-CSDN下载
（2）
基于Python-GDAL的遥感影像镶嵌脚本：
基于Python-GDAL的遥感影像镶嵌脚本_pythongdal库镶嵌-Python文档类资源-CSDN下载
（3）基于Arcpy以及GDAL的批量分类栅格图转shp矢量文件程序：
遥感影像分类结果栅格转矢量（Python）_python栅格转矢量,分类结果转矢量-Python文档类资源-CSDN下载
