Править
Вот надлежащий способ сделать это, и документация:
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)
Исходный вопрос
Я ищу информацию о том, как использовать osgeo.gdal.RasterizeLayer()
(docstring очень сжат, и я не могу найти его в C или C++ документами API. Я только нашел документ для привязки Java).
Я адаптировал модульный тест и попробовал его на .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 подходит для рендеринга географических данных (я не хочу сглаживания и очень устойчивого рендерера, который в состоянии потянуть очень большие и очень маленькие функции правильно, возможно от "грязных данных" (вырожденные полигоны, и т.д....), и иногда указываемый в больших координатах)?
Здесь, полигоны дифференцируются значением атрибута (цвета не имеют значения, я просто хочу иметь другой для каждого значения атрибута).
Из стандарта:
const char c = 'c';
char* pc;
const char** pcc = &pc; // not allowed
*pcc = &c;
*pc = 'C'; // would allow to modify a const object
-121--1275680- Вот еще одна реализация, которая рекурсивно распечатает все подклассы с отступом.
def findsubclass(baseclass, indent=0): if indent == 0: print "Subclasses of %s are:" % baseclass.__name__ indent = indent + 1 for c in baseclass.__subclasses__(): print "-"*indent*4 + ">" + c.__name__ findsubclass(c, indent)-121--2214929-
EDIT: Я думаю, я бы использовал привязки python qGIS: http://www.qgis.org/wiki/Python_Bindings
Это самый простой способ, который я могу придумать. Я помню, как что-то перекатывали до этого, но это некрасиво. qGIS было бы проще, даже если бы вам пришлось сделать отдельную установку Windows (чтобы заставить python работать с ней), а затем набор сервер XML-RPC, чтобы запустить его в отдельном процессе python.
Я могу заставить GDAL растрировать правильно, что тоже здорово.
Я давно не использовал gdal, но вот мое предположение:
burn _ значения
для ложного цвета, если вы не используете Z-значения. Все внутри полигона [255,0,0]
(красный), если используется burn = [1,2,3], burn _ values = [255,0,0]
. Я не уверен, что происходит с точками - они могут не заговариваться.
Если требуется использовать значения Z, используйте gdal.RasterureLayer (ds,bands,layer,burn_values, options = [«BURN _ VALUE _ FROM = Z»])
.
Я просто извлекаю это из тестов, на которые вы смотрели: http://svn.osgeo.org/gdal/trunk/autotest/alg/rasterize.py
Другой подход - вытащить многоугольные объекты, и нарисовать их с помощью формы, которая может не быть привлекательной. Или посмотрите на геоджанго (я думаю, что он использует открытые слои для печати в браузерах с помощью JavaScript).
Вам также нужно растрировать? Экспорт PDF может быть лучше, если вам действительно нужна точность.
На самом деле, я думаю, что использование Matplotlib (после извлечения и проецирования элементов) было проще, чем растеризация, и я мог получить гораздо больше контроля.
EDIT:
Здесь используется подход более низкого уровня:
http://svn.osgeo.org/gdal/trunk/gdal/swig/python/samples/gdal2grd.py \
Наконец, можно выполнить итерацию над многоугольниками (после преобразования их в локальную проекцию) и непосредственно построить их график. Но вам лучше не иметь сложных полигонов, или у вас будет немного горя. Если у вас сложные многоугольники... Возможно, лучше использовать формное и r-дерево из http://trac.gispython.org/lab , если вы хотите свернуть свой плоттер.
Геоджанго может быть хорошим местом для запроса.. они узнают гораздо больше меня. У них есть список рассылки? Вокруг также много экспертов по картированию питонов, но никто из них не беспокоится об этом. Я думаю, они просто планируют его в qGIS или GRASS или что-то еще.
Серьезно, я надеюсь, что кто-то, кто знает то, что они делают, может ответить.