GDAL：处理Landsat8影像数据

处理Landsat-8 影像

在当前案例中，处理卫星影像包括波段合成、影像拉伸效果、栅格计算、图像锐化效果等。

（1）波段合成

如下案例中，处理了Landsat8数据并创建了衍生产品。在landsat8目录下有5个landsat8的单波段数据，分别是：B2蓝、B3绿、B4红、B5近红外、B8全色波段。

首先进行波段的RGB合成，这里需要使用-separate指定具体的波段，如果直接使用gdal_merge.py脚本会直接加载所有的波段进行合并，运算较慢且可能产生错误。因此创建大文件时还是首选使用创建虚拟栅格gdalbuildvrt。

# 先创建虚拟栅格
cmd /r gdalbuildvrt -o rgb.vrt -separate RT_LC08_L1TP_137042_20190920_20190926_01_T1_2019-09-20_B4.tif RT_LC08_L1TP_137042_20190920_20190926_01_T1_2019-09-20_B3.tif RT_LC08_L1TP_137042_20190920_20190926_01_T1_2019-09-20_B2.tif
# 再合并
cmd /r gdal_translate rgb.vrt rgb.tif -co PHOTOMETRIC=RGB -co COMPRESS=DEFLATE

合并后的图像和栅格直方图如下所示，可以看到影像整体较灰暗，对比度较低，这是因为QGIS默认使用最大最小值进行对比度拉伸，效果并不好。

（2）直方图拉伸

下面使用-scale参数应用直方图拉伸增加对比度，由于大多数像素具有介于 0 和 0.3 之间的值，我们可以选择0是最小值、0.3是最大值，并应用对比度拉伸以使其从 0 变为 255，生成的图像将是一个 8 位彩色图像，其中输入像素值将线性缩放为目标值。但是这个结果仅仅用于可视化展示，不能用于数值分析，因为改变了正在界面展示的像素值。此外，还可以使用-exponent参数进行非线性拉伸，例如指数拉伸，选择介于 0 和 1 之间的指数值将增强低强度值，产生更亮的图像。

# 线性缩放
cmd /r gdal_translate -scale 0 0.3 0 255 -ot Byte rgb.tif rgb_stretch.tif
# 指数拉伸
cmd /r gdal_translate -scale 0 0.3 0 255 -exponent 0.5 -ot Byte rgb.tif rgb_stretch_exp.tif

（3）栅格计算

GDAL通过gdal_calc.py文件提供了栅格代数，输入栅格转为A-Z任意字母作为形参，使用--calc进行表达，此外，还支持了NumPy语法和参数。在运行之前，首先需要查看nodata信息，这对于除数运算等很重要，查看后发现nodata为-999。此外，Windows 用户需要指定gdal_calc.py的完整路径（%CONDA_PREFIX%\Scripts\gdal_calc.py ）并使用 Python 命令运行，Mac/Linux 操作系统用户可以直接键入脚本名称。

# 查看栅格信息
cmd /r gdalinfo -stats RT_LC08_L1TP_137042_20190920_20190926_01_T1_2019-09-20_B4.tif
# 计算NDVI
cmd /r python %CONDA_PREFIX%\Scripts\gdal_calc.py -A RT_LC08_L1TP_137042_20190920_20190926_01_T1_2019-09-20_B5.tif -B RT_LC08_L1TP_137042_20190920_20190926_01_T1_2019-09-20_B4.tif --outfile ndvi.tif --calc="(A-B)/(A+B)" --NoDataValue=-999

如果运行过程中遭遇如下错误，则需要使用 Conda 来安装 setuptools，代码： conda install -c conda-forge setuptools 。

Traceback (most recent call last):
  File "C:\Users\YIN\miniconda\envs\gdal\Scripts\gdal_calc.py", line 4, in <module>
    __import__('pkg_resources').run_script('GDAL==3.10.0', 'gdal_calc.py')
    ~~~~~~~~~~^^^^^^^^^^^^^^^^^
ModuleNotFoundError: No module named 'pkg_resources'

（4）全色锐化效果

Landsat-8 卫星可在 30m 空间分辨率生成图像 红色、绿色、蓝色波段和 在以 15m 空间分辨率 生成Panchromatic 波段。我们可以使用Panchromatic 波段以提高其他波段的分辨率，从而产生更清晰的图像和更多细节。这个过程称为全色锐化。

GDAL 附带一个gdal_pansharpen.py脚本实现 Brovey 算法来计算全色锐化输出。在下面的示例中，将 15m 分辨率的全色波段 （B8） 与RGB合成融合；并且可以应用与以前相同的对比度拉伸，并比较输出图。可以看到新生成的合成要清晰得多，并且可以更好地解析场景中的细节。

# 全色锐化
cmd /r python %CONDA_PREFIX%\Scripts\gdal_pansharpen.py RT_LC08_L1TP_137042_20190920_20190926_01_T1_2019-09-20_B8.tif rgb.tif pansharpened.tif -r bilinear -co COMPRESS=DEFLATE -co PHOTOMETRIC=RGB
# 对比度拉伸
cmd /r gdal_translate -scale 0 0.3 0 255 -exponent 0.5 -ot Byte -a_nodata 0 pansharpened.tif pansharpened_stretch.tif

本文由mdnice多平台发布