栅格化 GDAL 层
Posted
技术标签:
【中文标题】栅格化 GDAL 层【英文标题】:Rasterizing a GDAL layer 【发布时间】:2011-01-14 07:30:28 【问题描述】:编辑
这是正确的做法,documentation:
import random
from osgeo import gdal, ogr
RASTERIZE_COLOR_FIELD = "__color__"
def rasterize(pixel_size=25):
# Open the data source
orig_data_source = ogr.Open("test.shp")
# Make a copy of the layer's data source because we'll need to
# modify its attributes table
source_ds = ogr.GetDriverByName("Memory").CopyDataSource(
orig_data_source, "")
source_layer = source_ds.GetLayer(0)
source_srs = source_layer.GetSpatialRef()
x_min, x_max, y_min, y_max = source_layer.GetExtent()
# Create a field in the source layer to hold the features colors
field_def = ogr.FieldDefn(RASTERIZE_COLOR_FIELD, ogr.OFTReal)
source_layer.CreateField(field_def)
source_layer_def = source_layer.GetLayerDefn()
field_index = source_layer_def.GetFieldIndex(RASTERIZE_COLOR_FIELD)
# Generate random values for the color field (it's here that the value
# of the attribute should be used, but you get the idea)
for feature in source_layer:
feature.SetField(field_index, random.randint(0, 255))
source_layer.SetFeature(feature)
# Create the destination data source
x_res = int((x_max - x_min) / pixel_size)
y_res = int((y_max - y_min) / pixel_size)
target_ds = gdal.GetDriverByName('GTiff').Create('test.tif', x_res,
y_res, 3, gdal.GDT_Byte)
target_ds.SetGeoTransform((
x_min, pixel_size, 0,
y_max, 0, -pixel_size,
))
if source_srs:
# Make the target raster have the same projection as the source
target_ds.SetProjection(source_srs.ExportToWkt())
else:
# Source has no projection (needs GDAL >= 1.7.0 to work)
target_ds.SetProjection('LOCAL_CS["arbitrary"]')
# Rasterize
err = gdal.RasterizeLayer(target_ds, (3, 2, 1), source_layer,
burn_values=(0, 0, 0),
options=["ATTRIBUTE=%s" % RASTERIZE_COLOR_FIELD])
if err != 0:
raise Exception("error rasterizing layer: %s" % err)
原始问题
我正在寻找有关如何使用osgeo.gdal.RasterizeLayer()
的信息(文档字符串非常简洁,我在C 或C++ API 文档中找不到它。我只找到了java bindings 的文档)。
我改编了 unit test 并在由多边形组成的 .shp 上进行了尝试:
import os
import sys
from osgeo import gdal, gdalconst, ogr, osr
def rasterize():
# Create a raster to rasterize into.
target_ds = gdal.GetDriverByName('GTiff').Create('test.tif', 1280, 1024, 3,
gdal.GDT_Byte)
# Create a layer to rasterize from.
cutline_ds = ogr.Open("data.shp")
# Run the algorithm.
err = gdal.RasterizeLayer(target_ds, [3,2,1], cutline_ds.GetLayer(0),
burn_values=[200,220,240])
if err != 0:
print("error:", err)
if __name__ == '__main__':
rasterize()
它运行良好,但我得到的只是一个黑色的 .tif。
burn_values
参数是什么?可以使用RasterizeLayer()
栅格化具有不同颜色特征的图层吗?基于属性的值?
如果不能,我应该使用什么? AGG 是否适合渲染地理数据(我想要 no 抗锯齿和非常强大的渲染器,能够正确绘制非常大和非常小的特征,可能来自“脏数据”(退化多边形等) ..),有时以大坐标指定)?
这里,多边形是通过属性值来区分的(颜色无关紧要,我只想为属性的每个值设置不同的颜色)。
【问题讨论】:
感谢 Luper,今天这对我很有帮助! gdal 的文档很难找到正确的信息... 嗨@Luper,太好了!我正在寻找这个!您是否允许将您的示例代码(部分)包含在 GPLv3 许可的开源项目中,因为我正确地注明了您的姓名并链接到这个问题? @andreas-h 确定没问题。 @andreas-h 代码的最终形式可以在here 找到。它也是 GPLv3。 @LuperRouch 太好了,感谢您的链接! 【参考方案1】:编辑:我想我会使用 qGIS python 绑定:http://www.qgis.org/wiki/Python_Bindings
这是我能想到的最简单的方法。我记得以前用手卷过一些东西,但它很难看。 qGIS 会更容易,即使您必须进行单独的 Windows 安装(让 python 使用它)然后设置 XML-RPC 服务器以在单独的 python 进程中运行它。
我可以让 GDAL 正确光栅化,这也很棒。
我有一段时间没有使用 gdal,但这是我的猜测:
burn_values
用于假色,如果您不使用 Z 值。如果您使用burn=[1,2,3],burn_values=[255,0,0]
,则多边形内的所有内容都是[255,0,0]
(红色)。我不确定点会发生什么 - 他们可能不会绘图。
如果您想使用 Z 值,请使用 gdal.RasterizeLayer(ds,bands,layer,burn_values, options = ["BURN_VALUE_FROM=Z"])
。
我只是从你正在查看的测试中提取这个:http://svn.osgeo.org/gdal/trunk/autotest/alg/rasterize.py
另一种方法 - 将多边形对象拉出,并使用 shapely 绘制它们,这可能不吸引人。或者查看geodjango(我认为它使用openlayers 使用javascript 绘制到浏览器中)。
另外,您需要光栅化吗?如果您真的想要精确,导出 pdf 可能会更好。
实际上,我认为我发现使用 Matplotlib(在提取和投影特征之后)比光栅化更容易,而且我可以获得更多的控制权。
编辑:
这里有一个较低级别的方法:
http://svn.osgeo.org/gdal/trunk/gdal/swig/python/samples/gdal2grd.py\
最后,您可以迭代多边形(在将它们转换为局部投影之后),并直接绘制它们。但是你最好不要有复杂的多边形,否则你会有点悲伤。如果您有复杂的多边形...如果您想滚动自己的绘图仪,最好使用来自http://trac.gispython.org/lab 的 shapely 和 r-tree。
Geodjango 可能是一个询问的好地方.. 他们知道的比我多得多。他们有邮寄名单吗?周围也有很多 python 映射专家,但他们似乎都不担心这一点。我猜他们只是在 qGIS 或 GRASS 中绘制它。
说真的,我希望知道自己在做什么的人可以回复。
【讨论】:
是的,我需要栅格化(我编辑了问题以显示我想要做什么)。也许有像“BURN_VALUE_FROM_Z”这样的选项可以使用属性? 另外,我不明白为什么在我的测试中最终会出现黑色图像。 在#gdal 上的人们的帮助下,我已经能够使用 RasterizeLayer。问题出在转换/范围转移中,您必须使源范围和目标范围匹配,否则一切都会被剪掉。这实际上在文档中进行了解释,我不知道我是如何错过它的:D 无论如何感谢您的研究,我会接受您的回答并使用固定代码编辑我的问题。以上是关于栅格化 GDAL 层的主要内容,如果未能解决你的问题,请参考以下文章