Как энтузиаст кодирования я в настоящее время создаю свой собственный «геопространственный» инструмент, как я бы хотел, чтобы он работал.Однако с самого начала я уже сталкиваюсь с проблемой.Мой инструмент должен работать, используя 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()))
Большое спасибо человеку, который может дать мне подсказку в правильном направлении: -).
Ура!