Mask image based on no-data pixels in another image with different projection - GDAL or Python

Question

I've got a land cover image which covers the whole globe, and a satellite image which covers a very small part of the globe. These two images have different projections, co-ordinate systems and pixel sizes. Also, the satellite image is not fully 'square' within the image, it is on a slant, and there are significant no-data areas around the edge of the actual image.

I want to extract the pixels in the land cover image which intersect the satellite image.

All of this needs to be completely automatic - either via a Python script or a GDAL command-line. In the past when I've had a polygon outline of the image I've used a gdalwarp command like

gdalwarp -overwrite -of GTiff -cutline polygon.shp -crop_to_cutline landcover.tif output.tif

But I can't find a good automatic way to produce a polygon of the footprint of my satellite image, given than I don't want the footprint of the whole image, but just the part of the image that isn't No Data. So it doesn't look like that method will work. Theoretically I should be able to do a simple mask, but I can't work out how to do this when my images are of different resolutions.

Does anyone have any ideas how I should go about doing this?

Answer 1

You could create a mask from the image, polygonize this and then use the same gdalwarp command you normally use. If you wanted to do this in Python you could use the bandmath function in RSGISLib (rsgislib.org). For example:

# Import RSGISLib
import rsgislib
from rsgislib import imagecalc
from rsgislib.imagecalc import BandDefn
import subprocess

# 1) Create mask
inImage = 'input_image.tiff'
outMask = 'image_mask.tiff'

gdalFormat = 'GTiff'
dataType = rsgislib.TYPE_8UINT
expression = 'b1 > 0?1:0' # Only include values greater than zero.
bandDefns = []
bandDefns.append(BandDefn('b1', inImage, 1)) # Use band one (variable b1) of in image
imagecalc.bandMath(outMask, expression, gdalFormat, dataType, bandDefns)

# 2) Generate polygon from mask
subprocess.call("gdal_polygonize.py -mask image_mask.tiff image_mask.tiff -b 1 -f 'ESRI Shapefile' image_Footprint.shp",shell=True)

There are other ways you could create an image mask if you didn't want to use RSGISLib, you could use RIOS (https://bitbucket.org/chchrsc/rios/wiki/applier/simpleexample) to accomplish the same thing.

Answer 2

This is a great question. I typed up a possible answer (with pictures) in an IPython Notebook here: http://nbviewer.ipython.org/gist/sgillies/9713809. Rasterio, the better Python interface for GDAL, is used to get the mask, polygonize it, and write out a working GeoJSON cutline.