首页 学习区

实习期间需要批处理一批 Sentinel-2 哨兵影像,内容包括波段变换与影像裁剪等。自己琢磨了一段时间,完美达成了需求,记录一下心得。

读取影像并获取波段

打开下载好的哨兵影像的文件夹,直接使用gdal.Open()即可读取其中的.xml文件,打开包含不同传感器的数据集。

image-20210726231042705

然后可以使用GetSubDatasets()函数来获取不同分辨率传感器下的全部波段。至于具体的波段对应,还需要进行进一步的调试,详细见下方代码。

from osgeo import gdal
from gdalconst import GRIORA_Bilinear
from glob import glob
import os

ds = gdal.Open(src) # src即为.xml文件
ds_list = ds.GetSubDatasets()
    
ds_10m = gdal.Open(ds_list[0][0])  # 其中含有B4
ds_20m = gdal.Open(ds_list[1][0])  # 其中含有B8A、B12

重采样与波段组合

对于哨兵影像而言,不同的波段来自不同的传感器,自然也有着有不同的分辨率。要将属于不同分辨率的波段合成在一张图中,首先要做的就是重采样!我比较喜欢使用ReadAsArray()函数直接在读取栅格时进行重采样,一步到位非常方便。

band_B4 = ds_10m.GetRasterBand(1)
band_B8A = ds_20m.GetRasterBand(4)
band_B12 = ds_20m.GetRasterBand(6)
bd_list = [band_B12, band_B8A, band_B4]

xsize_10m = ds_10m.RasterXSize
ysize_10m = ds_10m.RasterYSize

rst_B4 = band_B4.ReadAsArray(0, 0, xsize_10m, ysize_10m)
# 使用`ReadAsArray`直接进行重采样
rst_B8a = band_B8A.ReadAsArray(buf_xsize = xsize_10m, buf_ysize = ysize_10m, 
                               resample_alg = GRIORA_Bilinear)
rst_B12 = band_B12.ReadAsArray(buf_xsize = xsize_10m, buf_ysize = ysize_10m, 
                               resample_alg = GRIORA_Bilinear)
rst_list = [rst_B12, rst_B8a, rst_B4]

关于波段组合,我参考了这篇文章,这里以短波红外为例,选取B12、B8A、B4三个波段。

image-20210801200002114

所谓的“波段组合”就是要在用python处理的过程中按照R、G、B的顺序依次写入波段。另外也要注意结合实际情况创建栅格文件的套路,具体代码如下:

tempfile = "somepath"
temp_tif = driver.Create(tempfile, xsize_10m, ysize_10m, 3, gdal.GDT_Float32)
temp_tif.SetProjection(ds_10m.GetProjection())  # 不同波段的投影相同
temp_tif.SetGeoTransform(ds_10m.GetGeoTransform())  # 统一设置为10m分辨率
for i in range(3):
    temp_tif.GetRasterBand(i+1).WriteArray(rst_list[i])
    temp_tif.GetRasterBand(i+1).SetDescription(bd_list[i].GetDescription())
    
if(temp_tif.GetProjection()[-5:-3]=='48'):
    clip_shp = "D:/Desktop/dian_48N.shp"
elif(temp_tif.GetProjection()[-5:-3]=='47'):
    clip_shp = "D:/Desktop/dian_47N.shp"
else:
    print("出现了额外投影,为: "+temp_tif.GetProjection()[-8:-3])
        
temp_tif.FlushCache()
temp_tif = None

影像裁剪

使用gdal.Warp进行影像裁剪,这里需要一个面属性的.shp文件作为裁剪范围,格外注意投影要与哨兵影像统一!这里我是手动选了两块不同投影tempfile,在Arcmap里手画的范围,这样就要求多一步判断投影,相关代码写在了上一段内容里。此外,还需要注意设置cropToCutline=True以去除影像裁剪范围以外的黑边。

gdal.Warp(out_tif, # 保存的完整路径
          tempfile,  # 待裁剪影像
          format='GTiff', 
          cutlineDSName = in_shp, # 用于裁剪的shp文件
          cropToCutline=True) # 去黑边

if os.path.exists(tempfile):
    os.remove(tempfile)
else:
    print("no tempfile exists!")

我批量处理影像的方法简单而粗暴:用glob.glob确定原始哨兵影像的目录并存储到列表里,依次提取出三波段合成的tempfile,从tempfile中裁剪出成果tif图,再删除tempfile。这样做的原因主要是节省空间,毕竟一张tempfile的大小就有1.6个G,而自己要跑几百张图。经实践,这样运行一分钟能处理两张图,跑几百张图也挺快的。




文章评论

目录