osgeo：GDAL

from osgeo import gdal, gdalconst, osr, ogr
from pathlib import Path

一、GDAL驱动

• 系统支持的数据格式

# gdal.GetDriverByName('GTiff')
drv_count = gdal.GetDriverCount()
for idx in range(drv_count):
    driver = gdal.GetDriver(idx)
    # print("%10s: %s" % (driver.ShortName, driver.LongName))
# GTiff: GeoTIFF; ESRI Shapefile: ESRI Shapefile; ENVI: ENVI .hdr Labelled
# SAFE: Sentinel-1 SAR SAFE Product

# 访问索引图像信息
!gdalinfo ./GDAL_data/Test.tif

二、读取遥感影像的信息

• dataset.GetDescription() 获得栅格的描述信息
• dataset.RasterCount 获得栅格数据集的波段数
• dataset.RasterXSize 栅格数据的宽度(X方向上的像素个数)
• dataset.RasterYSize 栅格数据的高度(Y方向上的像素个数)
• dataset.GetGeoTransform() 栅格数据的六参数。
• dataset.GetProjection() 栅格数据的投影, 用wkt表示
  

2.1 直接读取

# 读图像文件
def read_img(filename):
    dataset = gdal.Open(filename)  # 打开文件
    # print(dataset.GetDescription())
    im_width = dataset.RasterXSize  # 栅格矩阵的列数
    im_height = dataset.RasterYSize  # 栅格矩阵的行数
    im_bands = dataset.RasterCount  # 波段数
    im_geotrans = dataset.GetGeoTransform()  # 仿射矩阵，左上角像素的大地坐标和像素分辨率
    # （左上角x坐标0， 水平分辨率1，旋转参数2，左上角y坐标3，旋转参数4，垂直分辨率5）
    #  根据GDAL的六参数模型将给定的投影或地理坐标转为影像图上坐标（行列号）
    #  Xgeo = GT0 + Xpixel*GT1 + Ypixel*GT2
    #  Ygeo = GT3 + Xpixel*GT4 + Ypixel*GT5
    im_proj = dataset.GetProjection()  # 地图投影信息，字符串表示
    im_data = np.zeros((im_bands, im_height, im_width))
    for i in range(im_bands):
        im_data[i, :, :] = dataset.GetRasterBand(i+1).ReadAsArray(0, 0, im_width, im_height)

    del dataset  # 关闭对象dataset，释放内存
    return im_proj, im_geotrans, im_data, im_width, im_height, im_bands

im_proj, im_geotrans, im_data, im_width, im_height, im_bands = read_img('./GDAL_data/Test.tif')
im_geotrans, (im_width, im_height, im_bands)

((375870.0, 10.0, 0.0, 4480780.0, 0.0, -10.0), (2923, 2924, 7))

2.2 分块读取

# 分块读取影像
# -*- coding: utf-8 -*-
"""
真彩色生成
@author: LCH
"""
from osgeo import gdal
import numpy as np
import os
import math
import tqdm
import time

def natural_color(inputpath, outpath):
    datatset = gdal.Open(inputpath)
    xsize, ysize = datatset.RasterXSize, datatset.RasterYSize

    oneband = datatset.GetRasterBand(1)
    dataTYpe = gdal.GetDataTypeName(oneband.DataType)

    band_count = 3
    # 判断栅格数据的数据类型
    if 'int8' == dataTYpe:
        datatype = gdal.GDT_Byte
    elif 'int16' == dataTYpe:
        datatype = gdal.GDT_UInt16
    else:
        datatype= gdal.GDT_Float32

    # 驱动器
    driver = gdal.GetDriverByName("GTiff")
    out_tif = driver.Create(outpath, xsize, ysize, band_count, datatype)
    out_tif.SetProjection(datatset.GetProjection())
    out_tif.SetGeoTransform(datatset.GetGeoTransform())  # 地理坐标变换

    # 分块运算
    nBlockSize = 2048
    i, j = 0, 0
    b = xsize * ysize
    # 进度条参数
    XBlockcount = math.ceil(xsize / nBlockSize)
    YBlockcount = math.ceil(ysize / nBlockSize)

    try:
        with tqdm.tqdm(total=XBlockcount * YBlockcount, iterable='iterable',desc = '按整幅影像分块读取', mininterval=10) as pbar:
            while i < ysize:
                while j < xsize:
                    # 保存分块大小
                    nXBK, nYBK = nBlockSize, nBlockSize

                    # 最后不够分块的区域，有多少读取多少
                    if i + nBlockSize > ysize:
                        nYBK = ysize - i
                    if j + nBlockSize > xsize:
                        nXBK = xsize - j

                    # 分块读取影像
                    Image = datatset.ReadAsArray(j, i, nXBK, nYBK)
                    Image_3 = Image[:3,:,:]
                    out_tif.WriteArray(Image_3 , j, i)

                    j = j + nXBK
                    time.sleep(1)
                    pbar.update(1)
                j = 0
                i = i + nYBK
    except KeyboardInterrupt:
        pbar.close()
        raise
    pbar.close()
    out_tif.BuildOverviews('average', [2, 4, 8, 16, 32, 64])
    out_tif.FlushCache()  # 写入硬盘
    out_tif = None  # 关闭tif文件
    print('转换成功')


if __name__ == "__main__":
    file_name = r"G:\提取地表水\GF-1\GF1_WFV2_E114.1_N22.6_20210727_L1A0005780095_ortho.tif"
    output_path = r"G:\提取地表水\GF-1\GF1_WFV2_E114.1_N22.6_20210727_L1A0005780095_ortho_RGB1.tif"
    natural_color(file_name , output_path)

2.3 其他操作

ds = gdal.Open('./GDAL_data/Test.tif')
dir(ds)[:3] + ['... ...'] + dir(ds)[-3:]

# 获取栅格数据下的子数据集：HDF, NetCDF, SAFE等格式
ds = gdal.Open('./GDAL_data/Test.safe').GetSubDatasets()

gdal.GetDataTypeName(3), gdal.GetDataTypeByName('Int16')

(‘Int16’, 3)

三、获取波段信息及读取波段数据

• band.ReadAsArray() # 读取图像数据(以数组的形式)
• band.ReadRaster() # 读取图像数据(以二进制的形式)
  • xoff：列读取的起点，默认值为0;
  • yoff：行读取的起点，默认值为0;
  • win_xsize：要读取的列数，默认读取所有列 ;
  • win_ysize：要读取的行数，默认读取所有行 ;
  • buf_xsize：输出数组里的列数，默认为win_xsize的值，如果值不同于win_xsize，数据将会重采样 ;
  • buf_ysize：输出数组里的行数，默认为win_xsize的值，如果值不同于win_ysize，数据将会重采样 ;
  • buf_obj：用于存放数组，而不是创建数组的NumPy数组。
• WriteArray(array, xoff=0, yoff=0, band_list=None, interleave=‘band’, resample_alg=0, callback=None, callback_data=None)

band = ds.GetRasterBand(1)
# dir(band)[:6] + ['... ...'] + dir(band)[-3:]
# 获取波段大小
band.XSize, band.YSize, band.DataType, gdal.GetDataTypeName(band.DataType)

(2923, 2924, 2, ‘UInt16’)

# 获取波段数据的属性
band.GetNoDataValue(), band.GetMaximum(), band.GetMinimum()

(None, None, None)

band.ComputeRasterMinMax()

(0.0, 10088.0)

band.ReadAsArray(100, 100, 5, 5, 5, 5)

array([[1234, 1244, 1192, 1446, 2060],
[1232, 1207, 1198, 1393, 1961],
[1248, 1214, 1196, 1286, 1543],
[1251, 1248, 1225, 1312, 1518],
[1244, 1244, 1221, 1296, 1550]], dtype=uint16)

band.ReadRaster(100, 100, 5, 5, 5, 5)

bytearray(b'\xd2\x04\xdc\x04\xa8\x04\xa6\x05\x0c\x08\xd0\x04\xb7\x04\xae\x04q\x05\xa9\x07\xe0\x04\xbe\x04\xac\x04\x06\x05\x07\x06\xe3\x04\xe0\x04\xc9\x04 \x05\xee\x05\xdc\x04\xdc\x04\xc5\x04\x10\x05\x0e\x06')

# 指数运算时，避免分母为0
import numpy as np
data1 = ds.GetRasterBand(3).ReadAsArray(100, 100, 10, 10)
data2 = ds.GetRasterBand(2).ReadAsArray(100, 100, 10, 10)
mask = np.greater(data1 + data2, 0)  # np.where(data1 + data2 == 0, -99, (data2 - data1) / (data2 + data1))
# nodata = np.full((Green_arr.shape[0], Green_arr.shape[1]), -2, dtype=np.float32)
# ndwi = np.true_divide(numerator, denominator, out=nodata, where=denominator != 0.0)
express = np.choose(mask, (-99, (data2 - data1) / (data2 + data1)))

四、创建影像

• driver.Create(outName, width, height, bandNum, DataType)
  • outName.setGeoTransform()
  • outName.setProjection()
• driver.CreateCopy(outName, src_ds, strict=0, TILED=YES, COMPRESS=PACKBITS)
  • src_ds = gdal.Open(src_fileName)

driver = gdal.GetDriverByName('GTiff')
dst_ds = driver.Create()
dst_ds.SetGeoTransform([444720, 30, 0, 3751320, 0, -30])
# 美国常用的空间参考：NAD27
srs = osr.SpatialReference()
srs.SetUTM(11, 1)
srs.SetWellKnownGeogCS("NAD27")
dst_ds.SetProjection(srs.ExportToWkt())
raster = np.zeros((512, 512), dtype=np.uint8)
dst_ds.GetRasterBand(1).WriteArray(raster)
dst_ds = None

src_filename = "./GDAL_data/Test.tif"
dst_filename = "./GDAL_data/Test_RGB.tif"
driver = gdal.GetDriverByName('GTiff')
src_ds = gdal.Open(src_filename)
dst_ds = driver.CreateCopy(dst_filename, src_ds, 0)

def Write_Tiff(img_arr, geomatrix, projection, path):
    #     img_bands, img_height, img_width = img_arr.shape
    if "int8" in img_arr.dtype.name:
        datatype = gdal.GDT_Byte
    elif "int16" in img_arr.dtype.name:
        datatype = gdal.GDT_UInt16
    else:
        datatype = gdal.GDT_Float32

    if len(img_arr.shape) == 3:  # 多波段写入
        img_bands, img_height, img_width = img_arr.shape
        driver = gdal.GetDriverByName("GTiff")
        dataset = driver.Create(path, int(img_width), int(img_height), int(img_bands), datatype)
        
        if (dataset is not None) and (geomatrix != "") and (projection != ""):
            dataset.SetGeoTransform(geomatrix)  # 写入仿射变换参数
            dataset.SetProjection(projection)  # 写入投影
        for i in range(img_bands):
            dataset.GetRasterBand(i + 1).WriteArray(img_arr[i])
        del dataset

    elif len(img_arr.shape) == 2:  # 单波段写入
        # img_arr = np.array([img_arr])
        img_height, img_width = img_arr.shape
        img_bands = 1
        # 创建文件
        driver = gdal.GetDriverByName("GTiff")
        dataset = driver.Create(path, int(img_width), int(img_height), int(img_bands), datatype)

        if dataset is not None:
            dataset.SetGeoTransform(geomatrix)  # 写入仿射变换参数
            dataset.SetProjection(projection)  # 写入投影
        dataset.GetRasterBand(1).WriteArray(img_arr)
        del dataset

五、Warp函数（重投影、镶嵌、重采样等）

• WarpOptions(options=None, format=None,
  • outputBounds：输出边界(minX, minY, maxX, maxY)
  • outputBoundsSRS：空间参考系统
  • xRes, yRes：输出分辨率,可用于重采样 $\Rightarrow$ 自动计算矩阵的大小
  • width, height：栅格宽高,可用于重采样 $\Rightarrow$ 自动计算像元大小
  • srcSRS, dstSRS：源和输出SRS，可用于重投影、镶嵌等
  • outputType, workingType: 输出类型/工作类型(GDT_Byte等)
  • resampleAlg:重采样模式(near/bilinear/cubic/average)
  • transformerOptions: 转换参数
• Warp：用于坐标系转换、投影变换、图像合并与镶嵌、地理范围裁剪、更改分辨率、矢量裁剪等方面（扭曲一个或多个数据集），关键的参数在于options
  • destNameOrDestDS — Output dataset name or object
  • srcDSOrSrcDSTab — an array of Dataset objects or filenames, or a Dataset object or a filename
  • options — gdal.WarpOptions()

# 重投影
ds1 = gdal.Open('./GDAL_data/Test1.tif')
ds2 = gdal.Open('./GDAL_data/Test2.tif')
rojectedtmp = "./GDAL_data/Test_projected.tif"
ds = gdal.Warp(rojectedtmp, ds1, dstSRS="EPSG:4326")  # epsg可以通过https://epsg.io/查询，这里给出的是wgs84
# 镶嵌
outpath = './GDAL_data/Test_mosaic.tif'
options = gdal.WarpOptions(srcSRS=ds1, dstSRS=ds2, format='GTiff')
gdal.Warp(outpath, [ds1, ds2], options=options)
# 重采样
gdal.Warp(outpath, ds1, xRes=20, yRes=20)

# 批量镶嵌
def imgs_mosaic(inputdir, outdir):
    images = sorted(Path(inoutdir).glob('*.tif'))
    img_list = []
    for i in images:
        filename = Path(inoutdir) / i
        img_list = img_list.append(filename)
    for img1 in img_list:
        out_file1 = img1.split('_')[-5]
        for img2 in img_list:
            out_file2 = img2.split('_')[-5]
            if out_file1 == out_file2:
                if img1 != img2:
                    input_img1 = gdal.Open(img1, gdal.GA_ReadOnly)
                    inputProj = input_img1.GetProjection()
                    input_img2 = gdal.Open(img2, gdal.GA_ReadOnly)
                    out_file_name = outdir + '/' + out_file2 + '.tif'
                    if not Path(out_file_name).exists():
                        options = gdal.WarpOptions(srcSRS=inputProj, dstSRS=inputProj, format='GTiff')
                        gdal.Warp(out_file_name, [input_img1, input_img2], options=options)
                break
    print('已全部完成！')

六、GDAL和 Pillow 的互操作

6.1 数据读取

• gdal.Open(‘./GDAL_data/Test.tif’)
• PIL.Image.open(‘./GDAL_data/Test.tif’)

from PIL import Image
ds = gdal.Open('./GDAL_data/Test_RGB2.tif')
ds2 = Image.open('./GDAL_data/Test_RGB2.tif')
ds, ds2

6.2 数据转换

读取范围：

• GDAL（顶点X、顶点Y、宽、高）
• Pillow（顶点X、顶点Y、终点X，终点Y）

data = ds.ReadAsArray(30, 70, 5, 5)
list1 = [i for i in data]
datas = [np.reshape(i, (-1, 1)) for i in data]
datas = np.concatenate(datas, 1)
datas.tobytes(), ds2.crop((30, 70, 35, 75)).tobytes()

# 从波段看
r, g, b = ds2.crop((30, 70, 35, 75)).split()
band = ds.GetRasterBand(1)
r.tobytes(), band.ReadRaster(30, 70, 5, 5)