Python GDAL矢量轉柵格詳解

本文轉載自查看原文 2022-04-20 23:25 1133

前言：挺久沒有更新博客了，前段時間課程實驗中需要用代碼將矢量數據轉成柵格，常見的點柵格化方法通過計算將點坐標(X，Y)轉換到格網坐標(I，J)，線柵格化方法主要有DDA算法、Bresenham算法等，根據實現效果也可分為八方向和全路徑柵格化方法等，面柵格化方法主要有種子點填充、掃面線算法、邊界代數法等。詳細算法實現可參考GIS中將矢量數據轉換柵格數據算法和 GIS算法基礎（五）矢量數據向柵格數據的轉換（點，線算法實現）這兩篇博客。GDAL為用戶提供了矢柵轉換的方法，但網絡上相關資料比較少，官方文檔不明所以，因此本文將基於Python GDAL談談矢量轉柵格的用法。

使用Python GDAL的RasterizeLayer()方法可以比較容易地實現矢量圖層轉換為柵格，先看看官方API文檔給出的參數列表。

比較重要的參數有：

dataset 即輸出的柵格數據集
bands 輸出波段
layer 輸入待轉換的矢量圖層
pfnTransformer 幾何圖形坐標轉換圖像行列號函數
pTransformArg 幾何圖形坐標轉換圖像行列號參數
burn_values 設置每個輸出波段的像素值，需要輸入列表形式，如[0]
GDALProgressFunc callback 進度條回調函數

options 可選的柵格化參數（字符串參數對），對於該參數文檔官方文檔沒有詳細說明，可以參考C++ GDAL的類似API，如下：

"ATTRIBUTE":

指定輸入矢量數據的屬性字段中的字段值作為柵格值寫入柵格文件中，該值將輸出到所有輸出波段中。假如該值指定了，burn_Values參數的值將失效數可以設置為空。

"CHUNKYSIZE":

指定該運行操作的塊的高度。該值越大所需的計算時間越小。如果該值沒有設置或者設置為0則由GDAL的緩存大小根據公式：緩存所占的字節數/掃描函數的字節數得到。所以該值不會超出緩存的大小。

"ALL_TOUCHED":

設置為TRUE表示所有的像素接觸到矢量中線或多邊形，否則只是多邊形中心或被Bresenham算法選中的部分。默認值為FALSE。簡單來說，FALSE是八方向柵格化，TRUE是全路徑柵格化。

八方向柵格化：（圖片來自https://malagis.com/gis-vector-grid-data-conversion-algorithm.html）
全路徑柵格化：（圖片來自https://malagis.com/gis-vector-grid-data-conversion-algorithm.html)

2014-02-20_232654

"BURN_VALUE_FROM":
```
用於設置幾何圖形的Z值
```

"MERGE_ALG":

設置重寫或增加新值到柵格數據中。選擇REPLACE為 重寫，選擇ADD為添加一個新值到已存在的柵格數據中。默認值為REPLACE。

矢量轉柵格函數如下：

from osgeo import gdal,ogr,osr
from gdal import gdalconst

# 需要注意field，all_touch這些option的值必須為字符串
def vector2raster(inputfilePath, outputfile, templatefile,bands=[1],burn_values=[0],field="",all_touch="False"):
    # 輸入矢量文件
    inputfilePath = inputfilePath
    # 輸出柵格文件
    outputfile = outputfile
    # 柵格模板文件，確定輸出柵格的元數據（坐標系等，柵格大小，范圍等）
    templatefile = templatefile
    # 打開柵格模板文件
    data = gdal.Open(templatefile, gdalconst.GA_ReadOnly)
    # 確定柵格大小
    x_res = data.RasterXSize
    y_res = data.RasterYSize
    # 打開矢量文件
    vector = ogr.Open(inputfilePath)
    # 獲取矢量圖層
    layer = vector.GetLayer()
    # 查看要素數量
    featureCount = layer.GetFeatureCount()
    # print(featureCount)

    # 創建輸出的TIFF柵格文件
    targetDataset = gdal.GetDriverByName('GTiff').Create(outputfile, x_res, y_res, 1, gdal.GDT_Byte)
    # 設置柵格坐標系與投影
    targetDataset.SetGeoTransform(data.GetGeoTransform())
    targetDataset.SetProjection(data.GetProjection())
    # 目標band 1
    band = targetDataset.GetRasterBand(1)
    # 白色背景
    #NoData_value = -999
    NoData_value = 255
    band.SetNoDataValue(NoData_value)
    band.FlushCache()
    if field:
        # 調用柵格化函數。RasterizeLayer函數有四個參數，分別有柵格對象，波段，矢量對象，options
        # options可以有多個屬性，其中ATTRIBUTE屬性將矢量圖層的某字段屬性值作為轉換后的柵格值
        gdal.RasterizeLayer(targetDataset, bands,layer, burn_values=burn_values,options=["ALL_TOUCHED="+all_touch,"ATTRIBUTE="+field])
    else:
        gdal.RasterizeLayer(targetDataset, bands,layer, burn_values=burn_values,options=["ALL_TOUCHED="+all_touch])

對於ALL_TOUCHED屬性實現的八方向和全路徑柵格化，可以從下面兩幅圖看到結果略有不同，注意此時的ALL_TOUCHED屬性設置必須全為字符串。

函數調用：

vector2raster(inputpath,outputpath,maskpath,field="",all_touch="False")

八方向（ALL_TOUCHED=False）：

在這里插入圖片描述