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.

• Michael Stimson: Related gis.stackexchange.com/questions/390695/… The use of os.system should be avoided, look up subprocess.Popen as an alternative to os.system docs.python.org/3/library/subprocess.html#replacing-os-system but even better to use the GDAL object than to shell it.

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)

• LostinSpatialAnalysis: Apologies this took so long, but similar to your suggestion on another one of my posts, I tried this out, and it is working, but it is taking very long to run. A raster file is produced still, but it is 1,470 megabytes big, which seems to be too much for my simple NYC example. QGIS crashed when I try to open and view this output raster. The python code never actually finished. Is this based on how the resolution or pixel size was set up?