30

編集

これを行う適切な方法とドキュメントは次のとおりです。

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)

元の質問

使用方法に関する情報を探しています(docstring は非常に簡潔であり、C または C++ API ドキュメントで見つけることができません。Java バインディングosgeo.gdal.RasterizeLayer()のドキュメントしか見つかりませんでした)。

単体テストを適応させ、ポリゴンで作成された .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は地理データのレンダリングに適していますか (アンチエイリアシングがなく、非常に堅牢なレンダラーが必要で、おそらく「ダーティ データ」(縮退ポリゴンなど) から、非常に大きなフィーチャと非常に小さなフィーチャを正しく描画できる必要があります)。座標)?

ここでは、ポリゴンは属性の値によって区別されます (色は関係ありません。属性の値ごとに異なる色を使用したいだけです)。

4

1 に答える 1

10

編集:私はqGIS pythonバインディングを使用すると思います: http://www.qgis.org/wiki/Python_Bindings

それが私が考えることができる最も簡単な方法です。前に何かを手で転がしたのを覚えていますが、それは醜いです。qGISは、別のWindowsインストールを作成して(pythonを動作させるために)、別のpythonプロセスで実行するようにXML-RPCサーバーをセットアップする必要がある場合でも、より簡単です。

GDAL を適切にラスタライズすることができれば、それも素晴らしいことです。

しばらく gdal を使用していませんが、私の推測は次のとおりです。

burn_valuesZ 値を使用しない場合の偽色用です。を使用すると、ポリゴン内のすべてが[255,0,0](赤) になりますburn=[1,2,3],burn_values=[255,0,0]。ポイントがどうなるかはわかりません-プロットされない可能性があります。

Z 値を使用するgdal.RasterizeLayer(ds,bands,layer,burn_values, options = ["BURN_VALUE_FROM=Z"])場合に使用します。

私はあなたが見ていたテストからこれを引っ張っています: http://svn.osgeo.org/gdal/trunk/autotest/alg/rasterize.py

別のアプローチ - ポリゴン オブジェクトを引き出し、見栄えがよくない場合があります。または、geodjango を調べます (openlayers を使用して、JavaScript を使用してブラウザーにプロットすると思います)。

また、ラスタライズする必要がありますか?本当に精度が必要な場合は、pdf エクスポートの方が良いかもしれません。

実際、Matplotlib を使用すると (フィーチャを抽出して投影した後)、ラスタライズよりも簡単で、より多くの制御を行うことができることがわかったと思います。

編集:

低レベルのアプローチは次のとおりです。

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

最後に、(ローカル投影に変換した後) ポリゴンを反復処理し、直接プロットできます。ただし、複雑なポリゴンは使用しないほうがよいでしょう。複雑な多角形がある場合...独自のプロッターを作成したい場合は、 http: //trac.gispython.org/lab の形状の良い r-tree を使用するのがおそらく最善です。

Geodjango に尋ねるのが良いかもしれません..彼らは私よりも多くのことを知っているでしょう. 彼らはメーリングリストを持っていますか? Python マッピングの専門家もたくさんいますが、誰もこれについて心配していないようです。彼らはそれをqGISやGRASSなどでプロットしているだけだと思います。

真剣に、彼らが何をしているのか知っている誰かが答えてくれることを願っています。

于 2010-02-08T10:45:45.797 に答える