Convert an OGR File to a Raster

Question

i tried out this recipe from the GDAL/OGR cookbook:

from osgeo import gdal, ogr
# Define pixel_size and NoData value of new raster
pixel_size = 25
NoData_value = -9999
# Filename of input OGR file
vector_fn = 'test.shp'
# Filename of the raster Tiff that will be created
raster_fn = 'test.tif'
# Open the data source and read in the extent
source_ds = ogr.Open(vector_fn)
source_layer = source_ds.GetLayer()
x_min, x_max, y_min, y_max = source_layer.GetExtent()
# 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(raster_fn, x_res, y_res,        1, gdal.GDT_Byte)
target_ds.SetGeoTransform((x_min, pixel_size, 0, y_max, 0, -pixel_size))
band = target_ds.GetRasterBand(1)
band.SetNoDataValue(NoData_value)
# Rasterize
gdal.RasterizeLayer(target_ds, [1], source_layer, burn_values=[0])

I recieve a rectangular raster in the area were my polygon is located. I want to mask the raster exactly to the boundaries of a multipolygon shapefile. Is there an option to GetGeometry from the shapefile and use this as parameter instead of(x, y).

Have you tried using the ALL_TOUCHED=TRUE configuration suggested here:https://gis.stackexchange.com/a/31658/83630 — Logan Byers, Jan 21 '19 at 21:44
Perhaps you will need to assert your NoData value on the target_ds, at the moment it's probably all 0's (GDAL makes no guarantee of any particular value being set in the raster https://www.gdal.org/classGDALRasterBand.html#a6c7625813411cb5f0500e72d963249c2) and your burn value is also 0. Use band.Fill(NoData_value) https://www.gdal.org/classGDALRasterBand.html#a55bf20527df638dc48bf25e2ff26f353 to fill the raster with your NoData first then rasterize your layer with value 0. — Michael Stimson, Jan 21 '19 at 23:21

score 1 · Answer 1 · answered Mar 11 '19 at 15:14

Like @Michael Stimson indicated in the comments, create a NoData raster first and then Rasterize, see the code adaptation below:

from osgeo import gdal, ogr
import numpy as np ##NEW##
pixel_size = 25
NoData_value = -9999

vector_fn = 'test.shp'
raster_fn = 'test.tif'

source_ds = ogr.Open(vector_fn)
source_layer = source_ds.GetLayer()
x_min, x_max, y_min, y_max = source_layer.GetExtent()

x_res = int((x_max - x_min) / pixel_size)
y_res = int((y_max - y_min) / pixel_size)
target_ds = gdal.GetDriverByName('GTiff').Create(raster_fn, x_res, y_res, 1, gdal.GDT_Byte)
target_ds.SetGeoTransform((x_min, pixel_size, 0, y_max, 0, -pixel_size))
band = target_ds.GetRasterBand(1)
band.SetNoDataValue(NoData_value)

NullArr = np.full((y_res, x_res), NoData_value) ##NEW##
band.WriteArray(NullArr) ##NEW##

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

target_ds = None

Hope it helps

Convert an OGR File to a Raster

1 Answers1