如何使用GDAL VRT将图像堆叠为单个GeoTiff中的光栅带？

Question

我有一堆.tif图像，但我想将它们全部堆叠为1.tif图像。如何使用gdalVrt堆叠所有.Tif文件

from osgeo import gdal
from PIL import Image
import numpy as np
from numpy import asarray
import matplotlib.pyplot as plt
import os

os.listdir('../LC08_L1TP_137042_20210301_20210301_01_RT/')


band1=gdal.Open('LC08_L1TP_137042_20210301_20210301_01_RT_B1.TIF')
band2=gdal.Open('LC08_L1TP_137042_20210301_20210301_01_RT_B2.TIF')
band3=gdal.Open('LC08_L1TP_137042_20210301_20210301_01_RT_B3.TIF')
band4=gdal.Open('LC08_L1TP_137042_20210301_20210301_01_RT_B4.TIF') # near infra red
band5=gdal.Open('LC08_L1TP_137042_20210301_20210301_01_RT_B5.TIF') # near infra red
band6=gdal.Open('LC08_L1TP_137042_20210301_20210301_01_RT_B6.TIF')
band7=gdal.Open('LC08_L1TP_137042_20210301_20210301_01_RT_B7.TIF')
band8=gdal.Open('LC08_L1TP_137042_20210301_20210301_01_RT_B8.TIF')
band9=gdal.Open('LC08_L1TP_137042_20210301_20210301_01_RT_B9.TIF')
band10=gdal.Open('LC08_L1TP_137042_20210301_20210301_01_RT_B10.TIF')
band11=gdal.Open('LC08_L1TP_137042_20210301_20210301_01_RT_B11.TIF')
band12=gdal.Open('LC08_L1TP_137042_20210301_20210301_01_RT_BQA.TIF')

band1_array = band4.ReadAsArray()
band2_array = band4.ReadAsArray()
band3_array = band4.ReadAsArray()
band4_array = band4.ReadAsArray()
band5_array = band4.ReadAsArray()
band6_array = band4.ReadAsArray()
band7_array = band4.ReadAsArray()
band8_array = band5.ReadAsArray()
band9_array = band4.ReadAsArray()
band10_array = band4.ReadAsArray()
band11_array = band4.ReadAsArray()
band12_array = band4.ReadAsArray()

stack= (band1_array+band2_array+band3_array+band4_array+band5_array+band6_array+band7_array+band8_array+band9_array+band10_array+band11_array+band12_array)
print("The stacked images is \n \n", plt.imshow(stack))

我确信这不是正确的方法

Answer 1

步骤1.）使用选项“separate=True”创建虚拟光栅（VRT），以将图像堆叠为单独的条带：

from osgeo import gdal
ImageList = ['Band1.tif', 'Band2.tif', 'Band3.tif']  # or use sorted(glob.glob('*.tif')) if input images are sortable
VRT = 'OutputImage.vrt'
gdal.BuildVRT(VRT, ImageList, separate=True, callback=gdal.TermProgress_nocb)

步骤2.）将虚拟光栅（VRT）转换为GeoTiff：

InputImage = gdal.Open(VRT, 0)  # open the VRT in read-only mode
gdal.Translate('OutputImageName.tif', InputImage, format='GTiff',
               creationOptions=['COMPRESS:DEFLATE', 'TILED:YES'],
               callback=gdal.TermProgress_nocb)
del InputImage  # close the VRT