【GDAL基础教程】多张二维tif数据转三维tif数据

【GDAL基础教程】多张二维tif数据转三维tif数据

今天分享一下多张二维单波段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打开可以发现已经成功合并。

如果对你有帮助的话，请‘点赞’、‘收藏’，‘关注’，你们的支持是我更新的动力。欢迎关注【森气笔记】。