在之前这篇博客中介绍了利用GDAL读写带有地理信息的TIFF文件方法。 但最近又有个新的需求：有一个多光谱影像，不同波段有不同的波段名，且波段名就是该波段拍摄时间，要想获得拍摄时间必须从文件中读出波段名。 之前的代码保存的影像默认都是没有设置波段名，所以本篇博客主要介绍对多光谱影像波段名的读和写，以满足这个需求。
1.查看波段名
查看波段名非常简单，在装好GDAL后，利用gdalinfo命令即可，从名字就知道这是一个用来显示文件各种信息的命令。最简单的使用方法就是gdalinfo filename。在控制台中输入命令和文件名，即可以显示文件信息了。如下：
当然gdalinfo也有很多参数，主要是一些输出控制，输出什么、不输出什么这种。完整介绍可以参考官方文档，点击查看。
2.波段名读写
波段名其实是波段的Description属性，读写其实也非常简单，利用GetRasterBand()获取到每个波段后，调用成员函数SetDescription()就是写入，GetDescription()就是读取。Python代码如下：
# coding=utf-8
from osgeo import gdal
from gdalconst import *
def readImage(img_path):
band_data = []
band_name = []
# 以只读方式打开遥感影像
dataset = gdal.Open(img_path, GA_ReadOnly)
if dataset is None:
print("Unable to open image file.")
return band_data
else:
print("Open image file success.\n")
# 读取地理变换参数
param_geoTransform = dataset.GetGeoTransform()
print "GeoTransform info:\n", param_geoTransform, "\n"
# 读取投影信息
param_proj = dataset.GetProjection()
print "Projection info:\n", param_proj, "\n"
# 读取波段数及影像大小
bands_num = dataset.RasterCount
print("Image height:" + dataset.RasterYSize.__str__() + " Image width:" + dataset.RasterXSize.__str__())
print(bands_num.__str__() + " bands in total.")
# 依次读取波段数据
for i in range(bands_num):
# 获取影像的第i+1个波段
band_i = dataset.GetRasterBand(i + 1)
# 获取影像第i+1个波段的描述(名称)
name = band_i.GetDescription()
band_name.append(name)
# 读取第i+1个波段数据
data = band_i.ReadAsArray(0, 0, band_i.XSize, band_i.YSize)
band_data.append(data)
print("band " + (i + 1).__str__() + " read success.")
if name != "":
print "Name:", name
return band_data, param_geoTransform, param_proj, band_name
def writeImage(save_path, bands, geotrans=None, proj=None, names=None):
projection = [
# WGS84坐标系(EPSG:4326)
"""GEOGCS["WGS 84", DATUM["WGS_1984", SPHEROID["WGS 84", 6378137, 298.257223563, AUTHORITY["EPSG", "7030"]], AUTHORITY["EPSG", "6326"]], PRIMEM["Greenwich", 0, AUTHORITY["EPSG", "8901"]], UNIT["degree", 0.01745329251994328, AUTHORITY["EPSG", "9122"]], AUTHORITY["EPSG", "4326"]]""",
# Pseudo-Mercator、球形墨卡托或Web墨卡托(EPSG:3857)
"""PROJCS["WGS 84 / Pseudo-Mercator",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["X",EAST],AXIS["Y",NORTH],EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"],AUTHORITY["EPSG","3857"]]"""
]
if bands is None or bands.__len__() == 0:
return False
else:
# 认为各波段大小相等，所以以第一波段信息作为保存
band1 = bands[0]
# 设置影像保存大小、波段数
img_width = band1.shape[1]
img_height = band1.shape[0]
num_bands = bands.__len__()
# 设置保存影像的数据类型
if 'int8' in band1.dtype.name:
datatype = gdal.GDT_Byte
elif 'int16' in band1.dtype.name:
datatype = gdal.GDT_UInt16
else:
datatype = gdal.GDT_Float32
# 创建文件
driver = gdal.GetDriverByName("GTiff")
dataset = driver.Create(save_path, img_width, img_height, num_bands, datatype)
if dataset is not None:
# 写入仿射变换参数
if geotrans is not None:
dataset.SetGeoTransform(geotrans)
# 写入投影参数
if proj is not None:
if proj is 'WGS84' or \
proj is 'wgs84' or \
proj is 'EPSG:4326' or \
proj is 'EPSG-4326' or \
proj is '4326':
dataset.SetProjection(projection[0]) # 写入投影
elif proj is 'EPSG:3857' or \
proj is 'EPSG-3857' or \
proj is '3857':
dataset.SetProjection(projection[1]) # 写入投影
else:
dataset.SetProjection(proj) # 写入投影
# 逐波段写入数据
for i in range(bands.__len__()):
raster_band = dataset.GetRasterBand(i + 1)
# 设置没有数据的像素值为0
raster_band.SetNoDataValue(0)
if names is not None:
# 设置波段的描述(名称)
raster_band.SetDescription(names[i])
# 写入数据
raster_band.WriteArray(bands[i])
print("save image success.")
return True
if __name__ == '__main__':
# 读取影像
data, geo, prj, name = readImage("geoImage.tif")
# 修改名称
for i in range(len(name)):
name[i] = "This is band " + i.__str__().zfill(2)
# 输出影像
writeImage("out.tif", data, geo, prj, name)
# 读取输出影像
data, geo, prj, name = readImage("out.tif")
3.测试
首先准备了一个包含RGB3个波段的TIFF影像，波段是没有名字的，影像是蚌埠大学城和高铁南站的部分区域，如下。
GDALINFO信息如下:
利用代码首先读取数据，然后对每个波段添加名称信息并保存，保存后的文件GDALINFO信息如下：
可以看到每个波段已经有自己的名字了，还有我们刚刚设置的缺省值。利用代码读取波段名输出如下：
通过代码也是可以成功读取到各波段的名字的，这样便圆满完成了一开始提出的需求。
4.数据交互
需要注意的是，采用本文方法写入的波段名，在ENVI中是无法读取的，如下所示。
如果想要ENVI也可以读取波段名，需要按照ENVI的规则。ENVI是将波段名等头信息存储在同名的.hdr文件中的。所以我们只需要在输出TIFF文件的时候，顺带输出一个对应的.hdr文件，用ENVI就可以打开了。.hdr文件是一个纯文本文件，其内容格式如下：
总的来说基本就是我们在保存时输出的相关信息，只不过用ENVI能识别的格式写一下。所以我们可以增加一个输出.hdr文件的函数，并在writeImage()函数里调用它。这样每次输出的时候就会同时输出影像和.hdr文件了。完整的代码如下：
# coding=utf-8
from osgeo import gdal
from gdalconst import *
import time
def readImage(img_path):
band_data = []
band_name = []
# 以只读方式打开遥感影像
dataset = gdal.Open(img_path, GA_ReadOnly)
if dataset is None:
print("Unable to open image file.")
return band_data
else:
print("Open image file success.\n")
# 读取地理变换参数
param_geoTransform = dataset.GetGeoTransform()
print "GeoTransform info:\n", param_geoTransform, "\n"
# 读取投影信息
param_proj = dataset.GetProjection()
print "Projection info:\n", param_proj, "\n"
# 读取波段数及影像大小
bands_num = dataset.RasterCount
print("Image height:" + dataset.RasterYSize.__str__() + " Image width:" + dataset.RasterXSize.__str__())
print(bands_num.__str__() + " bands in total.")
# 依次读取波段数据
for i in range(bands_num):
# 获取影像的第i+1个波段
band_i = dataset.GetRasterBand(i + 1)
# 获取影像第i+1个波段的描述(名称)
name = band_i.GetDescription()
band_name.append(name)
# 读取第i+1个波段数据
data = band_i.ReadAsArray(0, 0, band_i.XSize, band_i.YSize)
band_data.append(data)
print("band " + (i + 1).__str__() + " read success.")
if name != "":
print "Name:", name
return band_data, param_geoTransform, param_proj, band_name
def writeImage(save_path, bands, geotrans=None, proj=None, names=None):
projection = [
# WGS84坐标系(EPSG:4326)
"""GEOGCS["WGS 84", DATUM["WGS_1984", SPHEROID["WGS 84", 6378137, 298.257223563, AUTHORITY["EPSG", "7030"]], AUTHORITY["EPSG", "6326"]], PRIMEM["Greenwich", 0, AUTHORITY["EPSG", "8901"]], UNIT["degree", 0.01745329251994328, AUTHORITY["EPSG", "9122"]], AUTHORITY["EPSG", "4326"]]""",
# Pseudo-Mercator、球形墨卡托或Web墨卡托(EPSG:3857)
"""PROJCS["WGS 84 / Pseudo-Mercator",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["X",EAST],AXIS["Y",NORTH],EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"],AUTHORITY["EPSG","3857"]]"""
]
if bands is None or bands.__len__() == 0:
return False
else:
# 认为各波段大小相等，所以以第一波段信息作为保存
band1 = bands[0]
# 设置影像保存大小、波段数
img_width = band1.shape[1]
img_height = band1.shape[0]
num_bands = bands.__len__()
# 设置保存影像的数据类型
if 'int8' in band1.dtype.name:
datatype = gdal.GDT_Byte
elif 'int16' in band1.dtype.name:
datatype = gdal.GDT_UInt16
else:
datatype = gdal.GDT_Float32
# 创建文件
driver = gdal.GetDriverByName("GTiff")
dataset = driver.Create(save_path, img_width, img_height, num_bands, datatype)
if dataset is not None:
# 写入仿射变换参数
if geotrans is not None:
dataset.SetGeoTransform(geotrans)
# 写入投影参数
if proj is not None:
if proj is 'WGS84' or \
proj is 'wgs84' or \
proj is 'EPSG:4326' or \
proj is 'EPSG-4326' or \
proj is '4326':
dataset.SetProjection(projection[0]) # 写入投影
elif proj is 'EPSG:3857' or \
proj is 'EPSG-3857' or \
proj is '3857':
dataset.SetProjection(projection[1]) # 写入投影
else:
dataset.SetProjection(proj) # 写入投影
# 逐波段写入数据
for i in range(bands.__len__()):
raster_band = dataset.GetRasterBand(i + 1)
# 设置没有数据的像素值为0
raster_band.SetNoDataValue(0)
if names is not None:
# 设置波段的描述(名称)
raster_band.SetDescription(names[i])
# 写入数据
raster_band.WriteArray(bands[i])
print("save image success.")
writeHdr(save_path, bands, geotrans, proj, names)
return True
def writeHdr(save_path, bands, geotrans=None, proj=None, names=None):
width = int(bands[0].shape[1])
height = int(bands[0].shape[0])
str_lines = "ENVI\ndescription={\n GEO-TIFF File Imported into ENVI "
localtime = "[" + time.strftime("%a %b %d %H:%M:%S %Y", time.localtime()) + "]}\n"
samples = "samples = " + width.__str__() + "\n"
lines = "lines = " + height.__str__() + "\n"
bands = "bands = " + len(bands).__str__() + "\n"
header_offset = "header offset = 0\n"
file_type = "file type = TIFF\n"
data_type = "data type = 1\n"
others = """interleave = bip
sensor type = Unknown
byte order = 0
read procedures = {idl_tiff_read_spatial, idl_tiff_read_spectral}
coordinate system string = {GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]}
wavelength units = Unknown
data ignore value = 0
map points = {"""
# 构造并计算影像四角点经纬度
map_points1 = "\n 1.0000, 1.0000, " + geotrans[0].__str__() + ", " + geotrans[3].__str__() + ",\n"
map_points2 = " " + (width + 1.0).__str__() + ", 1.0000, " + geotrans[0].__str__() + ", " + (
height * geotrans[5] + geotrans[3]).__str__() + ",\n"
map_points3 = " 1.0000, " + (height + 1.0).__str__() + ", " + (
geotrans[0] + width * geotrans[1]).__str__() + ", " + geotrans[3].__str__() + ",\n"
map_points4 = " " + (width + 1.0).__str__() + ", " + (height + 1.0).__str__() + ", " + (
geotrans[0] + width * geotrans[1]).__str__() + ", " + geotrans[3].__str__() + "}\n"
# 依次输出波段名
if names is not None:
band_names = "band names = {\n"
for i in range(len(names)):
band_names += names[i] + ", "
band_names = band_names[:-2] + "}\n"
# 拼接字符串
str_lines += localtime
str_lines += samples
str_lines += lines
str_lines += bands
str_lines += header_offset
str_lines += file_type
str_lines += data_type
str_lines += others
str_lines += map_points1
str_lines += map_points2
str_lines += map_points3
str_lines += map_points4
str_lines += band_names
# 保存文件
path = save_path.split(".")[0] + ".hdr"
output = open(path, 'w')
output.write(str_lines)
output.close()
print "save hdr file success."
if __name__ == '__main__':
# 读取影像
data, geo, prj, name = readImage("geoImage.tif")
# 修改名称
for i in range(len(name)):
name[i] = "This is band " + i.__str__().zfill(2)
# 输出影像
writeImage("out.tif", data, geo, prj, name)
对生成的out.tif首先用gdalinfo方式查看信息，如下：
可以看到波段名已经保存好了。同时看到多了一个.hdr文件，如下。
有了这个文件，再用ENVI打开就会有波段名，如下所示：
完整代码和测试的数据都放到了Github，感兴趣可以看看，点击查看。
5.工具推荐
TuiView是一个轻量级、功能强大的GIS栅格属性表处理软件，官网是这里。之所以推荐它是因为相比于ENVI，它可以直接读TIFF影像中的波段名。关于它的安装这里就不说了，官网主页里就有，很详细，从conda可以直接安装，十分方便。
6.参考资料
- [1]https://gdal.org/programs/gdalinfo.html
- [2]https://stackoverflow.com/questions/32609570/how-to-set-the-band-description-option-tag-of-a-geotiff-file-using-gdal-gdalw
- [3]https://blog.csdn.net/shi_weihappy/article/details/86559860
- [4]https://gdal.org/drivers/raster/gtiff.html?highlight=tifftag_imagedescription
- [5]http://tuiview.org/
- [6]https://gis.stackexchange.com/questions/247878/writing-band-names-into-header-envi-file-type-using-gdal
本文作者原创,未经许可不得转载,谢谢配合
