Rasterizar una capa de GDAL

Question

Editar

Aquí es la forma correcta de hacerlo, y la documentation:

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) 

pregunta original

Busco información sobre cómo para usar osgeo.gdal.RasterizeLayer() (la docstring es muy sucinta, y no puedo encontrarla en los documentos de la API C o C++. Solo encontré un documento para el java bindings).

adapté un unit test y probado en un .shp hecha de polígonos:

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() 

Funciona muy bien, pero todo lo que obtengo es un .tif negro.

¿Para qué sirve el parámetro burn_values? ¿Se puede usar RasterizeLayer() para rasterizar una capa con características coloreadas de forma diferente según el valor de un atributo?

Si no puede, ¿qué debo usar? Es AGG adecuado para renderizar datos geográficos (quiero no antialiasing y un renderizador muy robusto, capaz de dibujar características muy grandes y muy pequeñas correctamente, posiblemente a partir de "datos sucios" (polígonos degenerados, etc.), y algunas veces especificado en coordenadas grandes)?

Por ejemplo, yo quiero ir de esta:
http://i.imagehost.org/0232/from.png http://i.imagehost.org/0232/from.png

A esto:
http://f.imagehost.org/0012/to_4.png http://f.imagehost.org/0012/to_4.png

Aquí, los polígonos se diferencian por el valor de un atributo (los colores no importan, Solo quiero tener una diferente para cada valor del atributo).

Answer 1

EDIT: supongo que haría uso de qgis de Python: http://www.qgis.org/wiki/Python_Bindings

Esa es la manera más fácil de lo que puedo imaginar. Recuerdo haber hecho rodar algo antes, pero es feo. qGIS sería más fácil, incluso si tuviera que hacer una instalación de Windows por separado (para que Python trabaje con ella) luego configure un servidor XML-RPC para ejecutarlo en un proceso de python separado.

I puede obtener GDAL para rasterizar correctamente eso también es genial.

No he utilizado gdal por un tiempo, pero aquí está mi suposición:

burn_values es para el color falso si no se utiliza Z-valores. Todo dentro de su polígono es [255,0,0] (rojo) si usa burn=[1,2,3],burn_values=[255,0,0]. No estoy seguro de qué pasa con los puntos; es posible que no tracen.

Utilice gdal.RasterizeLayer(ds,bands,layer,burn_values, options = ["BURN_VALUE_FROM=Z"]) si desea utilizar los valores Z.

sólo estoy tirando esto desde las pruebas que estaba mirando: http://svn.osgeo.org/gdal/trunk/autotest/alg/rasterize.py

Otro enfoque - tire del polígono objetos out, y dibujar usando bien formada, que puede no ser atractivo.O mire en geodjango (creo que usa openlayers para trazar en navegadores usando JavaScript).

Además, ¿necesita rasterizar? Una exportación en pdf podría ser mejor, si realmente quiere precisión.

En realidad, creo que el uso de Matplotlib (después de extraer y proyectar las características) era más fácil que la rasterización, y podía obtener mucho más control.

EDIT:

Un enfoque de nivel inferior está aquí:

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

Por último, se puede iterar sobre los polígonos (después de su transformación en una proyección local), y trazar con ellos directamente. Pero es mejor que no tengas polígonos complejos, o tendrás un poco de dolor. Si tiene polígonos complejos ... probablemente sea mejor que use shapely y r-tree desde http://trac.gispython.org/lab si desea hacer rodar su propio trazador.

Geodjango podría ser un buen lugar para preguntar ... sabrán mucho más que yo. ¿Tienen una lista de correo? También hay muchos expertos en mapeo de python, pero ninguno parece preocuparse por esto. Supongo que lo trazaron en qGIS o GRASS o algo así.

En serio, espero que alguien que sepa lo que están haciendo pueda responder.