pixellisation une couche GDAL

Question

Modifier

Voici la bonne façon de le faire, et le 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) 

question originale

Je cherche des informations sur la façon dont pour utiliser osgeo.gdal.RasterizeLayer() (la docstring est très succinte, et je ne la trouve pas dans les docs d'API C ou C++. Je n'ai trouvé qu'un doc pour le java bindings).

Je me suis adapté et un unit test essayé sur un shp en polygones:

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

Il fonctionne très bien, mais tout ce que j'obtiens est un .tif noir.

À quoi sert le paramètre burn_values? Peut-on utiliser RasterizeLayer() pour pixelliser une couche avec des entités colorées différemment en fonction de la valeur d'un attribut?

Si ce n'est pas possible, que dois-je utiliser? Est AGG adapté pour le rendu de données géographiques (je veux non antialiasing et un moteur de rendu très robuste, capable de dessiner correctement de très grandes et très petites fonctionnalités, éventuellement à partir de "données corrompues" (polygones dégénérés, etc ...), et parfois spécifié en grandes coordonnées)?

Par exemple, je veux aller de ceci:
http://i.imagehost.org/0232/from.png http://i.imagehost.org/0232/from.png

à ceci:
http://f.imagehost.org/0012/to_4.png http://f.imagehost.org/0012/to_4.png

Ici, les polygones sont différenciés par la valeur d'un attribut (les couleurs ne comptent pas, Je veux juste en avoir un différent pour chaque valeur de l'attribut).

Answer 1

EDIT: Je pense que j'utiliser les liaisons python QGIS: http://www.qgis.org/wiki/Python_Bindings

Ce que je peux penser est la meilleure façon de. Je me souviens avoir roulé quelque chose avant, mais c'est moche. qGIS serait plus facile, même si vous deviez faire une installation Windows séparée (pour que python fonctionne avec), puis configurer un serveur XML-RPC pour l'exécuter dans un processus python séparé.

Vous pouvez obtenir GDAL pour pixelliser correctement, c'est génial aussi.

Je ne l'ai pas utilisé gdal pendant un certain temps, mais voici ma conjecture:

burn_values est pour fausse couleur si vous ne l'utilisez pas les valeurs Z. Tout à l'intérieur de votre polygone est [255,0,0] (rouge) si vous utilisez burn=[1,2,3],burn_values=[255,0,0]. Je ne suis pas sûr de ce qui arrive aux points - ils ne peuvent pas comploter. Pour utiliser les valeurs Z, utilisez gdal.RasterizeLayer(ds,bands,layer,burn_values, options = ["BURN_VALUE_FROM=Z"]).

Je suis juste tirer cela des tests que vous regardiez: http://svn.osgeo.org/gdal/trunk/autotest/alg/rasterize.py

Une autre approche - tirer le polygone des objets sur, et d'en tirer les utiliser bien faite, qui peut ne pas être attrayant.Ou regardez dans geodjango (je pense qu'il utilise Openlayers pour tracer dans les navigateurs en utilisant JavaScript).

De même, avez-vous besoin de pixelliser? Une exportation pdf pourrait être meilleure, si vous voulez vraiment la précision.

En fait, je pense avoir trouvé que l'utilisation de Matplotlib (après avoir extrait et projeté les fonctions) était plus facile que la pixellisation, et que je pouvais obtenir beaucoup plus de contrôle.

EDIT:

Une approche de niveau inférieur est ici:

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

Enfin, vous pouvez parcourir les polygones (après les transformer en une projection locale), et les tracer directement. Mais vous feriez mieux de ne pas avoir de polygones complexes, sinon vous aurez un peu de chagrin. Si vous avez des polygones complexes ... il vaut probablement mieux utiliser http://trac.gispython.org/lab si vous voulez lancer votre propre traceur.

Geodjango pourrait être un bon endroit à demander .. ils sauront beaucoup plus que moi. Ont-ils une liste de diffusion? Il y a aussi beaucoup d'experts de la cartographie python, mais aucun d'entre eux ne semble s'inquiéter à ce sujet. Je suppose qu'ils ont juste tracé dans qGIS ou GRASS ou quelque chose. Sérieusement, j'espère que quelqu'un qui sait ce qu'ils font peut répondre.