8

Comparing results of zonal statistics in QGIS and R, I noticed that the results they both yield do somehow differ. Manually checking, I also found that both QGIS and R seem to partly produce erroneous results. I wonder whether there is a robust solution to calculating mean values of raster cells within polygons.

My use case: I have a Raster Layer with cell size = 1 x 1 km² and a Polygon Layer with quadratic polygons of size = 250 x 250 m². Both are projected in EPSG:5677.

I used raster -> zonal stats with the "mean" option in QGIS ("QGIS.Result") and tried

R.Result1<-extract(Raster, Polygons, method="simple", fun=mean)

as well as

R.Result2<-extract(Raster, Polygons, method="simple", fun=mean, weights=TRUE, normalizeWeights=TRUE)

in R with the raster package. The results look like this:

ID   QGIS.Result  R.Result1    R.Result2
1   -63.9505964  -64.0000000  -64.60000000
2   -67.0000000  -67.0000000  -67.00000000
3   -67.0000000  -67.0000000  -67.00000000
4   -62.3765259  -64.0000000  -61.85714286
...  ...          ...          ...
19  -99.9498123  -97.5999908  -99.19999186
20  -92.7999878  -92.7999878  -92.79998779
21  -92.7999878  -92.7999878  -92.79998779
22  -102.5493486 -92.2999954  -99.19999186
23  -97.5452805  -89.3999939  -92.79998779
24  -97.5452805  -89.3999939  -92.79998779
25  -87.4895458  -87.0000000  -87.33333333
26  -86.0000000  -86.0000000  -86.00000000
27  -86.0000000  -86.0000000  -86.00000000
...  ...          ...          ...

Results differ whenever a Polygon intersects with more than one raster cell - note in the example below that:

  • all three methods yield the same result for polygons 20,21,26,27 which all do intersect with only one raster cell and different results for the other polygons which do intersect with more than one raster cell.

  • the QGIS result gives a smaller value for polygon 22 than for polygon 19 although it should be the other way round. R.Result2 gives the same result for both polygons. Only R.Result1 performs better here.

  • polygon 25 intersects with Raster values 88 and 86. R.Result1 is the mean of both: 87. However, as polygon 25 contains a larger amount of value 88 than of value 86, the result should be biased towards 88 (like in QGIS.Result and R.Result2).

enter image description here

Thus, all three methods seem not to provide what I need: a robust, coverage-based mean of the raster layer for each polygon. Is there a way to obtain such (preferably in R)?

The used data can be found here: https://www.dropbox.com/s/bctjidc4dz2j396/testdata.7z?dl=0

  • R.Result1 and R.Result2 have the same expression, is this a mistake? – Albert Oct 24 '17 at 14:08
  • thank you, @Albert. Yes it's been a mistake, I corrected it. – yenats Oct 26 '17 at 7:26
  • I tried to install your exactextractr package on R 3.4.4 (Debian jessie). I compiled it on Rstudio with the command mentioned on Github (devtools::install_github('isciences/exactextractr') It outputs the following error and stops : * installing source package ‘exactextractr’ ... ** libs Error in .shlib_internal(args) : C++14 standard requested but CXX14 is not defined Is there a way to overcome this error ? For information, gcc 4.9.2 is installed on debian 8. Thanks in advance – Olive gak Aug 22 '18 at 18:46
  • If you have a new question, please ask it by clicking the Ask Question button. Include a link to this question if it helps provide context. - From Review – Dan C Aug 22 '18 at 19:49
4

A Robust Zonal stats in QGIS can be implemented with PyQGIS. Following code was run with your layers by using a filter to select only ID_sub.pat between 19 and 27 (showed in your image).

import processing

registry = QgsMapLayerRegistry.instance()

Polygons = registry.mapLayersByName('Polygons')
request = QgsFeatureRequest().setFilterExpression (u'"ID_sub.pat" >= 19 AND "ID_sub.pat" <= 27')

iterator = Polygons[0].getFeatures(request)

feats_Polygons = [ feat for feat in iterator ]

Raster = registry.mapLayersByName('Raster')
rprovider = Raster[0].dataProvider()

xmin, ymin, xmax, ymax = Raster[0].extent().toRectF().getCoords()

extent = str(xmin)+ ',' + str(xmax)+ ',' +str(ymin)+ ',' +str(ymax) 

parameters = {"TYPE":1, 
              "EXTENT":extent, 
              "HSPACING":1000, 
              "VSPACING":1000, 
              "CRS":Raster[0].crs().authid(), 
              "OUTPUT":None}

path = processing.runalg('qgis:creategrid', 
                         parameters)

grid = QgsVectorLayer(path['OUTPUT'],
                      'grid',
                      'ogr')

feats_grid = [ feat for feat in grid.getFeatures() ]

for featg in feats_grid:
    for featp in feats_Polygons:
        if featg.geometry().intersects(featp.geometry()):
            geom = featg.geometry().intersection(featp.geometry())
            pt = geom.centroid().asPoint()
            value = rprovider.identify(pt,
                                       QgsRaster.IdentifyFormatValue).results()[1]

            print value, geom.area()/62500.0, featp.id()+1

After running code at Python Console of QGIS it was obtained this result:

-88.0 0.11367261864 22
-88.0 0.744772910144 25
-86.0 0.0389546011361 22
-86.0 0.152627219776 23
-86.0 0.152627219776 24
-86.0 0.255227089856 25
-86.0 1.0 26
-86.0 1.0 27
-102.399993896 0.744772910144 19
-102.399993896 0.631100291505 22
-92.799987793 0.255227089856 19
-92.799987793 1.0 20
-92.799987793 1.0 21
-92.799987793 0.21627248872 22
-92.799987793 0.847372780224 23
-92.799987793 0.847372780224 24  

where at first column are raster values, at second column weight factor and at third column id Polygons feature. You can observe that for polygons 20,21,26,27, its factors is always one (as expected). In other cases is proportional to intersection of Polygons feature with grid representing raster.

Results were corroborated with Value Tool and QuickWKT QGIS plugins and they were as expected.

4

QGIS and R's raster package use different methods to estimate zonal statistics. Briefly:

QGIS compares the centroid of each raster cell to the polygon boundary, initially considering cells to be wholly within or outside of the polygon based on the centroid. However, if fewer than two cell centroids fall within the polygon, an exact vector-based calculation is performed instead. That calculation is essentially the same as xunilk's answer.

R's raster package also uses a centroid test to determine if cells are inside or outside of the polygon. When weights=TRUE, it disaggregates the raster by a factor of 10 before performing the analysis, which reduces the error incurred by ignoring partially covered cells but reduces performance substantially.

If you're working in R, I would recommend the exactextractr package (I am the author). It provides a replacement for raster::extract, using an exact vector-based calculation that is hundreds of times faster than the original in most cases.

Here's an example using the test data provided:

library(sf)
library(raster)
library(exactextractr)

polygons <- st_read('testdata/Polygons.shp')
grid <- raster('testdata/Raster.tif')

polygons$result <- exact_extract(grid, polygons, fun=weighted.mean)

example <- polygons[c(1:4, 19:27), 'result']
st_geometry(example) <- NULL
print(example)

#       result
# 1  -63.95060
# 2  -67.00000
# 3  -67.00000
# 4  -62.37653
# 19 -99.94981
# 20 -92.79999
# 21 -92.79999
# 22 -98.04804
# 23 -91.76212
# 24 -91.76212
# 25 -87.48955
# 26 -86.00000
# 27 -86.00000
  • however I try, I cannot install exactextractr on my windows machine: 'ERROR: compilation failed for package 'exactextractr'' – yenats Aug 28 '18 at 20:00
  • @yenats Do you have the R build tools installed? (cran.r-project.org/bin/windows/Rtools) Feel free to post a GitHub issue with the full output you're seeing from devtools::install_github. – dbaston Aug 28 '18 at 20:11
  • yes, I have the R build tools installed. 'devtools::install_github('isciences/exactextractr')' would also install several packages; however at the end of the process it rolls everything back. – yenats Aug 29 '18 at 10:08
  • @yenats I'd need to see the full command output to be helpful. Feel free to email it to me (contact info in profile) if you don't want to post it on GitHub. – dbaston Aug 29 '18 at 19:44

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.