Python GDAL学习笔记(一)

一、读取tif格式影像

导入GDAL模块并查看版本/位置

from osgeo import gdal
print("GDAL's version is:" + gdal.__version__)
print(gdal)

GDAL's version is:3.0.4
<module 'osgeo.gdal' from 'E:\\Anaconda3-2019\\lib\\site-packages\\osgeo\\gdal.py'>

打开影像

dataset = gdal.Open('F:\GDAL learning\Landsat8.tif', gdal.GA_ReadOnly)
print(dataset)

<osgeo.gdal.Dataset; proxy of <Swig Object of type 'GDALDatasetShadow *' at 0x0000026EA51416C0> >

查看行列号和波段数

num_bands = dataset.RasterCount
print('Number of bands in image: {n}\n'.format(n=num_bands))

rows = dataset.RasterYSize
cols = dataset.RasterXSize
print('Image size is : {r} rows x {c} colums\n'.format(r=rows,c=cols))

Number of bands in image: 2

Image size is : 8900 rows x 8563 colums

查看影像文件描述和其他属性

文件描述，元数据，影像驱动，投影，仿射变换参数

desc = dataset.GetDescription()
metadata = dataset.GetMetadata()
print('Raster Description : {desc}'.format(desc=desc))
print('Raster metadata ：')
print(metadata)
print('\n')

driver = dataset.GetDriver()
print('Raster Driver : {d}\n'.format(d=driver.ShortName))

proj = dataset.GetProjection()
print('Image projection : ')
print(proj + '\n')

gt = dataset.GetGeoTransform()
print('Image geo-transform : {gt}\n'.format(gt=gt))

Raster Description : F:\GDAL learning\Landsat8.tif
Raster metadata ：
{'AREA_OR_POINT': 'Area', 'TIFFTAG_XRESOLUTION': '1', 'TIFFTAG_YRESOLUTION': '1'}


Raster Driver : GTiff

Image projection : 
PROJCS["unnamed",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Albers_Conic_Equal_Area"],PARAMETER["latitude_of_center",0],PARAMETER["longitude_of_center",105],PARAMETER["standard_parallel_1",25],PARAMETER["standard_parallel_2",47],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]

Image geo-transform : (847767.0584, 30.0, 0.0, 4527946.9968, 0.0, -30.0)

打开数据集波段并获取统计信息

数据类型，基本统计值

red = dataset.GetRasterBand(1)
print(red)

datatype = red.DataType
print('Band datatype : {dt}'.format(dt=datatype))

datatype_name = gdal.GetDataTypeName(datatype)
print('Band datatypename : {dt}'.format(dt=datatype_name))

band_max, band_min, band_mean, band_stddev = red.GetStatistics(0,1)
print('Band range : {minmum} - {maxmum}'.format(maxmum=band_max,minmum=band_min))
print('Band mean, stddev:{m},{s}'.format(m=band_mean, s=band_stddev)）

<osgeo.gdal.Band; proxy of <Swig Object of type 'GDALRasterBandShadow *' at 0x0000026EA51545D0> >
Band datatype : 2
Band datatypename : UInt16
Band range : 55487.0 - 0.0
Band mean, stddev:4252.3194931158,3968.5033117154

利用numpy将波段数据写入数组

import numpy as np
print(np.__version__)

help(red.ReadAsArray)

red_data = red.ReadAsArray()
print(red_data)
print('Red band mean is {m}'.format(m=red_data.mean()))
print('size is : {sz}'.format(sz=red_data.shape))

Red band mean is 4252.319493115796
size is : (8900, 8563)

将波段数据写入数组的两种方法对比

image = np.zeros((dataset.RasterYSize, dataset.RasterXSize, dataset.RasterCount))
for b in range(dataset.RasterCount):
    band = dataset.GetRasterBand(b + 1)
    image[:, :, b] = band.ReadAsArray()

print(image.dtype)    

from osgeo import gdal_array
image_datatype = dataset.GetRasterBand(1).DataType

image_correct = np.zeros((dataset.RasterYSize, dataset.RasterXSize, dataset.RasterCount), 
                         dtype=gdal_array.GDALTypeCodeToNumericTypeCode(image_datatype))
for b in range(dataset.RasterCount):
    band = dataset.GetRasterBand(b + 1)
    image[:, :, b] = band.ReadAsArray()
print("比较数据类型：")
print('之前的读取方法：{dt}'.format(dt=image.dtype))
print('这种读取方法：{dt}'.format(dt=image_correct.dtype))

比较数据类型：
之前的读取方法：float64
这种读取方法：uint16

二、操作影像

计算NDVI

from osgeo import gdal, gdal_array
import osr
import numpy as np

dataset = gdal.Open('F:\GDAL learning\Landsat8.tif', gdal.GA_ReadOnly)
image_datatype = dataset.GetRasterBand(1).DataType
image = np.zeros((dataset.RasterYSize,dataset.RasterXSize,dataset.RasterCount),
                 dtype=gdal_array.GDALTypeCodeToNumericTypeCode(image_datatype))
for b in range(dataset.RasterCount):
    band = dataset.GetRasterBand(b+1)
    image[:, :, b] = band.ReadAsArray()

print('Red band mean : {r}'.format(r=image[:, :, 0].mean()))
print('NIR band mean : {nir}'.format(nir=image[:, :, 1].mean()))

ndvi = (image[:, :, 1] - image[:, :, 0]) / (
        image[:, :, 1] + image[:, :, 0]).astype(np.float64)

nan = np.isnan(ndvi)
ndvi = ndvi[~nan]

print('ndvi mean : {nm}'.format(nm=np.mean(ndvi)))
print('ndvi max :{nmax}'.format(nmax=np.percentile(ndvi, 95)))
print('ndvi min : {nmin}'.format(nmin=np.percentile(ndvi, 5)))

ndvi = (image[:, :, 1] - image[:, :, 0]) / (
        image[:, :, 1] + image[:, :, 0]).astype(np.float64)
ndvi[np.isnan(ndvi)] = 99

三、写出tif影像

driver = gdal.GetDriverByName( 'GTiff' )
ndvi_filename = 'F:\GDAL learning\Landsat8_ndvi.tif'
ndvi_ds = driver.Create( ndvi_filename, dataset.RasterXSize, 
                        dataset.RasterYSize, 1, gdal.GDT_Float64)

ndvi_ds.SetGeoTransform( dataset.GetGeoTransform() )
ndvi_ds.SetProjection( dataset.GetProjection() )
ndvi_ds.GetRasterBand(1).WriteArray( ndvi )

outBand = ndvi_ds.GetRasterBand(1)
outBand.FlushCache()
outBand.SetNoDataValue(99)
del ndvi_ds  #写出文件后一定要关闭文件，不然在python关闭之前文件为空。

四、运行结果

PS：没有进行定标和大气校正