今天分享一下多张二维单波段tif数据合并为一张三维多波段tif数据的脚本，话不多说，详见代码。
原数据
# -*- encoding: utf-8 -*-
'''
@File : tif_3d.py
@Time : 2022/04/26 21:10:58
@Author : HMX
@Version : 1.0
@Contact : kzdhb8023@163.com
'''
# here put the import lib
from osgeo import gdal
import os
import glob
import numpy as np
import time
t1 = time.time()
tiflist = []
filepath = r'E:\Project\XINAN\WWLLN_NEW\CG\YEAR'
for file in glob.glob(os.path.join(filepath,'*.tif')):
print(file)
tiflist.append(file)
inds = gdal.Open(tiflist[0])
inband = inds.GetRasterBand(1)
indriver = gdal.GetDriverByName('GTiff')
outds = indriver.Create(r'D:\公众号\NO.10\out.tif',inband.XSize,inband.YSize,len(tiflist),inband.DataType)
outds.SetProjection(inds.GetProjection())
outds.SetGeoTransform(inds.GetGeoTransform())
ds = inds.ReadAsArray()
m,n = ds.shape[0],ds.shape[1]
ds_new = np.zeros((len(tiflist),m,n))
for r in range(len(tiflist)):
rds = gdal.Open(tiflist[r]).ReadAsArray()
nodatavalue = rds[0,0]
rds[rds==nodatavalue] = np.nan
ds_new[r,:,:] = rds[:,:]
outband = outds.GetRasterBand(r+1)
outband.WriteArray(ds_new[r,:,:])
outband.SetDescription(os.path.basename(tiflist[r]))
inds = None
rds = None
t2 = time.time()
print('共计用时：{:.2f}s'.format(t2-t1))
结果
利用gis打开可以发现已经成功合并。
如果对你有帮助的话，请‘点赞’、‘收藏’，‘关注’，你们的支持是我更新的动力。欢迎关注我的公众号【森气笔记】。
