[GIS] Mask raster with shapefile in Python using GDAL/OGR

gdalogrpythonrastershapefile

I'm building a script in Python to mask out land and keep the ocean. As seen in the following I have a shapefile for the coastlines (one feature only).

enter image description here

I want to mask out the the land to black or null in Python using GDAL/OGR.

I'm fairly new in using the gdal librairies and I've tried the following without luck:

    inraster = 'raster.tif'
    inshape = 'shape.shp'
    subprocess.call(['gdalwarp', inraster, outraster, '-cutline', inshape,
                 '-crop_to_cutline', '-cwhere'])

Best Answer

The gdal_rasterize utility can be used for masking too. In your case, make sure you are using a polygon, not a line, otherwise only pixels intersecting the line will be masked. Here is an example:

# get subprocess to call GDAL util
import subprocess

# define paths to raster and vector
inraster = 'raster.tif'
inshape = 'shape.shp'

# make gdal_rasterize command - will burn value 0 to raster where polygon intersects 
cmd = 'gdal_rasterize -burn 0 '+ inshape + ' ' + inraster

# run command
subprocess.call(cmd, shell=True)
Related Question