目錄

• • • 一、讀取tif格式影像
    • • 導入GDAL模塊并查看版本/位置
      • 打開影像
      • 查看行列號和波段數
      • 查看影像文件描述和其他屬性
      • 打開數據集波段并獲取統計信息
      • 利用numpy將波段數據寫入數組
      • 將波段數據寫入數組的兩種方法對比
    • 二、操作影像
    • • 計算NDVI
    • 三、寫出tif影像
    • 四、運行結果

一、讀取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: 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 : ze8trgl8bvbq\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 : 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.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).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_array import osr import numpy as npdataset = 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：沒有進行定標和大氣校正

總結

以上是生活随笔為你收集整理的Python GDAL学习笔记(一)的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。