Rasterize một lớp GDAL

Question

Sửa

Đây là cách thích hợp để làm điều đó, và 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) 

gốc câu hỏi

Tôi đang tìm thông tin về cách để sử dụng osgeo.gdal.RasterizeLayer() (chuỗi tài liệu rất ngắn gọn và tôi không thể tìm thấy nó trong tài liệu API C hoặc C++. Tôi chỉ tìm thấy tài liệu cho số java bindings).

tôi chuyển thể một unit test và thử nó trên một .shp làm bằng đa giác:

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

Nó chạy tốt, nhưng tất cả tôi có được là một .tif đen.

Thông số burn_values là gì? Có thể sử dụng RasterizeLayer() để rasterize một lớp với các đối tượng có màu khác nhau dựa trên giá trị của thuộc tính không?

Nếu không, tôi nên sử dụng thông tin gì? Là AGG thích hợp để hiển thị dữ liệu địa lý (tôi muốn không khử răng cưa và trình kết xuất đồ họa mạnh mẽ, có thể vẽ chính xác các tính năng rất lớn và nhỏ, có thể từ "dữ liệu bẩn" (đa giác thoái hóa, v.v ...) và đôi khi được chỉ định ở tọa độ lớn)?

Ví dụ tôi muốn đi từ này:
http://i.imagehost.org/0232/from.png http://i.imagehost.org/0232/from.png

Để này:
http://f.imagehost.org/0012/to_4.png http://f.imagehost.org/0012/to_4.png

Ở đây, đa giác được phân biệt bởi giá trị của một thuộc tính (màu sắc không quan trọng, Tôi chỉ muốn có một cái khác nhau cho mỗi giá trị của thuộc tính).

Answer 1

EDIT: Tôi đoán tôi muốn sử dụng QGIS bindings python: http://www.qgis.org/wiki/Python_Bindings

Đó là cách đơn giản nhất tôi có thể nghĩ đến. Tôi nhớ tay đang lăn cái gì đó trước đây, nhưng nó xấu xí. qGIS sẽ dễ dàng hơn, ngay cả khi bạn phải cài đặt Windows riêng biệt (để có được python để làm việc với nó), sau đó thiết lập một máy chủ XML-RPC để chạy nó trong một quá trình python riêng biệt.

Tôi có thể giúp GDAL rasterize đúng cách tuyệt vời.

tôi đã không sử dụng GDAL một thời gian, nhưng đây là suy đoán của tôi:

burn_values là cho màu giả nếu bạn không sử dụng Z-giá trị. Mọi thứ bên trong đa giác của bạn là [255,0,0] (màu đỏ) nếu bạn sử dụng burn=[1,2,3],burn_values=[255,0,0]. Tôi không chắc điều gì xảy ra với các điểm - họ có thể không âm mưu.

Sử dụng gdal.RasterizeLayer(ds,bands,layer,burn_values, options = ["BURN_VALUE_FROM=Z"]) nếu bạn muốn sử dụng giá trị Z.

Tôi chỉ kéo này từ các bài kiểm tra bạn đang tìm kiếm tại địa chỉ: http://svn.osgeo.org/gdal/trunk/autotest/alg/rasterize.py

cách tiếp cận khác - kéo đa giác đối tượng ra, và vẽ chúng bằng cách sử quyến rũ, mà có thể không hấp dẫn.Hoặc nhìn vào geodjango (Tôi nghĩ rằng nó sử dụng openlayers để âm mưu vào các trình duyệt bằng cách sử dụng JavaScript).

Ngoài ra, bạn có cần phải quét lại không? Xuất pdf có thể tốt hơn, nếu bạn thực sự muốn độ chính xác.

Thực ra, tôi nghĩ rằng tôi đã tìm thấy bằng cách sử dụng Matplotlib (sau khi trích xuất và chiếu các đối tượng) dễ hơn rasterization và tôi có thể kiểm soát nhiều hơn.

EDIT:

Một cách tiếp cận mức độ thấp hơn là ở đây:

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

Cuối cùng, bạn có thể lặp qua các đa giác (sau khi chuyển đổi chúng thành một chiếu địa phương), và âm mưu chúng trực tiếp. Nhưng bạn tốt hơn không có đa giác phức tạp, hoặc bạn sẽ có một chút đau buồn. Nếu bạn có các đa giác phức tạp ... có lẽ bạn nên sử dụng shapely và r-tree từ http://trac.gispython.org/lab nếu bạn muốn cuộn bản vẽ của riêng bạn.

Geodjango có thể là một nơi tốt để hỏi .. họ sẽ biết nhiều hơn tôi. Họ có danh sách gửi thư không? Ngoài ra còn có rất nhiều chuyên gia lập bản đồ python xung quanh, nhưng không ai trong số họ có vẻ lo lắng về điều này. Tôi đoán họ chỉ vẽ nó trong qGIS hoặc GRASS hoặc một cái gì đó.

Nghiêm túc, tôi hy vọng rằng ai đó biết những gì họ đang làm có thể trả lời.