将GDAL图层转换成栅格图像

31

编辑

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

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 是否适用于渲染地理数据(我不希望有抗锯齿效果,需要一个非常强大的绘图器,能够正确地绘制非常大和非常小的要素,可能来自“脏数据”(退化多边形等),有时指定在大坐标系下)?

这里,多边形是通过属性值进行区分的(颜色并不重要,我只想为每个属性值设置不同的颜色)。


1
谢谢Luper,今天对我非常有帮助!gdal的文档很难找到正确的信息... - yosukesabai
图片链接已经失效了,你能修复一下吗? - Yuchen
不好意思,我没有那些图片了,也在互联网档案馆里没找到它们 :( - Luper Rouch
嗨@Luper,有没有一种方法可以在不使用CopyDataSource进行复制的情况下完成这个操作?由于某种原因,这一步在我的程序中出现了段错误(gdal 2.0.1,python 2.7.6)。谢谢! - TM5
@TM5 你可以尝试,但是如果我没记错的话,如果你修改了原始数据源,它也会在磁盘上进行修改。 - Luper Rouch
显示剩余3条评论
1个回答

11

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

这是我能想到的最简单的方法。我记得以前手动编写过一些东西,但很丑陋。即使你需要安装一个单独的Windows版本(使python与之配合工作),然后设置一个XML-RPC服务器来在一个单独的python进程中运行它,qGIS也会更容易。

如果您可以正确地将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”可以使用属性吗? - Luper Rouch
1
另外,我不明白为什么我的测试结果是一张黑色的图片。 - Luper Rouch
2
我能够借助#gdal群里的人们成功使用RasterizeLayer。问题出在转换/范围偏移上,必须确保源和目标范围匹配,否则会被裁剪掉全部内容。实际上,文档中已经对此进行了解释,我不知道我是怎么错过它的 :D无论如何,还是感谢你的研究,我会接受你的答案并编辑我的问题,以修复代码。 - Luper Rouch

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接