eine GDAL Schicht

Question

bearbeiten

Hier Rastern ist der richtige Weg, es zu tun, und die 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) 

Original Frage

Ich bin der Suche nach Informationen darüber, wie zu verwenden osgeo.gdal.RasterizeLayer() (der Docstring ist sehr kurz und ich kann es nicht in der C oder C++ API-Dokumente finden. Ich fand nur ein Dokument für die java bindings).

ich angepasst, um ein unit test und versuchte es auf einem SHP von Polygonen:

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

Es fein läuft, aber alles, was ich bekommen ist ein schwarzer .tif.

Was ist der burn_values Parameter für? Kann RasterizeLayer() verwendet werden, um eine Ebene mit Features zu rastern, die basierend auf dem Wert eines Attributs unterschiedlich gefärbt sind?

Wenn es nicht kann, was soll ich verwenden? Ist AGG geeignet zum Rendern geografischer Daten (Ich möchte keine Antialiasing und ein sehr robuster Renderer, in der Lage, sehr große und sehr kleine Features richtig zu zeichnen, möglicherweise von "schmutzigen Daten" (degenerierte Polygone, etc ...), und manchmal angegeben in großen Koordinaten)?

Zum Beispiel habe ich von diesem gehen wollen:
http://i.imagehost.org/0232/from.png http://i.imagehost.org/0232/from.png

Um dies:
http://f.imagehost.org/0012/to_4.png http://f.imagehost.org/0012/to_4.png

Hier werden die Polygone durch den Wert eines Attributs unterschieden werden (die Farben sind nicht wichtig, Ich möchte nur einen anderen für jeden Wert des Attributs haben).

Answer 1

EDIT: Ich glaube, ich qgis Python-Anbindung verwenden würde: http://www.qgis.org/wiki/Python_Bindings

, dass der einfachste Weg, ich denken kann. Ich erinnere mich daran, wie ich vorher etwas rollte, aber es ist hässlich. qGIS wäre einfacher, selbst wenn Sie eine separate Windows-Installation durchführen müssten (damit Python damit arbeiten kann), richten Sie dann einen XML-RPC-Server ein, um ihn in einem separaten Python-Prozess auszuführen.

Ich kann GDAL richtig zu rastern, das ist auch toll.

Ich habe nicht gdal für eine Weile benutzt, aber hier ist meine Vermutung:

burn_values für falsche Farbe ist, wenn Sie keine Z-Werte verwenden. Alles in Ihrem Polygon ist [255,0,0] (rot), wenn Sie burn=[1,2,3],burn_values=[255,0,0] verwenden. Ich bin mir nicht sicher, was mit Punkten passiert - sie könnten nicht plotten.

Verwenden Sie gdal.RasterizeLayer(ds,bands,layer,burn_values, options = ["BURN_VALUE_FROM=Z"]), wenn Sie die Z-Werte verwenden möchten.

Ich ziehe gerade dies aus den Tests an Sie suchen: http://svn.osgeo.org/gdal/trunk/autotest/alg/rasterize.py

Ein anderer Ansatz - ziehen das Polygon aus Objekte und ziehen sie mit wohlgeformt, die nicht attraktiv sein kann.Oder schauen Sie sich Geodjango an (ich denke, es verwendet Openlayer, um in Browsern JavaScript zu verwenden).

Müssen Sie auch rasterisieren? Ein PDF-Export ist möglicherweise besser, wenn Sie wirklich Präzision wünschen.

Eigentlich fand ich Matplotlib (nach dem Extrahieren und Projizieren der Features) war einfacher als Rasterung, und ich könnte viel mehr Kontrolle bekommen.

EDIT:

Ein unterer Ebene Ansatz ist hier:

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

Schließlich können Sie die Polygone iterieren (nach ihnen in einen lokalen Vorsprung verwandeln), und sie direkt zeichnen. Aber du solltest lieber keine komplexen Polygone haben oder du wirst ein bisschen Kummer haben. Wenn Sie komplexe Polygone haben ... sind Sie wahrscheinlich am besten mit formschönen und r-Baum von http://trac.gispython.org/lab, wenn Sie Ihren eigenen Plotter rollen wollen.

Geodjango könnte ein guter Ort sein zu fragen .. sie werden viel mehr als ich wissen. Haben sie eine Mailingliste? Es gibt auch viele Python Mapping-Experten, aber keiner von ihnen scheint sich darüber Sorgen zu machen. Ich schätze, sie zeichnen es einfach in qGIS oder GRASS oder so.

Ernsthaft hoffe ich, dass jemand, der weiß, was sie tun, antworten kann.