I'm new to QGIS (2.18.11). I've got OS Terrain 50 elevation data, in around 150 ASC tiles (10 km x 10 km) covering my whole area of interest (north-east Scotland). There is just a single band. I've got a grid of 1 km squares for my area of interest (made using TomBio plugin).

I'd like to extract mean elevation values per 1 km square (I may look for further stats later but mean is fine for starters). I've used Zonal Statistics to do this on one tile at a time and it adds a field to the Table of Attributes of the TomBio grid. That's all fine, and just what I need - except, is there a quicker way to do it for all 150 tiles at the same time?

2 Answers 2


You could make a virtual raster (vrt) which will be quicker than merging them all (which is another option but will result in a huge file - so you will need to probably set the merge process to output a BigTiff). A VRT is a much lighter-weight solution for your needs than a merge. To build a VRT go Raster->Miscellaneous->Build Raster Catalogue Then do your stats on the vrt. For more information on the VRT see the documentation (this will let you tweak the settings)


Using R.


#set QGIS directory
qgis_env <-set_env(root ='C:/Program Files/QGIS 2.18') 

-Function and parametrs.

#entrada = folder with all rasters in .tif
#saida = folder witch zonal polygon will be dropped
#bacia = polygons

zonalTop = function(entrada, saida, bacia){

  nomesRasters = list.files(entrada, pattern="*.tif$")
  lerRasterCompl = list.files(entrada, pattern="*.tif$", full.names = TRUE)
  rasters = lapply(lerRasterCompl, raster)

  for (i in 1:length(rasters)){

    run_qgis(alg = "qgis:zonalstatistics",
             INPUT_RASTER = rasters[[i]],
             RASTER_BAND = 1,
             INPUT_VECTOR = bacia,
             GLOBAL_EXTENT = 1,
             OUTPUT_LAYER = paste(saida, tools::file_path_sans_ext(nomesRasters[[i]]),".shp", sep = ""))

  nomesDbf = list.files(saida, pattern="*.dbf$")
  lerDbfCompl = list.files(saida, pattern="*.dbf$", full.names = TRUE)
  Dbfs = lapply(lerDbfCompl, read.dbf)

  for (i in 1:length(Dbfs)){

    names(Dbfs[[i]]) = paste(tools::file_path_sans_ext(nomesDbf[i]), names(Dbfs[[i]]),  sep="_")


  Dbfs = clean_names(do.call("cbind", Dbfs))



- Runnin Function Example

anch10k = zonalTop("E:/artigo_declividade/morfomometriatif/anchi/",

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.