目录 一、读取tif格式影像导入GDAL模块并查看版本/位置打开影像查看行列号和波段数查看影像文件描述和其他属性打开数据集波段并获取统计信息利用numpy将波段数据写入数组将波段数据写入数组的两种方法对比 二、操作影像计算NDVI 三、写出tif影像四、运行结果
一、读取tif格式影像 导入GDAL模块并查看版本/位置 from osgeo import gdalprint("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 learningLandsat8.tif', gdal.GA_ReadOnly)print(dataset) <osgeo.gdal.Dataset; proxy of <Swig Object of type 'GDALDatasetShadow *' at 0x0000026EA51416C0> > 查看行列号和波段数 num_bands = dataset.RasterCountprint('Number of bands in image: {n}n'.format(n=num_bands))rows = dataset.RasterYSizecols = dataset.RasterXSizeprint('Image size is : {r} rows x {c} columsn'.format(r=rows,c=cols)) Number of bands in image: 2Image 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 learningLandsat8.tifRaster metadata ：{'AREA_OR_POINT': 'Area', 'TIFFTAG_XRESOLUTION': '1', 'TIFFTAG_YRESOLUTION': '1'}Raster Driver : GTiffImage 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.DataTypeprint('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 : 2Band datatypename : UInt16Band range : 55487.0 - 0.0Band mean, stddev:4252.3194931158,3968.5033117154 利用numpy将波段数据写入数组 import numpy as npprint(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.319493115796size 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_arrayimage_datatype = dataset.GetRasterBand(1).DataTypeimage_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_arrayimport osrimport numpy as npdataset = gdal.Open('F:GDAL learningLandsat8.tif', gdal.GA_ReadOnly)image_datatype = dataset.GetRasterBand(1).DataTypeimage = 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 learningLandsat8_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：没有进行定标和大气校正