2

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.

1

1 Answer 1

4

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)
1
  • 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? Jan 10, 2022 at 21:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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