Python遥感影像数据读写

主要参考：https://blog.csdn.net/theonegis/article/details/80089375

1 影像基本信息读取

from osgeo import gdal

# 打开栅格影像
ds = gdal.Open('data\\Pic_Test.tif')

# =============================================================================
# 影像信息提取
# =============================================================================

# 投影
proj = ds.GetProjection()
 # 返回的是六个参数的tuple
trans = ds.GetGeoTransform() 
#波段数
nbands = ds.RasterCount
#行列数
nx = ds.RasterXSize
ny = ds.RasterYSize

# 获取数据集的元数据信息
metadata = ds.GetMetadata_Dict()
for key, value in metadata.items():
    print(f'{key} -> {value}')

b_Type = []
b_TypeName = []
b_NoData = []
b_MINMAX = []

for b in range(nbands):
    # 注意GDAL中的band计数是从1开始的
    band = ds.GetRasterBand(b + 1)
    # 波段数据的一些信息
    b_Type.append(band.DataType)  # DataType属性返回的是数字
    b_TypeName.append(gdal.GetDataTypeName(b_Type[b]))  #数据类型
    b_NoData.append(band.GetNoDataValue())  # NoData值
    b_MINMAX.append( band.ComputeRasterMinMax())  # 最大最小值

2 Dataset转为Numpy的ndarray

# 1 在数据集层面转换
image = ds.ReadAsArray()

print(f'数据的尺寸：{image.shape}')
# 输出结果为：数据的尺寸：(3, 4800, 4800)
# 这说明ReadAsArray方法将每个波段都转换为了一个二维数组

# 获得第一个波段的数据
band1 = image[0]

# 2 在波段层面的转换


for b in range(nbands):
    # 注意GDAL中的band计数是从1开始的
    band = ds.GetRasterBand(b + 1)
    band = band.ReadAsArray()
    print(f'波段大小：{band.shape}')

3 格式转换

# =============================================================================
# 格式转换
# GTiff -raster- (rw+vs): GeoTIFF
# HFA -raster- (rw+v): Erdas Imagine Images (.img)
# ENVI -raster- (rw+v): ENVI .hdr Labelled
# =============================================================================
    
# 方法一：Translate()函数
ds_env1 = gdal.Translate('data\\Pic_Test-env', ds, format='ENVI')

#方法二：CreateCopy()方法
driver = gdal.GetDriverByName('ENVI')
ds_env2 = driver.CreateCopy('data\\Pic_Test-env2', ds)



# 关闭数据集
ds = None
ds_env1 = None
ds_env2 = None

4 文件读写

4.1函数

from osgeo import gdal

#读图像文件
def read_img(filename):
    dataset=gdal.Open(filename)       #打开文件

    im_width = dataset.RasterXSize    #栅格矩阵的列数
    im_height = dataset.RasterYSize   #栅格矩阵的行数

    im_geotrans = dataset.GetGeoTransform()  #仿射矩阵
    im_proj = dataset.GetProjection() #地图投影信息
    im_data = dataset.ReadAsArray() #将数据写成数组，对应栅格矩阵

    del dataset 
    return im_proj,im_geotrans,im_data,im_width,im_height

#写文件，以写成tif为例
# =============================================================================
# 其他格式
# GTiff -raster- (rw+vs): GeoTIFF
# HFA -raster- (rw+v): Erdas Imagine Images (.img)
# ENVI -raster- (rw+v): ENVI .hdr Labelled
# =============================================================================
def write_img(filename,im_proj,im_geotrans,im_data):
    #gdal数据类型包括
    #gdal.GDT_Byte, 
    #gdal .GDT_UInt16, gdal.GDT_Int16, gdal.GDT_UInt32, gdal.GDT_Int32,
    #gdal.GDT_Float32, gdal.GDT_Float64

    #判断栅格数据的数据类型
    if 'int8' in im_data.dtype.name:
        datatype = gdal.GDT_Byte
    elif 'int16' in im_data.dtype.name:
        datatype = gdal.GDT_UInt16
    else:
        datatype = gdal.GDT_Float32

    #判读数组维数
    if len(im_data.shape) == 3:
        im_bands, im_height, im_width = im_data.shape
    else:
        im_bands, (im_height, im_width) = 1,im_data.shape 

    #创建文件
    driver = gdal.GetDriverByName("GTiff")            #数据类型必须有，因为要计算需要多大内存空间
    dataset = driver.Create(filename, im_width, im_height, im_bands, datatype)

    dataset.SetGeoTransform(im_geotrans)              #写入仿射变换参数
    dataset.SetProjection(im_proj)                    #写入投影

    if im_bands == 1:
        dataset.GetRasterBand(1).WriteArray(im_data)  #写入数组数据
    else:
        for i in range(im_bands):
            dataset.GetRasterBand(i+1).WriteArray(im_data[i])

    del dataset

4.2 读写影像

from osgeo import gdal
import numpy as np
import IMG_RW as ig

im_proj,im_geotrans,image,im_width,im_height = ig.read_img('data\\Pic_Test.tif')

#计算NDVI##############################################################
bnd_red = image[0].astype(float)  # 红波段
bnd_nir = image[3].astype(float)  # 近红外波段
idx_ndvi = (bnd_nir - bnd_red) / (bnd_nir + bnd_red)  # NDVI

#创建影像 法一：函数
out1_file = 'data\\NDVI.tif'
ig.write_img(out1_file,im_proj,im_geotrans,idx_ndvi)

#创建影像 法二：copy模板（out1_file）
out2_file = 'data\\NDVI-2.tif'
driver = gdal.GetDriverByName('GTiff')
dataset = driver.CreateCopy(out2_file, gdal.Open(out1_file))
dataset.GetRasterBand(1).WriteArray(idx_ndvi)  #写入数组数据


#波段组合#################################################################

image321 = np.array((image[2], image[1], image[0]),dtype = image.dtype)
#法一：
out3_file = 'data\\Pic_Test321.tif'
ig.write_img(out3_file,im_proj,im_geotrans,image321)

#法二：
out4_file = 'data\\Pic_Test321-2.tif'
driver = gdal.GetDriverByName('GTiff')
dataset2 = driver.CreateCopy(out4_file, gdal.Open(out3_file))
for i in range(3):
    dataset2.GetRasterBand(i+1).WriteArray(image321[i])