之前在这篇博客中简单介绍了GDAL的使用，但并不是很完整。本文给出一个相对规范、比较完整，具有一定通用性的GDAL读写tif影像的代码，以便以后使用。代码以图像裁剪和波段融合/分离两个小功能为例进行介绍。
1.影像裁剪代码
# coding=utf-8
from osgeo import gdal
from gdalconst import *
import numpy as np
def writeTif(bands, path):
if bands is None or bands.__len__() == 0:
return
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(path, img_width, img_height, num_bands, datatype)
if dataset is not None:
for i in range(bands.__len__()):
dataset.GetRasterBand(i + 1).WriteArray(bands[i])
print("save image success.")
def readTifImage(img_path):
data = []
# 以只读方式打开遥感影像
dataset = gdal.Open(img_path, GA_ReadOnly)
if dataset is None:
print("Unable to open image file.")
return data
else:
print("Open image file success.")
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个波段数据
band_data = band_i.ReadAsArray(0, 0, band_i.XSize, band_i.YSize)
data.append(band_data)
print("band " + (i + 1).__str__() + " read success.")
return data
def readTifImageWithWindow(img_path, start_x, start_y, x_range, y_range):
data = []
# 以只读方式打开遥感影像
dataset = gdal.Open(img_path, GA_ReadOnly)
if dataset is None:
print("Unable to open image file.")
return data
else:
print("Open image file success.")
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个波段数据
band_data = band_i.ReadAsArray(start_x, start_y, x_range, y_range)
data.append(band_data)
print("band " + (i + 1).__str__() + " read success.")
return data
bands_data = readTifImageWithWindow("D1923-561-32_rc.tif", 0, 8000, 5000, 4000)
writeTif(bands_data, "cut_img.tif")
代码比较简单，读写功能包装成了函数方便调用。整个脚本的功能是对tif影像进行裁剪提取ROI。裁剪后的影像如下。
这样便实现了单波段遥感影像的读取与存储。
2.多波段融合代码
# coding=utf-8
from osgeo import gdal
from gdalconst import *
import numpy as np
def writeTif(bands, path):
if bands is None or bands.__len__() == 0:
return
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(path, img_width, img_height, num_bands, datatype)
if dataset is not None:
for i in range(bands.__len__()):
dataset.GetRasterBand(i + 1).WriteArray(bands[i])
print("save image success.")
def readTifImage(img_path):
data = []
# 以只读方式打开遥感影像
dataset = gdal.Open(img_path, GA_ReadOnly)
if dataset is None:
print("Unable to open image file.")
return data
else:
print("Open image file success.")
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个波段数据
band_data = band_i.ReadAsArray(0, 0, band_i.XSize, band_i.YSize)
data.append(band_data)
print("band " + (i + 1).__str__() + " read success.")
return data
def readTifImageWithWindow(img_path, start_x, start_y, x_range, y_range):
data = []
# 以只读方式打开遥感影像
dataset = gdal.Open(img_path, GA_ReadOnly)
if dataset is None:
print("Unable to open image file.")
return data
else:
print("Open image file success.")
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个波段数据
band_data = band_i.ReadAsArray(start_x, start_y, x_range, y_range)
data.append(band_data)
print("band " + (i + 1).__str__() + " read success.")
return data
# 蓝波段
band4_data = readTifImageWithWindow("D1923-561-04_rc_out.tif", 1300, 13000, 3000, 4000)
# 绿波段
band11_data = readTifImageWithWindow("D1923-561-11_rc_out.tif", 1300, 13000, 3000, 4000)
# 红波段
band22_data = readTifImageWithWindow("D1923-561-22_rc_out.tif", 1300, 13000, 3000, 4000)
# 近红外波段
band32_data = readTifImageWithWindow("D1923-561-32_rc.tif", 1300, 13000, 3000, 4000)
# 构造一个list用于存放各波段，存放顺序为B、G、R、NIR
join_data = []
join_data.extend(band4_data)
join_data.extend(band11_data)
join_data.extend(band22_data)
join_data.extend(band32_data)
writeTif(join_data, "join.tif")
利用波段配准后的数据，波段融合成一个包含B、G、R、NIR，4个波段的tif文件。合成后用ENVI打开如下。波段列表中一共有4个波段，说明合并成功了。
依次点击，装载RGB波段效果如下。
NIR、G、B波段效果如下。
在近红外波段下植被呈现红色。
3.波段拆分
# coding=utf-8
from osgeo import gdal
from gdalconst import *
import numpy as np
def writeTif(bands, path):
if bands is None or bands.__len__() == 0:
return
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(path, img_width, img_height, num_bands, datatype)
if dataset is not None:
for i in range(bands.__len__()):
dataset.GetRasterBand(i + 1).WriteArray(bands[i])
print("save image success.")
def readTifImage(img_path):
data = []
# 以只读方式打开遥感影像
dataset = gdal.Open(img_path, GA_ReadOnly)
if dataset is None:
print("Unable to open image file.")
return data
else:
print("Open image file success.")
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个波段数据
band_data = band_i.ReadAsArray(0, 0, band_i.XSize, band_i.YSize)
data.append(band_data)
print("band " + (i + 1).__str__() + " read success.")
return data
def readTifImageWithWindow(img_path, start_x, start_y, x_range, y_range):
data = []
# 以只读方式打开遥感影像
dataset = gdal.Open(img_path, GA_ReadOnly)
if dataset is None:
print("Unable to open image file.")
return data
else:
print("Open image file success.")
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个波段数据
band_data = band_i.ReadAsArray(start_x, start_y, x_range, y_range)
data.append(band_data)
print("band " + (i + 1).__str__() + " read success.")
return data
# 装载数据
bands_data = readTifImage("join.tif")
# 数据拆分
# 注意虽然是一个波段，但是依然要以list的形式传入
band_b = [bands_data[0]]
band_g = [bands_data[1]]
band_r = [bands_data[2]]
band_nir = [bands_data[3]]
writeTif(band_b, "band_b.tif")
writeTif(band_g, "band_g.tif")
writeTif(band_r, "band_r.tif")
writeTif(band_nir, "band_nir.tif")
运行效果如下，下图是G波段单独拆分出来的结果。
4.功能整合脚本
最后将这几个功能写到了一块变成了一个脚本，方便使用。对于日常使用，就不用再打开很大的ENVI了，直接运行脚本即可，可以节约点时间。
# coding=utf-8
from osgeo import gdal
from gdalconst import *
import sys
import os
def writeTif(bands, path):
if bands is None or bands.__len__() == 0:
return
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(path, img_width, img_height, num_bands, datatype)
if dataset is not None:
for i in range(bands.__len__()):
dataset.GetRasterBand(i + 1).WriteArray(bands[i])
print("save image success.")
def readTifImage(img_path):
data = []
# 以只读方式打开遥感影像
dataset = gdal.Open(img_path, GA_ReadOnly)
if dataset is None:
print("Unable to open image file.")
return data
else:
print("Open image file success.")
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个波段数据
band_data = band_i.ReadAsArray(0, 0, band_i.XSize, band_i.YSize)
data.append(band_data)
print("band " + (i + 1).__str__() + " read success.")
return data
def readTifImageWithWindow(img_path, start_x, start_y, x_range, y_range):
data = []
# 以只读方式打开遥感影像
dataset = gdal.Open(img_path, GA_ReadOnly)
if dataset is None:
print("Unable to open image file.")
return data
else:
print("Open image file success.")
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个波段数据
band_data = band_i.ReadAsArray(start_x, start_y, x_range, y_range)
data.append(band_data)
print("band " + (i + 1).__str__() + " read success.")
return data
def findAllFiles(root_dir, filter):
print("Finding files ends with \'" + filter + "\' ...")
separator = os.path.sep
paths = []
names = []
files = []
# 遍历
for parent, dirname, filenames in os.walk(root_dir):
for filename in filenames:
if filename.endswith(filter):
paths.append(parent + separator)
names.append(filename)
for i in range(paths.__len__()):
files.append(paths[i] + names[i])
print (names.__len__().__str__() + " files have been found.")
paths.sort()
names.sort()
files.sort()
return paths, names, files
if __name__ == '__main__':
if sys.argv.__len__() >= 2:
if sys.argv[1] == 'help' or 'HELP':
print("Help Info:")
print("Resize Mode:[img_path] [start_x] [start_y] [x_range] [y_range]")
print("Join Mode:[img_dir] [img_type]")
print("Separate Mode:[img_path] [output_dir]")
pass
elif sys.argv[1] == 'resize' or 'RESIZE':
img_path = sys.argv[2]
start_x = int(sys.argv[3])
start_y = int(sys.argv[4])
x_range = int(sys.argv[5])
y_range = int(sys.argv[6])
bands_data = readTifImageWithWindow(img_path, start_x, start_y, x_range, y_range)
writeTif(bands_data, "cut_img.tif")
elif sys.argv[1] == 'join' or 'JOIN':
img_dir = sys.argv[2]
img_type = sys.argv[3]
paths, names, files = findAllFiles(img_dir, img_type)
bands_data = []
for i in range(files.__len__()):
band_data = readTifImage(files[i])
bands_data.extend(band_data)
print("joined " + (i + 1).__str__() + " bands.")
print(bands_data.__len__().__str__() + " bands in total.")
writeTif(bands_data, "join.tif")
elif sys.argv[1] == 'separate' or 'SEPARATE':
img_path = sys.argv[2]
output_dir = sys.argv[3]
bands_data = readTifImage(img_path)
for i in range(bands_data.__len__()):
writeTif([bands_data[i]], output_dir + os.path.sep + "band_" + i.__str__().zfill(2) + ".tif")
print("saved " + (i + 1).__str__() + "/" + bands_data.__len__().__str__())
else:
print("Unknown mode, input scriptname.py help to get help infomation.")
else:
print("Unknown mode, input scriptname.py help to get help infomation.")
以上便介绍了GDAL对遥感影像的读取和写入，并基于此实现了遥感影像的裁剪和波段融合/拆分。
本文作者原创,未经许可不得转载,谢谢配合
