python批量实现对遥感影像的裁剪、重采样、投影【GDAL库】

利用python的第三方库GDAL来显示对遥感影像的批量裁剪、投影和重采样

GDAL库简介

GDAL 是读写大量的栅格空间数据格式的广泛应用的开源库。GDAL/OGR 使用面向对象的 C++ 语言编写,这令该库在支持百余种格式的同时,还具有很高的执行效率。 GDAL/OGR 同时还提供多种主流编程语言的绑定,除了 C 和 C++语言之外,用户还可以在 Perl、python、VB6、Ruby、Java、C# 等语言中调用 GDAL,这令 GDAL 的应用变得非常广泛。GDAL 项目维护了使用 SWIG 生成的 Python 的 GDAL/OGR 绑定。【这段复制官网介绍】
总之，使用GDAL可以帮助我们更高效快捷的处理我们的遥感影像。

GDAL下载

GDAL并非是一个纯粹的第三方库，所以平常的pip安装就不太行了，但是如果你是用的anaconda的话，可以使用python conda install GDAL进行安装。
在基础环境下，我们需要去GDAL官方网站GDAL官网找到自己电脑对应环境的GDAL来进行下载。

下载好之后，同时按住windows+R并输入cmd来进入命令窗口，利用cd命令跳转到我们下载的GDAL库所在文件夹之后，使用命令pip install GDAL-3.4.3-cp310-cp310-win_amd64.whl后面根据自己下载来进行安装。
这样我们就安装好了GDAL。接下来是使用它来实现我们需要的功能。

GDAL库的使用

导入库

import glob
import netCDF4 as nc
import numpy as np
from osgeo import gdal,osr,ogr,gdalconst
import os

利用GDAL库下的Warp处理影像【批量】

参考链接

https://www.jianshu.com/p/821c741ff169

import glob
import netCDF4 as nc
import numpy as np
from osgeo import gdal,osr,ogr,gdalconst
import os


Input_folder = r'输入路径'
Output_folder = r'输出路径'
input_shape = r"裁剪矢量面"

#获取投影参数
InputImage = r'带投影的栅格数据'
Raster = gdal.Open(InputImage, gdal.GA_ReadOnly)
Projection = Raster.GetProjectionRef()

# 读取所有tif数据
data_list = glob.glob(Input_folder + '/*.tif')     

for i in range(len(data_list)):
    data = data_list[i]
    output_raster= Output_folder + os.sep + data.split('\\')[-1].split('_')[0] + '_' + "slip" +".tif"    #文件命名
    # tif输入路径，打开文件
    input_raster = data
    # 矢量文件路径，打开矢量文件
    input_raster=gdal.Open(input_raster)
    ds = gdal.Warp(output_raster, # 裁剪图像保存完整路径（包括文件名）
                   input_raster, # 待裁剪的影像
                   format='GTiff', # 保存图像的格式
                   cutlineDSName = input_shape, # 矢量文件的完整路径
                   cropToCutline = True,
                   copyMetadata=True,
                   creationOptions=['COMPRESS=LZW',"TILED=True"],
                   dstSRS=Projection,
                   resampleAlg=gdalconst.GRA_Cubic,   #三次卷积采样
                   xRes=1000,   #重采样像元1000m
                   yRes=1000,
                   dstNodata=None)
    print(f'{data} is ok!!!')
print("完成！")

总结

相比较于Arcgis中使用Arcpy来说，这个速度是相当好的，也比较方便，但是在我们实际的处理数据过程中也发现了一些问题，gdal下的warp在对遥感影像进行裁剪和重采样时，是先对数据进行矢量裁剪之后再进行重采样的，这就会造成很大的问题，输出的影像往往不是我们需要的影像。这个我还没有更好的解决办法，只能选择一个比较大的矢量面再进行二次裁剪得到我需要的区域影像。