在网上找资料的时候，发现了一位大神的个人网站(http://zhaoxuhui.top/）,上面分享了许多有关地信知识的教程，非常inspiring，大家有兴趣的话可以去关注一下。但是，对于”基于Gdal包在python中合成多波段图像并给波段命名”的相关问题，我个人认为大神讲的有些绕。所以，此篇是对大神文章的进一步解释和总结（分别设定三个难度依次递增的情景）。我本身也只是一个初学者，疏漏之处在所难免，希望各位内行不要笑话，多多提建议，帮助小张提高自己！！！

本篇文章将不对gdal的基础语法进行过多的解释。不熟悉gdal的朋友们可以先看看犹他州里大学所编撰的这个手册(https://www.osgeo.cn/python_gdal_utah_tutorial/index.html#)。大概扫几眼后，就可以对gdal有个初步的认识。废话不再多说，进入正题。

情景:

一．现有多个单波段图像，需要将其合成为一个多波段图像并给每个波段命名

二．现有一个没有给每个波段命名的多波段图像，需要在原图形的基础上，给原图形的各波段命名

三．现有一个没有给每个波段命名的多波段图像，不能改变原图形，复制一个新图形后，再给新图形的各波段命名

———————————————————————————————————————————
一．现有多个单波段图像，需要将其合成为一个多波段图像并给每个波段命名
思路其实很简单。首先，生成一个“尺寸空间”等信息和单波段图像相同的“空的多波段图像”，然后逐个将单波段图像存储进“空的多波段图像”，同时给波段命名。例如，我们有有23张由modis13q1提取出的单波段ndvi图像（图一），现在想把他们合成为一张有名字的多波段的ndvi图像。让我们在代码中理解这个思路。

23张2001年的单波段NDVI图像

打开python，载入此操作需要的包

#####载入此操作需要的包#####
import glob
import os
import numpy as np
import pandas as pd
import os
import gdal
from osgeo import gdal
from gdalconst import *

定义一个返回为“栅格文件”，以及“栅格文件的一个band”的function

#返回的“栅格文件”以及”栅格文件的一个band”的信息将被用来创建“空的多波段图像”
#此function是直接用大神的，但其实可能很多人都用过类似的 
def get_dataset_band(bandfile):   
    input_dataset = gdal.Open(bandfile)
    input_band = input_dataset.GetRasterBand(1)
    return input_dataset, input_band

将单波段NDVI数据，以及数据的名字分别储存在名为bandfile，bandnames的list里

inws= "F:\\Modis\\NDVI" #请自行替换为自己单波段数据储存的位置    
NDVIS=glob.glob(os.path.join(inws, "*.tif")) 
bandfile = []  #bandfile为储存单波段NDVI数据的list
bandnames=[]  #bandnames为储存单波段NDVI数据名称的list，当然原名称太长，我们之后要改造
for ndvi in NDVIS:
       bandfile.append(ndvi)
       names=os.path.basename(ndvi).split('.')[0]+"_"+os.path.basename(ndvi).split('.')[1][1:8]+"_"+os.path.basename(ndvi).split('.')[2]+".tif" #对名字进行改造，这样的改造结果如'MOD13Q1_2011353_h25v04.tif'
       bandnames.append(names)

基于单波段的NDVI图像的信息创建空的多波段图像栅格文件

inputdata=[]
for i in range(len(bandfile)):
              inputdata.append(get_dataset_band(bandfile[i]))
inputdataset_1, inputband_1 = inputdata[0]  #inputdataset_1为栅格文件，inputband_1为栅格文件下波段的信息，之后我们要基于inputdataset_1和inputband_1这两个中间变量的信息（如尺寸，坐标等）来定义**空的多波段图像**栅格文件
#创建驱动
file_driver=gdal.GetDriverByName('Gtiff')
#定义**空的多波段图像**的存储路径
output_path="F:\\STSG\\reconstructedNDVI"
#定义**空的多波段图像**的名称
output_name="NDVI"+"_"+location+"_"+year
#创建**空的多波段图像**栅格文件
output_dataset = file_driver.Create(output_path+"\\"+output_name+".tif", inputband_1.XSize, inputband_1.YSize, 23, inputband_1.DataType)
#基于前述的栅格文件来定义**空的多波段图像**空间信息等等
output_dataset.SetProjection(inputdataset_1.GetProjection())
output_dataset.SetGeoTransform(inputdataset_1.GetGeoTransform())

将23个单波段的NDVI逐个写入空的多波段图像栅格文件，并顺手建立一个同名的hdr头文件

for i in range(len(bandfile)):
     inputband_data = inputdata[i][1].ReadAsArray()
     raster_band =output_dataset.GetRasterBand(i + 1)
     raster_band.SetDescription(bandnames[i])  #写入各个波段的名字
     # 写入数据
     raster_band.WriteArray(inputband_data)
     output_dataset.BuildOverviews('average',[2,4,8,16,32])  #建立pyramid
     #截至到这里，数据以及名字都已经被写进**空的多波段图像**栅格文件里了，但是这样写入的名字无法在ENVI软件中被读取，所以还需要再顺手给这个**多波段图像**栅格文件写个hdr头文件
     # 写hdr头文件，bands_num（波段数）和bandnemes（波段名称）不能少；其他无所谓
     h1='ENVI'
     h2='description ={'
     h3='File Imported into ENVI. }'
     h4='img_width = '+str(inputband_1.XSize)    
     #inputband_1.XSize来求每一波段的列数，前面定义过
     h5='img_height = '+str(inputband_1.YSize)   
     #inputband_1.YSize来求每一波段的行数，前面定义过
     h6='bands = '+str(23)                               #23为波段数
     h7='header offset = 0'
     h8='file type =TIFF'
     h9='sensor type = Unknown'
     h10='band names = {'
     h11=','.join(bandnames)+'}'   #用这个join语句，删掉list两侧的[]，以及波段名称两端的双引号""
     h=[h1,h2,h3,h4,h5,h6,h7,h8,h9,h10,h11]
     doc = open(output_path+"\\"+output_name+".hdr", 'w')   
     for i in range (0,11):
              print(h[i], end='', file=doc)
              print('\n', end='', file=doc)
     doc.close()       
     #输出为同名的hdr头文件，并且此文件要与多波段栅格文件在一个文件夹里面。

让我们在ENVI中检验一下程序运行的结果，看是不是波段名称已经被写进去了

在ENVI中检验一下程序运行的结果，看是不是波段名称已经被写进去了

顺便看一看生成的hdr头文件的内容

看一看生成的hdr头文件的内容

总结
其实只是非常基础的一些地信知识，以及非常基础的python语句。希望可以帮助到像小张一样搞地貌水文环境等的学生们。写到这里感觉有点累了。今天就到此为止吧。（主要是急着去追综艺）情景二和三以后再更新。感谢阅读哈哈哈。

基于Gdal包在python中合成多波段图像，并给波段命名 [1]

基于Gdal包在python中合成多波段图像，并给波段命名 [1]

情景: