命令行记录-矢量向栅格转换
作者:互联网
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、运行结果
标签:target,lyr,矢量,options,栅格,命令行,ogr,ds,gdal 来源: https://www.cnblogs.com/vividautumn/p/11557300.html