GDAL读写多光谱影像的波段名称

Sep 5,2019   10124 words   37 min


在之前这篇博客中介绍了利用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

本文作者原创,未经许可不得转载,谢谢配合

返回顶部