python 矢量數據轉柵格數據

from osgeo import gdal,osr,ogr
#定義投影
sr = osr.SpatialReference('LOCAL_CS["arbitrary"]')
#在內存中創建一個shape文件的圖層，含有兩個多邊形
source_ds = ogr.GetDriverByName('Memory').CreateDataSource( 'shapefile' )
source_lyr = source_ds.CreateLayer('poly', srs=sr, geom_type=ogr.wkbPolygon )
source_lyr.CreateField(ogr.FieldDefn('TCODE',ogr.OFTReal))
wkt_geom = ['POLYGON((1020 1030 40,1020 1045 30,1050 1045 20,1050 1030 35,1020 1030 40))',
 'POLYGON((1010 1046 85,1015 1055 35,1055 1060 26,1054 1048 35,1010 1046 85))']
#柵格值
celsius_field_values = [50,200]
for i in range(len(wkt_geom)):
    feat = ogr.Feature(source_lyr.GetLayerDefn())
    feat.SetGeometryDirectly(ogr.Geometry(wkt=wkt_geom[i]))
    feat.SetField('TCODE', celsius_field_values[i])
    source_lyr.CreateFeature(feat)
#在內存中，創建一個 100*100 大小的1波段的空白圖像
#‘’代表不往磁盤上寫的話，文件名可以是空
target_ds = gdal.GetDriverByName('MEM').Create('', 100, 100, 1, gdal.GDT_Byte )
target_ds.SetGeoTransform( (1000,1,0,1100,0,-1) )
target_ds.SetProjection( sr.ExportToWkt())
#調用柵格化函數。RasterizeLayer函數有四個參數，分別有柵格對象，波段，矢量對象，TCODE的屬性值將為柵格值
err = gdal.RasterizeLayer( target_ds, [1], source_lyr,options= ["ATTRIBUTE=TCODE"])
#將內存中的圖像，存儲到硬盤文件上
gdal.GetDriverByName('GTiff').CreateCopy('rasterized_poly.tif', target_ds)
del target_ds
del source_ds