J'essaie de graver un fichier de formes dans un raster en utilisant RasterizeLayer de GDAL. Je pré-crée une raster de zone d'intérêt à partir d'un fichier de formes différent, étant donné une taille de pixel spécifique. Cet AOI sert ensuite de base à toutes les rastérisations suivantes (même nombre de colonnes et de lignes, même projection et géotransformation).
Le problème se produit, cependant, lorsque je vais graver les formes dans leur propre raster, en fonction de la même taille de pixel et des mêmes projections. Le lien ci-dessous (pas assez de représentants pour publier l'image) montre le fichier de formes original en beige et le rose foncé où RasterizeLayer a gravé des données. Le rose clair correspond aux valeurs de nodata pour les données raster rose foncé. Le gris est l'AOI sur la base duquel la gravure du fichier de formes a été terminée.
Compte tenu de l'étendue des polygones de fichiers de formes, je m'attendrais à voir des valeurs de gravure dans les deux coins inférieurs, ainsi que les deux pixels sous les données qui s'affichent. Mais ce n'est évidemment pas le cas.
Voici le code que j'ai utilisé pour les générer. Toutes les formes ont été créées à l'aide de QGIS et ont toutes été créées dans la même projection. (Il convient de noter que le quadrillage dans l'image montrée était juste pour donner une idée de la taille de pixel que j'utilisais.)
from osgeo import ogr
from osgeo import gdal
aoi_uri = 'AOI_Raster.tif'
aoi_raster = gdal.Open(aoi_uri)
def new_raster_from_base(base, outputURI, format, nodata, datatype):
cols = base.RasterXSize
rows = base.RasterYSize
projection = base.GetProjection()
geotransform = base.GetGeoTransform()
bands = base.RasterCount
driver = gdal.GetDriverByName(format)
new_raster = driver.Create(str(outputURI), cols, rows, bands, datatype)
new_raster.SetProjection(projection)
new_raster.SetGeoTransform(geotransform)
for i in range(bands):
new_raster.GetRasterBand(i + 1).SetNoDataValue(nodata)
new_raster.GetRasterBand(i + 1).Fill(nodata)
return new_raster
shape_uri = 'activity_3.shp'
shape_datasource = ogr.Open(shape_uri)
shape_layer = shape_datasource.GetLayer()
raster_out = 'new_raster.tif'
raster_dataset = new_raster_from_base(aoi_raster, raster_out, 'GTiff',
-1, gdal.GDT_Int32)
band = raster_dataset.GetRasterBand(1)
nodata = band.GetNoDataValue()
band.Fill(nodata)
gdal.RasterizeLayer(raster_dataset, [1], shape_layer, burn_values=[1])
Est-ce un bogue dans GDAL, ou est-ce que RasterizeLayer brûle des données en se basant sur autre chose que la simple présence ou absence d'un polygone dans une zone de pixels spécifiée?
Les fichiers que j'utilisais se trouvent ici .
la source
Réponses:
J'ai joué avec GDALRasterizeLayers cette semaine et j'ai une assez bonne idée de ce qu'il fait. Par défaut, il pixellise un pixel si le centre du pixel se trouve dans le polygone. S'il n'y a rien au centre, il ne sera pas tramé, même s'il y a des parties d'un polygone dans les limites des pixels. Pour permettre à la pixellisation de fonctionner comme vous le souhaitez, essayez l'option "ALL_TOUCHED":
la source
['ALL_TOUCHED=TRUE']
, mais malheureusement, que seules les couches de polygones fixes. Mes couches de fichiers de formes ponctuelles sont toujours super bancales et affichent un pixel de l'endroit où elles sont placées.