python 矢量轉栅格

    雖然利用arcpy.PolygonToRaster_conversion可以将矢量檔案轉為栅格，但是arcpy的速度令人堪憂，有沒有更快更好的辦法呢？對于經常用python處理資料的人來說，會優先考慮用python來解決。百度一下“python 将矢量轉栅格”，會有别人分享的一些代碼，但是有些代碼我自己使用後，發現一些bug。于是在他們的基礎上，我對代碼進行了優化和修改，封裝成一個函數，經過自己資料的測試沒問題，分享給大家。

import gdal
import ogr
import gdalconst

def Rasterize(input_shp,input_tif,output_tif,field,filed_type,NoValue=-9999):
    """
    input_shp:需要轉為栅格的矢量檔案（矢量檔案路徑）
    input_tif:模闆栅格，用于讀取地理變換資訊、栅格大小，将其應用于新的栅格上
    output_tif:輸出栅格檔案（栅格檔案路徑）
    field:字元串，栅格值的字段
    filed_type:栅格值類型，一般選擇gdal.GDT_Int16,gdal.GDT_Int32,gdal.GDT_Float32,gdal.GDT_Float64等幾種類型
    NoValue:整型或浮點型，矢量空白區轉換後的值
    """
    data = gdal.Open(input_tif, gdalconst.GA_ReadOnly)
    geo_transform = data.GetGeoTransform()
    proj=data.GetProjection()
    
    open_shp = ogr.Open(input_shp)
    shp_ly = open_shp.GetLayer()
    x_min, x_max, y_min, y_max = shp_ly.GetExtent()
    
    pixel_size = geo_transform[1]
    x_res = int((x_max - x_min) / pixel_size)
    y_res = int((y_max - y_min) / pixel_size)

    target_ds = gdal.GetDriverByName('GTiff').Create(output_tif, x_res, y_res, 1, filed_type)

    target_ds.SetGeoTransform((x_min, pixel_size, 0.0, y_max, 0.0, -pixel_size))
    target_ds.SetProjection(proj)
    band = target_ds.GetRasterBand(1)
    band.SetNoDataValue(NoValue)
    band.FlushCache()
    
    if field is None:
        gdal.RasterizeLayer(target_ds, [1], shp_ly, None)
    else:
        OPTIONS = ['ATTRIBUTE=' + field]
        gdal.RasterizeLayer(target_ds, [1], shp_ly, options=OPTIONS)
    
    target_ds = None