[GIS] Get Raster Values from a Polygon Overlay in Opensource GIS Solutions

polygonqgisrrasterstatistics

I have two layers. A polygon-shape-layer with many tiles and a raster-layer containing CORINE 2006 land cover with many categories in a colourmap.
I want to obtain for each polygon in the shapelayer a sum of each land cover category of the raster-layer.

For example there is a polygon with id '2' and i want to Attributes like this for this polygon (in percent or square meters):

Arable land : 15 %
Forest: 11 %
Streets:2 %
(… and so one)

I tried to do it in grass, qgis (no function), saga (just sums up every to a total value) r(total sum), but i still found no solution. Most plugins (zonal statistics in qgis) only support 0-1 raster layers. v.rast.stats didn't help either. Iam open to any good and smart solution!. Maybe i even used a wrong approach or made mistakes.

In Arcgis this task is quite easy, if am remember right, but i am still missing a good solution for your everyday linux user.

I am running a debian linux system and this why i can only use programs for this OS.

EDIT:
Because this question still has so many views and visitors:
I wrote a QGIS-plugin, which also is capable of calculating the landcover of raster layer. I have'nt coded a polygon overlay yet, but it definitely planed. Find the plugin here and install the Scipy library first.

Best Answer

Use 'extract' to overlay polygon features from a SpatialPolygonsDataFrame (which can be read from a shapefile using maptools:readShapeSpatial) on a raster, then use 'table' to summarise. Example:

> c=raster("cumbria.tif") # this is my CORINE land use raster
> summary(spd)
Object of class SpatialPolygonsDataFrame
[...]
> nrow(spd)  # how many polygons:
[1] 2
> ovR = extract(c,spd)
> str(ovR)
List of 2
 $ : num [1:542] 26 26 26 26 26 26 26 26 26 26 ...
 $ : num [1:958] 18 18 18 18 18 18 18 18 18 18 ...

So my first polygon covers 542 pixels, and my second covers 958. I can summarise each of them:

> lapply(ovR,table)
[[1]]

 26  27 
287 255 

[[2]]

  2  11  18 
 67  99 792

So my first polygon is 287 pixels of class 26, and 255 pixels of class 27. Easy enough to sum and divide and multiply by 100 to get percentages.

Related Solutions

[GIS] combine raster and polygon values in R

The following R example essentially performs what ArcGIS users call Zonal Statistics. This should be a good building block for your analysis. The main function performing this analysis is extract() from the raster package.

require(raster)

# Create some sample raster data
r <- raster(ncol=36, nrow=18)
r[] <- 1:ncell(r)
plot(r)

#Create some sample polygons
cds1 <- rbind(c(-180,-20), c(-160,5), c(-60, 0), c(-160,-60), c(-180,-20))
cds2 <- rbind(c(80,0), c(100,60), c(120,0), c(120,-55), c(80,0))
polys <- SpatialPolygons(list(Polygons(list(Polygon(cds1)), 1), 
                              Polygons(list(Polygon(cds2)), 2)))
plot(polys)

# Extract the raster values underlying the polygons
v <- extract(r, polys)
v

# simplify to display mean values
output = unlist(lapply(v, function(x) if (!is.null(x)) mean(x, na.rm=TRUE) else NA ))
print(output)

Edit:

Here is a simple calculation of the zonal mean using your own shapefile and single band raster:

enter image description here

Which results in the mean pixel value for each polygon.

enter image description here

require(raster)
require(maptools)

# Read the polygon shapefile
poly = readShapePoly("C:/temp/poly.shp")
plot(poly)

# Read the single band raster
raster = raster("C:/temp/subset.tif")

# Extract the raster values underlying the polygons
v <- extract(raster, poly, fun = mean)
output = data.frame(v)
print(output)

[GIS] Step-by-step: How to extract Raster values from Polygon overlay with Q-GIS or R

One problem I see is that if you do not have all classes present in each polygon then the resulting data does not match the dimension of the data frame in the polygon @data slot. You can avoid a looping function entirely. Here is a quick fix that will add NA values where you do not have a class present.

require(raster)

# Create interger class raster
r <- raster(ncol=36, nrow=18)
  r[] <- round(runif(ncell(r),1,10),digits=0)

# Create two polygons
cds1 <- rbind(c(-180,-20), c(-160,5), c(-60, 0), c(-160,-60), c(-180,-20))
cds2 <- rbind(c(80,0), c(100,60), c(120,0), c(120,-55), c(80,0))
polys <- SpatialPolygonsDataFrame(SpatialPolygons(list(Polygons(list(Polygon(cds1)), 1), 
                                  Polygons(list(Polygon(cds2)), 2))),data.frame(ID=c(1,2)))

# Extract raster values to polygons                             
( v <- extract(r, polys) )

# Get class counts for each polygon
v.counts <- lapply(v,table)

# Calculate class percentages for each polygon
( v.pct <- lapply(v.counts, FUN=function(x){ x / sum(x) } ) )

# Create a data.frame where missing classes are NA
class.df <- as.data.frame(t(sapply(v.pct,'[',1:length(unique(r)))))  

# Replace NA's with 0 and add names
class.df[is.na(class.df)] <- 0   
  names(class.df) <- paste("class", names(class.df),sep="")

# Add back to polygon data
polys@data <- data.frame(polys@data, class.df)
  head(polys@data)

Best Answer

Related Solutions

[GIS] combine raster and polygon values in R

[GIS] Step-by-step: How to extract Raster values from Polygon overlay with Q-GIS or R

Related Question