Растеризация слоя GDAL

Edit

Вот правильный способ сделать это, и документация:

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 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 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 очень кратка, и я не могу найти ее в документах API C или С++. Я нашел только документ для привязки java).

Я адаптировал unit test и попробовал его на .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 подходит для рендеринга географических данных (я хочу не сглаживание и очень надежный рендерер, способный рисовать очень большие и очень маленькие функции, возможно, из-за "грязных данных" (вырожденных полигонов и т.д.), а иногда и в больших координатах)?

Например, я хочу перейти от этого:
http://i.imagehost.org/0232/from.png http://i.imagehost.org/0232/from.png

К этому:
http://f.imagehost.org/0012/to_4.png http://f.imagehost.org/0012/to_4.png

Здесь полигоны дифференцируются по значению атрибута (цвета не имеют значения, я просто хочу иметь другое значение для каждого значения атрибута).

Ответы

Ответ 1

EDIT: Я предполагаю, что буду использовать привязки python qGIS: http://www.qgis.org/wiki/Python_Bindings

Это самый простой способ, о котором я могу думать. Я помню, как раньше что-то покачивала, но это уродливо. qGIS было бы проще, даже если вам пришлось сделать отдельную установку Windows (чтобы заставить python работать с ней), затем настройте сервер XML-RPC для запуска его в отдельном процессе python.

Я могу заставить 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 (после извлечения и проецирования функций) было проще, чем растеризация, и я мог получить намного больше контроля.

EDIT:

Ниже приведен подход ниже:

http://svn.osgeo.org/gdal/trunk/gdal/swig/python/samples/gdal2grd.py\

Наконец, вы можете выполнять итерацию по многоугольникам (после преобразования их в локальную проекцию) и непосредственно строить их. Но лучше не иметь сложных полигонов, иначе у вас будет немного горя. Если у вас сложные полигоны... вам, вероятно, лучше всего использовать фигурные и r-деревья из http://trac.gispython.org/lab, если вы хотите свернуть свой собственный плоттер.

Geodjango может быть хорошим местом, чтобы спросить.. они будут знать намного больше, чем я. У них есть список рассылки? Там также много экспертов по картированию питона, но никто из них, похоже, не беспокоится об этом. Думаю, они просто замышляют это в qGIS или GRASS или что-то в этом роде.

Серьезно, я надеюсь, что кто-то, кто знает, что они делают, может ответить.