实习期间需要批处理一批 Sentinel-2 哨兵影像,内容包括波段变换与影像裁剪等。自己琢磨了一段时间,完美达成了需求,记录一下心得。
读取影像并获取波段
打开下载好的哨兵影像的文件夹,直接使用gdal.Open()
即可读取其中的.xml
文件,打开包含不同传感器的数据集。
然后可以使用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三个波段。
所谓的“波段组合”就是要在用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,而自己要跑几百张图。经实践,这样运行一分钟能处理两张图,跑几百张图也挺快的。