[GIS] Improving Clip Raster workflow in QGIS

clipqgis-2raster

I've tried to come up with an open-source equivalent to the following process in ArcGIS 10.3:

  1. Spatial Analyst Toolbox > Extract by Mask
  2. target input: 1" SRTM DEM in EPSG 4019
  3. mask layer: single rectangular polygon, shapefile in EPSG 4283
  4. Environment Settings: Snap Raster to target input, output coordinate system EPSG 4283

This gets me a 'perfect' clip – no changes to pixel size or alignment, identical cell values in the output. The output doesn't align perfectly to the clipping polygon of course, but I don't need it to. Arc just extracts all the cells intersecting the polygon (although I note that the top-most and right-most columns are set to no-data, presumably because their upper-right corners were outside the bounding box).

I've worked out one solution using GRASS 7 and GDAL processing tools in QGIS 2.14.3:

  1. Create a mapset and import the raster to be clipped
  2. Run g.region.multiple.raster so that region settings match the imported dataset
  3. Use v.in.ogr to import the bounding box shapefile
  4. Use v.to.rast.constant to convert the bounding box to a raster with value 1 for all cells. The output is aligned with the input due to the region settings.
  5. Use g.region.zoom with the new raster boundary to shrink the region without altering cell size or alignment
  6. use r.out.gdal.geotiff to export the source dataset. Only the area within the updated region is exported.
  7. Use gdalwarp to re-project the output .tif in EPSG 4283

The output is nearly identical to the ESRI workflow, but where ESRI set the top-most and right-most columns to no-data, GRASS just didn't export them. Fine by me.

Can an equivalent procedure be run in QGIS without using GRASS?

Best Answer

Ok, having learnt a lot of stuff in the months since I asked this question, here are a couple of options:

In R, very similar to Bastien's answer:

library(raster)

rootdir <- 'C:/Users/obrl_soil/Downloads'
setwd(rootdir)
tinyraster <- raster('tinyrast.tif')
tinypolygon <- readOGR('tinypolygon.shp')

# alter your mask polygon to line up with the nearest pixel edges
tpoly_aligned <- alignExtent(tinypolygon, tinyraster, snap='near')

# clip and export in one hit
tinyclip <- crop(tinyraster, tpoly_aligned, filename='tinyclip.tif')

I noticed when checking the outputs in QGIS that writeRaster introduced some spurious data in the output tif, but crop didn't - e.g. input elevation for a pixel was 108.10545, output with crop was the same, but output with writeRaster was 108.10545349121094.

Alternate workflow - get the polygon bounding box coordinates and feed them into gdal_translate using -projwin. Just be careful which version of gdal_translate you use!. 1.11 is fine, and 2.1.2 should be when its released.

You can also do this in R, like

gtrans111 <- 'C:/Program Files/GDAL/gdal_translate.exe'
tinypolygon <- readOGR('tinypolygon.shp')
tpbb <- toString(c(tinypolygon@bbox[1], tinypolygon@bbox[4],
                   tinypolygon@bbox[3], tinypolygon@bbox[2]))
tpbb <- gsub(', ', ' ', tpbb)

system2(gtrans111, args= c('-projwin', tpbb, 'tinyraster.tif', 'tinyclip2.tif' ))

or even bypass creating any R objects by just inputting the polygon bounding coordinates directly 

system2(gtrans111, args= c('-projwin', '148.665 -20.88 148.67 -20.884', 
                           'tinyraster.tif', 'tinyclip3.tif'))

The output extent is not identical to the raster package methods, but does align correctly. All of these methods are easy to loop across multiple datasets, which is the real advantage of R. You can do something similar in Python too.

Related Question