GDAL层的栅格化

Question

编辑

下面是正确的方法，以及文档

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()的信息( docstring非常简洁，而且我在C或C++ API文档中找不到它。我只为java绑定找到了一个文档)。

我改编了一个单元测试，并在一个由多边形组成的.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()可以用于根据属性值不同颜色的图层进行栅格化吗？

如果不能，我该用什么？阿格是否适合于绘制地理数据(我希望没有反混叠和一个非常健壮的渲染器，能够正确地绘制非常大和非常小的特性，可能是从“脏数据”(简并多边形等)，有时在大坐标中指定)？

在这里，多边形是通过属性值来区分的(颜色并不重要，我只是想对属性值的每个值有一个不同的颜色)。

Answer 1

编辑:我想我会使用qGIS python绑定：绑定

这是我能想到的最简单的方法。我记得以前手滚动过什么东西，但很难看。qGIS会更容易一些，即使您必须单独安装一个(让python使用它)，然后设置一个XML服务器在单独的python进程中运行它。

我，你可以让GDAL正确地进行栅格化，这也很棒。

我已经有一段时间没有使用gdal了，但我猜：

如果不使用Z-值，burn_values是用于假颜色的。如果使用[255,0,0]，多边形内的所有内容都是burn=[1,2,3],burn_values=[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

另一种方法-拉出多边形对象，并使用形状绘制，这可能是没有吸引力的。或者查看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或草之类的地方画出来的。

说真的，我希望知道自己在做什么的人能够回答。