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?
Best Answer
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:
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.