使用Python和OGR/GDAL栅格化WKT面(来自geopandas)
使用Python和OGR/GDAL栅格化WKT面(来自geopandas)
提问于 2019-03-12 23:49:23
作为一个编码爱好者，我目前正在创建我自己的“地理空间”工具，因为我希望它是如何工作的。然而，从一开始，我就已经面临一个问题。我的工具应该使用GeoPandas来提取信息，然后使用OGR/GDAL进行数据编辑，因为我希望它能快速工作。我喜欢分析很多大数据！
这个问题的代码片段应该栅格化一个GeoPandas多边形。我尝试使用以下路径来实现此目的。-使用geopandas提取WKT-polygon from a polygon -使用WKT-polygon创建一个OGR-feature -使用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()))非常感谢能给我正确方向提示的人:-)。
干杯!
回答 2
Stack Overflow用户
发布于 2019-03-13 15:30:56
嗨，我没有运行你的代码，但我可以给你一些建议。此时，您正在根据多边形的'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()Stack Overflow用户
发布于 2020-01-24 22:56:24
我也没有运行代码，但为了防止有人在寻找答案，我认为您在"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_dsx_resolution和y_resolution不应强制转换为int。分辨率应为x_ GeoTransform /y_y_size，而不是x_size/(-)分辨率。
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接：
https://stackoverflow.com/questions/55125643
复制相关文章
相似问题
腾讯云开发者
