[GIS] How to get extent of raster in QGIS Modeler

gdalqgisqgis-modeler

In my journey to solve a different question, I've decided to create a QGIS processing model to automatically correlate rasters to a projection (UTM) grid.

The inputs are a raster file and the UTM grid (shapefile).

The process is: get raster extent as a polygon, do a join by location with the raster extent on the UTM grid, get the correct zone, and pipe it into GDAL Warp method to output a correctly projected raster.

There is a function in the processing modeler called Raster layer bounds that I can use on the raster input, but cannot use as an input to anything else.

My question is: how do I get from the Raster layer bounds to a polygon in modeler? I see many tools for making a polygon into other things, but nothing to make a series of points into a polygon.

Perhaps I am going about this wrong and there is something that will help me extract the extent of the raster file more easily, but searching by "extent" yields nothing. There appears to be a very helpful Polygon from layer extent function, but alas, it is only for vector not raster.

Best Answer

You could create a custom script, copy the source code for the Polygon from layer extent tool and modify it slightly to take a raster as input. You can create one from:

Processing Toolbox > Scripts > Tools > Create new script

Then use something like the following:

##Example=name
##Layer=raster
##Output=output vector

from qgis.core import QgsField, QgsPoint, QgsGeometry, QgsFeature, QGis
from PyQt4.QtCore import QVariant

layer = processing.getObject(Layer)
fields = [
    QgsField('MINX', QVariant.Double),
    QgsField('MINY', QVariant.Double),
    QgsField('MAXX', QVariant.Double),
    QgsField('MAXY', QVariant.Double),
    QgsField('CNTX', QVariant.Double),
    QgsField('CNTY', QVariant.Double),
    QgsField('AREA', QVariant.Double),
    QgsField('PERIM', QVariant.Double),
    QgsField('HEIGHT', QVariant.Double),
    QgsField('WIDTH', QVariant.Double),
]

rect = layer.extent()
minx = rect.xMinimum()
miny = rect.yMinimum()
maxx = rect.xMaximum()
maxy = rect.yMaximum()
height = rect.height()
width = rect.width()
cntx = minx + width / 2.0
cnty = miny + height / 2.0
area = width * height
perim = 2 * width + 2 * height

rect = [QgsPoint(minx, miny), QgsPoint(minx, maxy), QgsPoint(maxx, maxy), QgsPoint(maxx, miny), QgsPoint(minx, miny)]
geometry = QgsGeometry().fromPolygon([rect])
feat = QgsFeature()
feat.setGeometry(geometry)
attrs = [minx, miny, maxx, maxy, cntx, cnty, area, perim, height, width]

writer = processing.VectorWriter(Output, None, fields, QGis.WKBPolygon, layer.crs())
feat.setAttributes(attrs)
writer.addFeature(feat)
del writer

Make sure to save the script into your /.qgis2/processing/scripts directory.

Example:

  1. A simple model with a raster layer being used as an input to the script:

    Example model

  1. Here is the input raster:

    Raster

  1. Here is the result of the model:

    Result

You could take the output polygon and use that as input for another tool.

Related Question