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'])
  • I think gdalwarp needs a polygon as input. You could try to create a polygon from your linestring (coast line). For that you could extent the line in this way that the resulting polygon covers the land area of your raster. Usage: gis.stackexchange.com/questions/211748/… – Stefan Mar 14 '18 at 11:05

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)
| improve this answer | |
  • gdal_rasterize does work with lines. gdal.org/gdal_rasterize.html – Jon Mar 15 '18 at 19:22
  • @Jon, yes, I will edit my answer to be more clear. It will work with lines, but it will only burn value for intersection with the line, which is not what Damuno wants. – Justin Braaten Mar 15 '18 at 19:27

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.

| improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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