python - pyplot - Rasterizando una capa GDAL

title plt python

Question

Editar

Aquí está 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

Estoy buscando información sobre cómo usar osgeo.gdal.RasterizeLayer() (la cadena de documentos es muy sucinta, y no puedo encontrarla en los documentos de la API C o C ++. Solo encontré un documento para los enlaces de Java ).

Adapté una prueba de unidad y la probé en un .shp hecho 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 bien, pero todo lo que obtengo es un .tif negro.

¿Para qué es 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 representar datos geográficos (no quiero antialiasing y un renderizador muy robusto, capaz de dibujar características muy grandes y muy pequeñas correctamente, posiblemente de "datos sucios" (polígonos degenerados, etc.), y algunas veces se especifican en grandes coordenadas)?

Por ejemplo quiero pasar de esto:
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 uno diferente para cada valor del atributo).

Answer 1

EDITAR: Supongo que usaría los enlaces qGIS python: http://www.qgis.org/wiki/Python_Bindings

Esa es la forma más fácil en la que puedo pensar. Recuerdo a mano rodar algo antes, pero es feo. qGIS sería más fácil, incluso si tuviera que realizar una instalación de Windows por separado (para que Python funcione con ella) y luego configurar un servidor XML-RPC para ejecutarlo en un proceso de Python por separado.

Si consigues que GDAL se rasterice correctamente, eso también es genial.

No he usado gdal por un tiempo, pero aquí está mi conjetura:

burn_values es para el color falso si no usas valores Z. Todo lo que está 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 lo que sucede con los puntos, podrían no trazar.

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

Solo estoy sacando esto de las pruebas que estaba viendo: http://svn.osgeo.org/gdal/trunk/autotest/alg/rasterize.py

Otro enfoque: sacar los objetos poligonales y dibujarlos usando una figura bien formada, que puede no ser atractiva. O busque en geodjango (creo que usa openlayers para trazar en navegadores usando JavaScript).

Además, ¿necesitas rasterizar? Una exportación de pdf podría ser mejor, si realmente quieres precisión.

En realidad, creo que el uso de Matplotlib (después de extraer y proyectar las funciones) fue más fácil que la rasterización, y podría tener mucho más control.

EDITAR:

Un enfoque de nivel inferior está aquí:

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

Finalmente, puede iterar sobre los polígonos (después de transformarlos en una proyección local) y trazarlos 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 de http://trac.gispython.org/lab si desea rodar su propio plotter.

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

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