1

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:gis.stackexchange.com/a/31658/83630 – Logan Byers Jan 21 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 gdal.org/…) and your burn value is also 0. Use band.Fill(NoData_value) gdal.org/… to fill the raster with your NoData first then rasterize your layer with value 0. – Michael Stimson Jan 21 at 23:21
0

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.