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)