With ArcGIS zonal statistics as a table I can easily calculate the mean value from a raster contained within a polygon from a shapefile. I now have a Mac, and therefore I cannot access ArcGIS anymore. Actually, I am trying to use open source software only (e.g. R and QGIS).

I thought that "zonal statistics" from QGIS would do the trick, however, an essential feature for my needs is missing: in ArcGIS, you can specify the column by which you want to aggregate the values. For example, if in your shapefile you have 100 polygons, you may not want to just calculate the mean value from the raster in each polygon. Rather, you may want to know the mean value of a group of polygons. In my case, one column ("veg_class") classifies all polygons as "tropical", "temperate" or "boreal" vegetation.

How can I, in QGIS or any other software, select the column to aggregate the statistics so, in the end, I obtain one mean per veg_class (i.e. 3) instead of per polygon?


You can do this in R with the help of raster and rgdal library



Load the shapefile and raster

shp <- readOGR("shape file directory", "shapegile_name")
in.ras <- raster("raster file")

Extract the raster value over the shapefile val <- extract(in.ras, shp)

Create a data frame to store the extracted value and the corresponding pixels group. df <- data.frame(matrix(NA, 0 , 2)) n <- length(val) for (i in 1:n){ v = val[i] m <- length(v) for ( j in 1:m){ df <- rbind(df, c(shp@data$group[i], j)) } }

Find the mean by group wise

df2 <- aggregate(df$val, by=list(df$group), mean)

Later you can save the data frame as csv or shapefile. It will be good if you find some other way instead of for loop.

| improve this answer | |
  • thanks, I appreciate your help. The extraction takes forever and finally the computer crashes before completion. I don't know if removing all unnecessary columns (all except the grouping column) in the shp (27Mb) would accelerate the process, although I don't know if that's even possible or how to do it. For the moment, the extraction process in R seems not suited for my purpose in terms of efficiency. – fede_luppi Oct 1 '16 at 15:43
  • To remove the other other columns in shape file, you may use this shp@data = shp@data$column_name. Then save this as separate, shape file using writeOGR. I don't know whether this will be reduce the time for extraction. I never tried extract over large area in R. – S. Thiyaku Oct 2 '16 at 4:23
  • Thanks, but it didn't help. It just seems R is not very efficient running this process so I need an alternative to extract – fede_luppi Oct 3 '16 at 8:56
  • 1
    If your polygon not covering the majority of raster, try to crop the raster before extract (crop(raster, shp)) or Check the zonal statistics code using python or in qgis zonal statistics will output count and mean, so do sum(mean * count)/sum(count) for every group in field calculator or qgis python console. I hope somebody else will come with better and easier solution – S. Thiyaku Oct 3 '16 at 13:09

Not a direct solution but a possible (and possibly tedious) workaround is to:

  1. Use the Split vector layer tool on your "veg_class" column which creates a shapefile for each classification.

  2. Then run the Zonal Statistics plugin (which allows you to choose which statistics you want calculated) or the tool from the Processing Toolbox or Raster menu on each vector layer.

  3. (Optional) Then merge the outputs of the zonal statistics into one shapefile using the Merge vector layers from the Processing Toolbox or Vector menu.

| improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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