How to rasterize polygons using GDAL commands in python?

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 would like to use a GDAL command line command to accomplish this, but I would like to carry this out in Jupyter Notebook using python, rather than through the OSGeo4W Shell.

I am calling this input shapefile boroughs.shp

This is what I am trying to use after navigating to the folder boroughs.shp is located in:

from osgeo import gdal
from osgeo import ogr
import numpy as np
import matplotlib.pyplot as plt
import os

gdal.UseExceptions()

os.system([gdal_rasterize -b1 -burn 1 -a_nodata 0 -ot Byte -tr 25 25 -l boroughs boroughs.shp boroughs.tif])

And this returns the following error:

  File "<ipython-input-7-921412e61b57>", line 9
    os.system([gdal_rasterize -b1 -burn 1 -a_nodata 0 -ot Byte -tr 25 25 -l boroughs boroughs.shp boroughs.tif])
                                        ^
SyntaxError: invalid syntax

Is there a simple fix to my syntax so that the burning option is correct, or is this approach to using a GDAL command using python in Jupyter Notebook not the correct approach at all? My goal is to simply output a raster file of the NYC borough polygons with 1 burned into where the polygons are, and 0 outside of the polygons, with the same extent as the input shapefile.

Answer 1

Your syntax is invalid because you need to pass a list of "strings" to os.system (note: you should always use functions from the subprocess module in preference to os.system).

You're also trying to create enormous 25x25 degree pixels. Either reproject your vector data or use a rough decimal to meter conversion (0.00025)

You're better off using gdal.Rasterize (python equivalent of gdal_rasterize):

from osgeo import gdal

# Define NoData value of new raster
NoData_value = -9999

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

# 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.00025  # about 25 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)