栅格化GDAL层

编辑

这是正确的方法,以及文档:

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绑定的文档).

我调整了一个单元测试并在由多边形组成的.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是否适合渲染地理数据(我不需要抗锯齿和非常强大的渲染器,能够正确地绘制非常大和非常小的特征,可能来自"脏数据"(退化多边形等等),有时指定大坐标)？

这里,多边形通过属性的值来区分(颜色无关紧要,我只想为属性的每个值使用不同的颜色).

编辑:我想我会使用qGIS python绑定:http://www.qgis.org/wiki/Python_Bindings

这是我能想到的最简单的方法.我记得以前手上滚动的东西,但它很难看.即使您必须单独安装Windows(让python使用它),然后设置XML-RPC服务器以在单独的python进程中运行它,qGIS也会更容易.

我可以让GDAL正确地光栅化,这也很棒.

我有一段时间没用过gdal了,但这是我的猜测:

burn_values如果您不使用Z值,则为假色.[255,0,0]如果使用,多边形内的所有内容都是(红色)burn=[1,2,3],burn_values=[255,0,0].我不确定点会发生什么 - 他们可能不会策划.

使用gdal.RasterizeLayer(ds,bands,layer,burn_values, options = ["BURN_VALUE_FROM=Z"]),如果你想使用的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中的匀称和r-tree .

Geodjango可能是个好地方......他们会比我知道更多.他们有邮件列表吗？还有很多python地图专家,但他们似乎都不担心这一点.我猜他们只是在qGIS或GRASS或其他东西中绘制它.

说真的,我希望有人知道他们在做什么可以回复.