GDAL的基本操作

　　上一節簡單介紹了GDAL，這一節將介紹一些GDAL的基本操作，如影像讀寫、波段提取、波段合成等。代碼均用python編寫。

1.遙感影像的讀寫

　　眾所周知，遙感影像是以柵格形式存儲的，GDAL中使用dataset表示一個柵格數據（使用抽象類GDALDataset表示），一個dataset包含了對於柵格數據的波段，空間參考以及元數據等信息。一張GeoTIFF遙感影像，一張DEM影像，或者一張土地利用圖，在GDAL中都是一個GDALDataset。

• 坐標系統（使用OGC WKT格式表示的空間坐標系統或者投影系統）
• 地理放射變換（使用放射變換表示圖上坐標和地理坐標的關系）
• GCPs（大地控制點記錄了圖上點及其大地坐標的關系，通過多個大地控制點可以重建圖上坐標和地理坐標的關系）
• 元數據（鍵值對的集合，用於記錄和影像相關的元數據信息）
• 柵格波段（使用GDALRasterBand類表示，真正用於存儲影像柵格值，一個柵格數據可以有多個波段）
• 顏色表（Color Table用於圖像顯示）

　　GDAL使用Open（）函數讀取影像數據，函數返回Dataset對象。GDAL目前支持約100種格式的柵格數據讀取，包括ERDAS Imagine、ENVI、GRASS、GeoTIFF、HDF4、HDF5、TIFF、JPEG、JPEG2000、PNG、GIF、BMP等。

import numpy as np
from osgeo import gdal
in_ds = gdal.Open(r"D:\data\test.tif")  # 打開樣本文件
xsize = in_ds.RasterXSize  # 獲取行列數
ysize = in_ds.RasterYSize
bands = in_ds.RasterCount  # 獲取波段數
geotransform = in_ds.GetGeoTransform() # 獲取投影信息
projection = in_ds.GetProjectionRef()
block_data = in_ds.ReadAsArray(0,0,xsize,ysize).astype(np.float32)# 獲取影像信息
B = block_data[0, :, :] #獲取第一波段數據
G = block_data[1,:, :] #獲取第二波段數據
R = block_data[2,:, :] #獲取第三波段數據

　　Dataset對象的RasterYSize、RasterXSize和RasterCount屬性分別返回柵格數據的行數、列數和波段數。Dataset對象的GetGeoTransform()方法返回柵格數據的坐標轉換參數，即行列坐標與空間坐標的轉換參數。Dataset對象的GetProjection()方法返回柵格數據的空間參照系統信息（WKT文本）。柵格數據通常有多個波段， Dataset對象的GetRasterBand(index)方法將返回某個波段的數據對象（Band對象），如ds.GetRasterBand(1)返回第一個波段的數據對象（注意：第一波段是1，不是0GetMinimum()和GetMaximum() 方法返回波段數據的最小值和最大值。ComputeStatistics()方法返回波段數據的統計信息，包括最小值、最大值、平均值和標准偏差，approx_ok參數表示是否抽樣統計，值為False表示統計所有柵格。ReadAsArray()函數返回的二維數組，astype()可以改變數組的數據類型。

driver = gdal.GetDriverByName('GTiff')
out_dataset=driver.Create("RGB.tif",xsize,ysize,3,gdal.GDT_Float32)
out_band1 = out_dataset.GetRasterBand(1)
out_band1.WriteArray(B)
out_band2 = out_dataset.GetRasterBand(2)
out_band2.WriteArray(G)
out_band3 = out_dataset.GetRasterBand(3)
out_band3.WriteArray(R)
out_dataset.SetGeoTransform(geotransform)  # 寫入仿射變換
out_dataset.SetProjection(projection)

　　將上一步得到的R、G、B三個波段數據存儲為新的數據（RGB.tif）。gdal.GDT_Float32代表數據類型，數據類型決定了柵格值的范圍，如數據類型為GDT_Byte，則柵格值的范圍是0～255；如要存儲的數據柵格值的范圍是0～65535，則數據類型應該是GDT_UInt16。　　

　　gdal常用的數據類型包括：

• • gdal.GDT_Byte           8bit正整型
  • gdal.GDT_UInt16       16bit正整型
  • gdal.GDT_Int16          16bit整型
  • gdal.GDT_UInt32        32bit 正整型
  • gdal.GDT_Int32           32bit整型 
  • gdal.GDT_Float32       32bit 浮點型
  • gdal.GDT_Float64       64bit 浮點型 

2.柵格數據行列號和地理坐標相互轉換

2.1行列坐標轉空間坐標

from osgeo import gdal
def Pixel2world(geotransform, line, column):
    originX = geotransform[0]
    originY = geotransform[3]
    pixelWidth = geotransform[1]
    pixelHeight = geotransform[5]
    x = column*pixelWidth + originX - pixelWidth/2
    y = line*pixelHeight + originY - pixelHeight/2
    return(x,y)

　　line和column為行列坐標（行列號）；pixelWidth和pixelHeigt為x方向比例尺和y方向比例尺；originX和originY為左上角的x和y坐標。

2.2空間坐標轉行列坐標：

from osgeo import gdal
def world2Pixel(geotransform, x, y):
    originX = geotransform[0]
    originY = geotransform[3]
    pixelWidth = geotransform[1]
    pixelHeight = geotransform[5]
    line = int((y-originY)/pixelHeight)+1
    column = int((x-originX)/pixelWidth)+1
    return (line,column)

 

import numpy as np