命令行记录-矢量向栅格转换

1、

（1）重点学习如何应用 gdal.RasterizeLayer 函数

 

• gdal.RasterizeLayer( target_ds, [3, 2, 1],source_lyr, burn_values = [10,10, 55], options = ["ALL_TOUCHED=TRUE"])注意：栅格化时，将波段 1,2,3 上，分别赋值 55，10,10，source_lyr为打开的文件，burn_values和options这两个属性是非必须的。
  

 

（2）关键是选项的选择问题，有4种赋值方式，如：

 

• options=["ATTRIBUTE=LUCODE"]，表示字段" LUCODE "的属性值将为栅格值，如果没有设定字段名，则赋值为0；
• options=["BURN_VALUE_FROM=Z"]，表示 3D 图形的栅格值为其 Z 值（高程值） ；
• options=["ALL_TOUCHED=TRUE"]，表示图形接触到的像素均将输出；
• options=["MERGE_ALG=ADD","MERGE_ALG=REPLACE"]， 表示多图形覆盖同一像素的取值方式

2、
from osgeo import gdal, ogr, osr

#定义投影
sr = osr.SpatialReference('LOCAL_CS["arbitrary"]')
#在内存中，没有真正放到图层里去,创建一个 Shape 文件的图层，栅格值来自于CELSIUS字段
source_ds = ogr.GetDriverByName('Memory').CreateDataSource( 'wrk' )
mem_lyr = source_ds.CreateLayer('poly', srs=sr,geom_type=ogr.wkbPolygon )
mem_lyr.CreateField(ogr.FieldDefn('CELSIUS',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)):
    # 创建了一个空的feature
    feat = ogr.Feature(mem_lyr.GetLayerDefn())
    # 首先设置几何部分
    feat.SetGeometryDirectly(ogr.Geometry(wkt=wkt_geom[i]))#首先设置几何部分
    feat.SetField('CELSIUS', celsius_field_values[i])
    mem_lyr.CreateFeature(feat)

#在内存中，创建一个 100*100 大小的 3 波段的空白图像
target_ds = gdal.GetDriverByName('MEM').Create('', 100, 100, 1, gdal.GDT_Byte )
target_ds.SetGeoTransform( (1010,1,0,1060,0,-1) )
target_ds.SetProjection( sr.ExportToWkt())

#调用栅格化函数
err = gdal.RasterizeLayer( target_ds, [1], mem_lyr, options= ["ATTRIBUTE=CELSIUS"])

#将内存中的图像，存储到硬盘文件上
poly01_ts=gdal.GetDriverByName('GTiff').CreateCopy('./raster_poly01.tif',target_ds)
del target_ds
del source_ds

3、这段代码运行过程中需要注意的点是，在将内存中的栅格数据copy到硬盘后，需要及时退出解释器，使生成的句柄释放（执行了copy命令，同时内存数据还需要刷新到硬盘上）。解决方法是：

（1）演示完后，执行exit()退出解释器；

（2）在copy语句左侧加变量接收反馈的句柄，然后删除句柄。

4、运行结果