Python GDAL矢量转栅格详解

前言：挺久没有更新博客了，前段时间课程实验中需要用代码将矢量数据转成栅格，常见的点栅格化方法通过计算将点坐标(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)

“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.RasterXSizey_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 1band = targetDataset.GetRasterBand(1)# 白色背景#NoData_value = -999NoData_value = 255band.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）：

全路径（ALL_TOUCHED=True）：

为了实现类似ArcGIS中矢量转栅格（每个要素赋不同值）的效果，可以先为用于测试的原始矢量数据添加一个属性字段用于赋值，例如下图的运动场面数据，以路网范围的栅格作为模板。

其属性表如下，pp_id为自行添加的字段，表示每个要素转换成栅格时的赋值：

函数调用：

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

得到的栅格化结果的每个要素具有不同的灰度值，对应pp_id字段的值。此时burn_values值会失效，根据pp_id字段值来对栅格赋值。

Python GDAL矢量转栅格详解相关推荐

python编程语法大全-Python编程入门——基础语法详解
今天小编给大家带来Python编程入门--基础语法详解. 关于怎么快速学python,可以加下小编的python学习群:611+530+101,不管你是小白还是大牛,小编我都欢迎,不定期分享干货每天 ...
python编程语法-Python编程入门——基础语法详解
今天小编给大家带来Python编程入门--基础语法详解. 一.基本概念 1.内置的变量类型: Python是有变量类型的,而且会强制检查变量类型.内置的变量类型有如下几种: #浮点 float_num ...
python编程if语法-Python编程入门基础语法详解经典
原标题:Python编程入门基础语法详解经典一.基本概念 1.内置的变量类型: Python是有变量类型的,而且会强制检查变量类型.内置的变量类型有如下几种: #浮点 float_number = ...
python编程语法-Python编程入门——基础语法详解（经典）
今天小编给大家带来Python编程入门--基础语法详解.温馨提示: 亮点在最后! 在这里还是要推荐下我自己建的Python开发学习群:301056051,群里都是学Python开发的,如果你正在学习P ...
用python绘制漂亮的图形-用python绘制图形的实例详解
1.环境系统:windows10 python版本:python3.6.1 使用的库:matplotlib,numpy 2.numpy库产生随机数几种方法import numpy as npnumpy ...
python怎么导入文件-Python文件如何引入？详解引入Python文件步骤
python基本语法--引入Python文件 1.新建python文件 :在同目录lib下创建mylib.py和loadlib.py两个文件 2.在mylib.py文件中创建一个Hello的类并且给 ...
python怎么画条形图-python绘制条形图方法代码详解
1.首先要绘制一个简单的条形图 import numpy as np import matplotlib.pyplot as plt from matplotlib import mlab from ...
python画直方图成绩分析-python plotly绘制直方图实例详解
计算数值出现的次数 import cufflinks as cf cf.go_offline() import numpy as np import pandas as pd set_slippage ...
python利器怎么编程-Python任务调度利器之APScheduler详解
任务调度应用场景所谓的任务调度是指安排任务的执行计划,即何时执行,怎么执行等.在现实项目中经常出现它们的身影:特别是数据类项目,比如实时统计每5分钟网站的访问量,就需要每5分钟定时从日志数据分析访问 ...

Python GDAL矢量转栅格详解

Python GDAL矢量转栅格详解相关推荐

最新文章

热门文章