python 栅格化 GDAL 层
声明：本页面是StackOverFlow热门问题的中英对照翻译，遵循CC BY-SA 4.0协议，如果您需要使用它，必须同样遵循CC BY-SA许可，注明原文地址和作者信息，同时你必须将它归于原作者（不是我）：StackOverFlow
原文地址: http://stackoverflow.com/questions/2220749/
Warning: these are provided under cc-by-sa 4.0 license. You are free to use/share it, But you must attribute it to the original authors (not me):
StackOverFlow
Rasterizing a GDAL layer
提问by Luper Rouch
Edit
编辑
Here is the proper way to do it, and the 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)
Original question
原始问题
I'm looking for information on how to use osgeo.gdal.RasterizeLayer()(the docstring is very succinct, and I can't find it in the C or C++ API docs. I only found a doc for the java bindings).
我正在寻找有关如何使用的信息osgeo.gdal.RasterizeLayer()（文档字符串非常简洁，我在 C 或 C++ API 文档中找不到它。我只找到了java bindings的文档）。
I adapted a unit testand tried it on a .shp made of polygons:
我改编了一个单元测试并在由多边形组成的 .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()
It runs fine, but all I obtain is a black .tif.
它运行良好，但我得到的只是一个黑色的 .tif。
What's the burn_valuesparameter for ? Can RasterizeLayer()be used to rasterize a layer with features colored differently based on the value of an attribute ?
什么是burn_values对的参数？可RasterizeLayer()用于栅格化具有基于属性值的不同颜色要素的图层？
If it can't, what should I use ? Is AGGsuitable for rendering geographic data (I want noantialiasing and a very robust renderer, able to draw very large and very small features correctly, possibly from "dirty data" (degenerate polygons, etc...), and sometimes specified in large coordinates) ?
如果不能，我应该用什么？是AGG适合绘制地理数据（我想没有抗锯齿和一个非常强大的渲染器，能够正确地绘制非常大和非常小的特征，可能是由“脏数据”（退化多边形，等...），有时在指定的大坐标）？
Here, the polygons are differentiated by the value of an attribute (the colors don't matter, I just want to have a different one for each value of the attribute).
在这里，多边形由属性值区分（颜色无关紧要，我只想为属性的每个值设置不同的颜色）。
采纳答案by wisty
EDIT: I guess I'd use qGIS python bindings: http://www.qgis.org/wiki/Python_Bindings
编辑：我想我会使用 qGIS python 绑定：http://www.qgis.org/wiki/Python_Bindings
That's the easiest way I can think of. I remember hand rolling something before, but it's ugly. qGIS would be easier, even if you had to make a separate Windows installation (to get python to work with it) then set up an XML-RPC server to run it in a separate python process.
这是我能想到的最简单的方法。我记得以前手滚过东西，但很丑。qGIS 会更容易，即使您必须进行单独的 Windows 安装（让 python 使用它）然后设置一个 XML-RPC 服务器以在单独的 python 进程中运行它。
I you can get GDAL to rasterize properly that's great too.
我可以让 GDAL 正确光栅化，这也很棒。
I haven't used gdal for a while, but here's my guess:
我有一段时间没有使用 gdal，但这是我的猜测：
burn_valuesis for false color if you don't use Z-values. Everything inside your polygon is [255,0,0](red) if you use burn=[1,2,3],burn_values=[255,0,0]. I'm not sure what happens to points - they might not plot.
burn_values如果不使用 Z 值，则用于假色。您的多边形内部的一切[255,0,0]（红色），如果你使用burn=[1,2,3],burn_values=[255,0,0]。我不确定点会发生什么 - 它们可能不会绘图。
Use gdal.RasterizeLayer(ds,bands,layer,burn_values, options = ["BURN_VALUE_FROM=Z"])if you want to use the Z values.
使用gdal.RasterizeLayer(ds,bands,layer,burn_values, options = ["BURN_VALUE_FROM=Z"])，如果你想使用的Z值。
I'm just pulling this from the tests you were looking at: http://svn.osgeo.org/gdal/trunk/autotest/alg/rasterize.py
我只是从您正在查看的测试中提取此内容：http: //svn.osgeo.org/gdal/trunk/autotest/alg/rasterize.py
Another approach - pull the polygon objects out, and draw them using shapely, which may not be attractive. Or look into geodjango (I think it uses openlayers to plot into browsers using JavaScript).
另一种方法 - 拉出多边形对象，并使用匀称绘制它们，这可能不吸引人。或者查看 geodjango（我认为它使用 openlayers 使用 JavaScript 绘制到浏览器中）。
Also, do you need to rasterize? A pdf export might be better, if you really want precision.
另外，你需要光栅化吗？如果您真的想要精度，pdf 导出可能会更好。
Actually, I think I found using Matplotlib (after extracting and projecting the features) was easier than rasterization, and I could get a lot more control.
实际上，我认为我发现使用 Matplotlib（在提取和投影特征之后）比光栅化更容易，而且我可以获得更多控制。
EDIT:
编辑：
A lower level approach is here:
较低级别的方法在这里：
http://svn.osgeo.org/gdal/trunk/gdal/swig/python/samples/gdal2grd.py\
http://svn.osgeo.org/gdal/trunk/gdal/swig/python/samples/gdal2grd.py\
Finally, you can iterate over the polygons (after transforming them into a local projection), and plot them directly. But you better not have complex polygons, or you will have a bit of grief. If you have complex polygons ... you are probably best off using shapely and r-tree from http://trac.gispython.org/labif you want to roll your own plotter.
最后，您可以迭代多边形（在将它们转换为局部投影之后），并直接绘制它们。但是你最好不要有复杂的多边形，否则你会有点悲伤。如果您有复杂的多边形……如果您想滚动自己的绘图仪，最好使用来自http://trac.gispython.org/lab 的shapely 和 r-tree 。
Geodjango might be a good place to ask .. they will know a lot more than me. Do they have a mailing list? There's also lots of python mapping experts around, but none of them seem to worry about this. I guess they just plot it in qGIS or GRASS or something.
Geodjango 可能是一个提问的好地方……他们会比我知道的多得多。他们有邮寄名单吗？周围也有很多 python 映射专家，但他们似乎都不担心这个。我猜他们只是在 qGIS 或 GRASS 或其他东西中绘制它。
Seriously, I hope that somebody who knows what they are doing can reply.
说真的，我希望知道自己在做什么的人可以回复。
