对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()（文档说明得很简洁，我在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适合用来渲染地理数据吗（我希望不使用抗锯齿，并且需要一个非常强大的渲染器，能够正确绘制非常大和非常小的特征，可能还要处理一些“脏数据”（退化多边形等），有时坐标值也很大）？

在这里，这些多边形是通过一个属性的值来区分的（颜色无所谓，我只是想为每个属性值设置不同的颜色）。

gdal shapefile rendering rasterization geospatial data attribute-based styling agglomerative geometry coordinate systems

对GDAL图层进行栅格化

1 个回答

撰写回答