Растеризация полигона WKT (из геопанд) с помощью Python и OGR / GDAL

Как энтузиаст программирования, я в настоящее время создаю свой собственный «геопространственный» инструмент, как мне бы хотелось, чтобы он работал. Однако с самого начала я уже сталкиваюсь с проблемой. Мой инструмент должен работать с использованием GeoPandas для извлечения информации, а затем OGR / GDAL для редактирования данных, поскольку я хочу, чтобы он работал быстро. Мне нравится много анализировать и большие данные!

Отрезанный код проблемы должен растрировать один полигон GeoPandas. Я пытаюсь сделать это по этому пути. - извлечь с помощью геопанд WKT-полигон из полигона - Создать объект OGR с помощью WKT-полигона - растеризовать его с помощью GDAL.

Проблема, с которой я столкнулся, заключается в том, что я получаю только растр, который состоит из 0, а не из 0 и 1 ...

код показан ниже:

import geopandas as gpd
import ogr, osr
import gdal
import uuid

tf = r'f:test2.shp'

def vector_to_raster(source_layer, output_path, x_size, y_size, options, data_type=gdal.GDT_Byte):
    '''
    This method should create a raster object by burning the values of a source layer to values.
    '''

    x_min, x_max, y_min, y_max = source_layer.GetExtent()
    print(source_layer.GetExtent())
    x_resolution = int((x_max - x_min) / x_size)
    y_resolution = int((y_max - y_min) / -y_size)  
    print(x_resolution, y_resolution)

    target_ds = gdal.GetDriverByName(str('GTiff')).Create(output_path, x_resolution, y_resolution, 1, data_type)
    spatial_reference = source_layer.GetSpatialRef()         
    target_ds.SetProjection(spatial_reference.ExportToWkt())
    target_ds.SetGeoTransform((x_min, x_size, 0, y_max, 0, -y_size))
    gdal.RasterizeLayer(target_ds, [1], source_layer, options=options)
    target_ds.FlushCache()
    return target_ds


#create geopandas dataframe
gdf = gpd.read_file(tf)

#grab projection from the gdf
projection = gdf.crs['init']

#get geometry from 1 polygon (now just the 1st one)
polygon = gdf.loc[0].geometry 

#grab epsg from projection
epsg = int(projection.split(':')[1])

#create geometry
geom = ogr.CreateGeometryFromWkt(polygon.wkt)

#create spatial reference
proj = osr.SpatialReference()
proj.ImportFromEPSG(epsg) 

#get driver
rast_ogr_ds = ogr.GetDriverByName('Memory').CreateDataSource('wrk')

#create polylayer with projection
rast_mem_lyr = rast_ogr_ds.CreateLayer('poly', srs=proj)

#create feature
feat = ogr.Feature(rast_mem_lyr.GetLayerDefn())

#set geometry in feature
feat.SetGeometryDirectly(geom) 

#add feature to memory layer
rast_mem_lyr.CreateFeature(feat)

#create memory location
tif_output = '/vsimem/' + uuid.uuid4().hex + '.vrt'

#rasterize
lel = vector_to_raster(rast_mem_lyr, tif_output, 0.001, -0.001,['ATTRIBUTE=Shape__Len', 'COMPRESS=LZW', 'TILED=YES', 'NBITS=4'])

# output should consist of 0's and 1's
print(np.unique(lel.ReadAsArray()))

Большое спасибо человеку, который может подсказать мне верное направление :-).

Ваше здоровье!


python polygon gdal ogr rasterize
person Mr unkown    schedule 12.03.2019    source источник

Ответы (2)


arrow_upward
1
arrow_downward

Привет, я не запускал ваш код, но могу дать пару предложений. В настоящий момент вы растрируете в соответствии с полем 'Shape__Len' вашего многоугольника и указываете выходной растр как GDT_Byte (который может содержать только значения от 0 до 255), убедитесь, что 'Shape__Len' соответствует типу данных, или создайте новое поле в многоугольнике, содержащее целые числа от 0 до 255 для растеризации, или измените тип выходных данных на GDT_Float32.

В качестве альтернативы, если вам нужны только 1 и 0, вы можете записать значение 1 там, где есть многоугольник:

gdal.RasterizeLayer(target_ds, [1], source_layer,burn_values=[1])

Также разумно, поскольку вы создаете некоторые инструменты для себя, отслеживать / управлять своими значениями NoData. Если вы хотите отображать только растеризованные полигоны, вы можете добавить следующие шаги:

target_ds = gdal.GetDriverByName(str('GTiff')).Create(output_path, x_resolution, y_resolution, 1, data_type)
spatial_reference = source_layer.GetSpatialRef()         
target_ds.SetProjection(spatial_reference.ExportToWkt())
target_ds.SetGeoTransform((x_min, x_size, 0, y_max, 0, -y_size))
NowBand = target_ds.GetRasterBand(1) # ADD
NowBand.SetNoDataValue(0) # ADD NoData Value of your choice (according to your specified data type)
NowBand.Fill(0) # ADD Fill the band with NoData as to only display your rasterized features
gdal.RasterizeLayer(target_ds, [1], source_layer, burn_values=[1]) # if you only want to burn 1 in your values
target_ds.FlushCache()
person Jascha Muller    schedule 13.03.2019

arrow_upward
0
arrow_downward

Я тоже не запускал код, но на всякий случай, если кто-то ищет ответ, я думаю, что вы допустили ошибку в своей функции "vector_to_raster":

def vector_to_raster(source_layer, output_path, x_size, y_size, options, data_type=gdal.GDT_Byte):
'''
This method should create a raster object by burning the values of a source layer to values.
'''

x_min, x_max, y_min, y_max = source_layer.GetExtent()
print(source_layer.GetExtent())
x_resolution = (x_max - x_min) / x_size
y_resolution = (y_max - y_min) / -y_size
print(x_resolution, y_resolution)

target_ds = gdal.GetDriverByName(str('GTiff')).Create(output_path, x_resolution, y_resolution, 1, data_type)
spatial_reference = source_layer.GetSpatialRef()         
target_ds.SetProjection(spatial_reference.ExportToWkt())
target_ds.SetGeoTransform((x_min, x_resolution, 0, y_max, 0, y_resolution))
gdal.RasterizeLayer(target_ds, [1], source_layer, options=options)
target_ds.FlushCache()
return target_ds

x_resolution и y_resolution не следует преобразовывать как int. GeoTransform должен принимать x_resolution / y_resolution вместо x_size / (-) y_size.

person Selim    schedule 24.01.2020