Calculating Image boundary / footprint of satellite images using open source tools?

Question

I need to create the polygon outlines of several single layer raster images, not the extent/bounding box, but the area without the nodata values, as shown here: Creating shapefile showing footprints of Rasters?.

In the answer to the question above, the Image Boundary plugin was mentioned, however, I don't find it in the QGIS 1.8.0 Lisboa, running on Ubuntu.

Is the tool still available?

If not, is there a way to do this using opensource tools, either gdal, R, QGIS, GRASS, or similar, not with ArcMap?

Answer 1

I finally found a way to do this:

step 1: gdalwarp -dstnodata 0 -dstalpha -of GTiff foo1 foo2

This does two important things: it sets the destination No Data (outside border) values to 0, and it creates an alpha band.

step 2: gdal_polygonize.py foo2 -b 2 -f "ESRI Shapefile" foo3

The second step uses the alpha band (band 2), created in step 1, and creates a shapefile from that band.

This can then easily be scripted in a bash script, if you have many images, to create exact outlines for.  

Answer 2

Image Boundary plugin did not work for me either, therefore I used the same approach with GDAL. Nevertheless it only worked for me after changing the first step to:

step 1: gdalwarp -srcnodata 0 -dstalpha -of GTiff foo1 foo2

I am working with Landsat8 band (where no data=0) and when using the -dstnodata function I get:

band1 with no data = 'no data'
band2 (Alpha band) = '255' for the entire scene/extent

whereas with -srcnodatafunction I get:

band1 with no data = 'no data'
band2 (alpha band) with no data = 'no data' and valid data = '255' which then allows to extract polygon for valid data only.

I couldn't fully understand the reason behind this behavior (how alpha is computed?) , but I hope this might help others facing the same problem.

Answer 3

I used gdal_translate as suggested by the GDAL project.

gdal_translate -b mask -of vrt -a_nodata 0 test.tif test.vrt
# Note the  -a_nodata 0 doesn't seem to work when the mask is input, so do another pass
gdal_translate -b 1 -of vrt -a_nodata 0 test.vrt test2.vrt
gdal_polygonize.py -q  -8 test2.vrt -b 1 -f "ESRI Shapefile" testdata.shp

Answer 4

You can add the old repository where have plugins out in Official Repository.

http://pyqgis.org/repo/contributed

The image boundary plugin have option for calculate valid pixel, but, the image need be the full scene, how CBERS or Landsat, where this process calculate the first 4 corners. The image need have nodata with ZERO value for area without imaging(scanned by satellite sensor).

Author of Image Boundary

Answer 5

The Image Boundary stayed obsolete. You can use the Footprint plugin. http://pt.slideshare.net/LuizMotta3/qgis-ibama-imagefootprint

Answer 6

gdal_footprint, which came out in 3.8.0 (Nov 13, 2023), should do this.

Example:

gdal_footprint -t_srs EPSG:4326 input.tif output.geojson

Answer 7

It is also possible to use rasterio and extract the shapes of all connected pixel outlines:

from shapely.geometry import shape
import geopandas as gpd
import rasterio
from rasterio import features

rasterfile= r"/home/bera/Desktop/gistest/four_polygons_rasterized.tif"
dataset = rasterio.open(rasterfile, mode="r")
arr = dataset.read()[0,:,:]

#If your input raster have many values, each unique value will become one polygon.
#Reclassify to either 1 or 0. You dont have to do this.
arr[(arr>0)]=1
arr[(arr<1)]=0

#Vectorize each raster blob
mask = arr>0 #Dont create polygons where the value is <=0
shapes = rasterio.features.shapes(source=arr, connectivity=4, mask=mask,
                                              transform=dataset.transform)

#Iterate over the shapes, each row is a tuple of the geometry and the raster value:
#({'type': 'Polygon', 'coordinates': [[(618821.3726, 6480531.5759), ....), 1.0)

rastervalues = [] #A list for the values
geometries = [] #And geometries
for geom, rasterval in shapes:
    geometries.append(shape(geom)) #Create a shapely polygon
    rastervalues.append(rasterval)

df = gpd.GeoDataFrame(data=rastervalues, geometry=geometries, 
                      crs=dataset.crs.to_epsg(), columns=["rasterval"])
#df.to_file(r"/home/bera/Desktop/gistest/raster_boundaries.gpkg")