Mask raster with shapefile in Python using GDAL/OGR

Question

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).

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'])

Answer 1

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)

Answer 2

The reason a polygon is required (rather than a line), is that a line can divide the raster, but how would gdal know which side of the line you wanted to keep?

You can close your line to make a polygon as suggested by Stefan; there is another approach if you want to do it more algorithmically:

Step 1. Burn your polyline into a raster with the same number of rows, columns, and extents using gdal_rasterize. You can get the resolution, etc. with gdal bindings in Python.

Step 2. Load your rasterized line tiff into a numpy array. Invert the array (so the line pixels are 0, all other pixels are 1). Use scipy's label function to label the areas on each side of the burned-in line; there will be two labels, one corresponding to each side.

Step 3. If you don't want to have to manually choose which side of the line to keep, you will need a metric that differentiates between the ocean and land pixels. Looking at your raster, that shouldn't be a problem. Something like the average pixel value within each of the two areas (I'm not sure what your raster is showing) should suffice. E.g. avg(ocean) is always less than avg(land), so select the side (i.e. label) with the smaller average value.

Step 4. Now you have a raster mask for the pixels you want--i.e. those with the label corresponding to the side you automatically chose with your index in Step 3. You can include the line pixels, too, if you want.

Might be overkill, but if you have dozens or more images to process it might be worth it. There may also be problems with this approach if your shoreline leaves an image and then re-enters it elsewhere; in such cases you will have more labels and will have to threshold the avg pixel value to determine which labeled areas to keep.