Rasterize a shapefile generate a black TIFF

Question

I run the following codes in jupyter notebook to rasterize the shapefile: building.shp, but it generates a black TIFF, which in QGIS shows only 0 and nan values.

Here are the source codes

Could anyone help me to find what the problem is?

import gdal
from osgeo import osr
from osgeo import ogr

raster_path = 'C:\\exercise_data\\training_manual_data\\epsg4326\\buildings.tif'
shapefile = 'C:\\exercise_data\\training_manual_data\\epsg4326\\buildings.shp'

# 1) opening the shapefile
source_ds = ogr.Open(shapefile)
source_layer = source_ds.GetLayer()



# 2) Creating the destination raster data source

pixelWidth = pixelHeight = 0.001 # depending how fine you want your raster
x_min, x_max, y_min, y_max = source_layer.GetExtent()
cols = int((x_max - x_min) / pixelHeight)
rows = int((y_max - y_min) / pixelWidth)

target_ds = gdal.GetDriverByName('GTiff').Create(raster_path, cols, rows, 1, gdal.GDT_Byte) 
target_ds.SetGeoTransform((x_min, pixelWidth, 0, y_min, 0, pixelHeight))
band = target_ds.GetRasterBand(1)
band.FlushCache()

# 4) Instead of setting a general burn_value, use optionsand set it to the attribute that contains the relevant unique value ["ATTRIBUTE=ID"]
gdal.RasterizeLayer(target_ds, [1], source_layer, burn_values=[1])

# 5) Adding a spatial reference
target_dsSRS = osr.SpatialReference()
target_dsSRS.ImportFromEPSG(4326)
target_ds.SetProjection(target_dsSRS.ExportToWkt())

Answer 1

I've had a similar problem when using the GDAL RasterizeLayer. I am not 100% sure what is happening behind the scenes when you are calling RasterizeLayer, and I am sure you can find out if you dig around in the source code. The problem is, I think, that RasterizeLayer is fussy on what fields(columns) in your attribute table it accepts as "valid" to create a raster from (e.g. it does not recognize a "string" column in your attribute table as a valid field for creating rasters). It seems to me (after downloading your shapefile) that there is no consistent\complete ID field in your dataset. What I would suggest is that you create a new field in your shapefile, something like "UID" and populate it with unique integer values (So that every shape has an ID). Do this with whatever GIS software you are using (QGIS or ARCMap).

Then replace this line:

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

with

gdal.RasterizeLayer(target_ds, [1], source_layer, options=["ATTRIBUTE=UID"])

Where the field "UID" is the one that you have just created.

---

I took another look and examined your code in a bit more detail. I managed to get it working. I made comments in the code and I explain it below. This is how I changed it:

import gdal
from osgeo import osr
from osgeo import ogr

raster_path = 'path on my pc'
shapefile = 'path on my pc'


# 1) opening the shapefile    

source_ds = ogr.Open(shapefile)
source_layer = source_ds.GetLayer()


# 2) Creating the destination raster data source

pixelWidth = pixelHeight = 0.0001 # depending how fine you want your raster ##COMMENT 1
x_min, x_max, y_min, y_max = source_layer.GetExtent()
cols = int((x_max - x_min) / pixelHeight)
rows = int((y_max - y_min) / pixelWidth)

target_ds = gdal.GetDriverByName('GTiff').Create(raster_path, cols, rows, 1, 
gdal.GDT_Float32) ##COMMENT 2

target_ds.SetGeoTransform((x_min, pixelWidth, 0, y_max, 0, -pixelHeight) ##COMMENT 3

# 5) Adding a spatial reference ##COMMENT 4
target_dsSRS = osr.SpatialReference()
target_dsSRS.ImportFromEPSG(4326)
target_ds.SetProjection(target_dsSRS.ExportToWkt())

band = target_ds.GetRasterBand(1) 
band.SetNoDataValue(-9999) ##COMMENT 5

gdal.RasterizeLayer(target_ds, [1], source_layer, options=["ATTRIBUTE = UID"]) ##COMMENT 6

target_ds = None ##COMMENT 7

COMMENT 1: I changed the pixel size and even with this size all the buildings were not rasterized, so you might want to make it even smaller

COMMENT 2: I changed the data type from 8 bit to 32 bit (GDT_Float32). This should not really make a difference, so try it out with the 8 bit (gdal.GDT_Byte).

COMMENT 3: I changed the y_min to y_max and I made the pixelHeight negative. This is because the ogr function GetExtent() and the gdal function SetGeoTransform or GetGeoTransform use different termanology for the x y limits. This is a bit confusing but ill try to explain. The variables for the Set/Get GeoTransform is x_upperlimit, x_size, x_skew, y_upperlimit, y_skew, y_size. The upperlimit variables refers to the top left point of your raster grid (I had to plot this out for myself to get the hang of it). The upper left point translates then to x_min and y_max on any Cartesian system, which is part of the variables of the ogr GetExtent function.

COMMENT 4: I just moved the spatial reference up a bit, it makes more sense to me here.

COMMENT 5: I just added a NoData value of -9999

COMMENT 6: As I described in my initial answer

COMMENT 7: I think this oneliner is the thing that solved the probelm. You need to remember to close the active data-set once you are done with it.

I hope this makes sense and that one of these changes will help you (especially comment 7).