Rasterizing polygons in shapefile using GDAL

Question

I have the following polygons shapefile consisting of polygons of the five boroughs of New York City from "NYC Open Data", linked here: https://data.cityofnewyork.us/City-Government/Borough-Boundaries/tqmj-j8zm

I want to take this shapefile and rasterize it to a .tif file using gdal in Python.

I am using the following python script example to accomplish this task, cited from "Python GDAL/OGR Cookbook 1.0 documentation", linked here: https://pcjericks.github.io/py-gdalogr-cookbook/raster_layers.html

This is what I am using:

from osgeo import gdal, ogr

gdal.UseExceptions()

# Define pixel_size and NoData value of new raster
pixel_size = 25
NoData_value = -9999

# Filename of input OGR file
vector_fn = 'C:/Users/MyName/Desktop/NYC_GIS/Boroughs_Test/geo_export_d92842a7-2631-4fe5-ab83-ee070743980c.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])

However, this returns the following error message:

RuntimeError                              Traceback (most recent call last)
<ipython-input-32-0e6d6c7404d4> in <module>
     21 x_res = int((x_max - x_min) / pixel_size)
     22 y_res = int((y_max - y_min) / pixel_size)
---> 23 target_ds = gdal.GetDriverByName('GTiff').Create(raster_fn, x_res, y_res, 1, gdal.GDT_Byte)
     24 target_ds.SetGeoTransform((x_min, pixel_size, 0, y_max, 0, -pixel_size))
     25 band = target_ds.GetRasterBand(1)

~\Anaconda3\lib\site-packages\osgeo\gdal.py in Create(self, *args, **kwargs)
   1892     def Create(self, *args, **kwargs):
   1893         """Create(Driver self, char const * utf8_path, int xsize, int ysize, int bands=1, GDALDataType eType, char ** options=None) -> Dataset"""
-> 1894         return _gdal.Driver_Create(self, *args, **kwargs)
   1895 
   1896 

RuntimeError: Attempt to create 0x0 dataset is illegal,sizes must be larger than zero.

How do I address this error so that the raster .tif file is correctly created and sent to my directory folder? I am not sure why this message is telling me I am creating a 0x0 raster file or what is wrong with the raster I am creating, since it seems like all I have to do to run this code is simply entering my input shapefile, which it seems like I did correctly.

UPDATE (with second attempt at code, incorporating pixel size and setting the projection):

from osgeo import gdal, ogr

# Define NoData value of new raster
NoData_value = -9999

# Filename of input OGR file
vector_fn = 'C:/Users/MyName/Desktop/NYC_GIS/Boroughs_Test/geo_export_d92842a7-2631-4fe5-ab83-ee070743980c.shp'

# Filename of the raster Tiff that will be created
raster_fn = 'C:/Users/MyName/Desktop/NYC_GIS/Boroughs_Test/test2.tif'

# Open the data source and read in the extent
source_ds = ogr.Open(vector_fn)
source_layer = source_ds.GetLayer()
proj = source_layer.GetSpatialRef()
x_min, x_max, y_min, y_max = source_layer.GetExtent()

# Create the destination data source
target_ds = gdal.GetDriverByName('GTiff').Create(raster_fn, 100, 100, 1, gdal.GDT_Byte)
target_ds.SetProjection(proj.ExportToWkt())
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])

This code runs and outputs a 0. I am not sure what this means. I then go to open test2.tif in QGIS, and all I see is a black square with all 0 values, not the rasterized polygon shapes I expect.

Answer 1

https://github.com/corteva/geocube

from geocube.api.core import make_geocube

vector_fn = 'C:/Users/MyName/Desktop/NYC_GIS/Boroughs_Test/geo_export_d92842a7-2631-4fe5-ab83-ee070743980c.shp'

out_grid = make_geocube(
    vector_data=vector_fn,
    measurements=["column_name"],
    resolution=(-25, 25),
    fill=-9999,
)
out_grid["column_name"].rio.to_raster("my_rasterized_column.tif")

Answer 2

You are telling GDAL to burn all 0 values with burn_values=[0].

Additionally, setting xRes and yRes to 100 won't work as your data uses lon, lat coords (1 degree is 111,000m at the equator so 100 degrees is huuuuge).

By the way, the newer gdal.Rasterize is easier to use than the old gdal.RasterizeLayer.

Particularly there's no need to pre-create the output dataset, just pass an output filename. For example:

from osgeo import gdal

# Define NoData value of new raster
NoData_value = -9999

# Filename of input OGR file
vector_fn = 'Boroughs.geojson'

# Filename of the raster Tiff that will be created
raster_fn = 'Boroughs.tif'

# Open the data source and read in the extent
source_ds = gdal.OpenEx(vector_fn)
pixel_size = 0.001  # about 11 metres(ish) use 0.001 if you want roughly 100m 

gdal.Rasterize(raster_fn, source_ds, format='GTIFF', outputType=gdal.GDT_Byte, creationOptions=["COMPRESS=DEFLATE"], noData=NoData_value, initValues=NoData_value, xRes=pixel_size, yRes=-pixel_size, allTouched=True, burnValues=1)

Answer 3

The parameters you have entered into the Create() function are incorrect. You can find them using the python help function:

from osgeo import gdal
driver = gdal.GetDriverByName('GTiff')
help(driver.Create)
Create(Driver self, char const * utf8_path, int xsize, int ysize, int bands=1, GDALDataType eType, char ** options=None)

In your case, you have entered the pixel resolution (xres, yres) inplace of the raster dimensions (xsize, ysize):

target_ds = gdal.GetDriverByName('GTiff').Create(raster_fn, x_res, y_res, 1, gdal.GDT_Byte)

To avoid having to manually calculate the raster dimensions, you could rasterise using a reference image as illustrated here.

You have also failed to define a projection system. Consequently, the output image will not be georeferenced. You can do this using the SetProjection() function as follows:

vector = ogr.Open(vector_fn)
layer = vector.GetLayer()
proj = layer.GetSpatialRef()
target_ds.SetProjection(proj.ExportToWkt())

Your burn value is set to zero. You should probably change this to 1:

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

Answer 4

I still use gdal 1.11 version and gdal_rastersize with this issue actually, it will need to have the extent of the shapefile (can be retrieved by python gdal):

    geom = feature.GetGeometryRef()
    extent = geom.GetEnvelope()

Then, gdal_rastersize command:

gdal_rasterize -burn 1 -a_nodata 0 -ot Byte -tr 0.025 0.025  polygon0.shp boroughs.tiff -a_srs EPSG:4326 -te 16.2475 43.6218 17.1124 44.3496