zonal.stats: zonal.stats

Description Usage Arguments Value Note Author(s) Examples

View source: R/zonal.stats.R

Description

Polygon zonal statistics of a raster

Usage

1
zonal.stats(x, y, stat, trace = TRUE, plot = TRUE)

Arguments

x

Polygon object of class SpatialPolygonsDataFrame

y

Raster object of class raster

stat

Statistic or function

trace

Should progress counter be displayed

plot

Should subset polygons/rasters be plotted (TRUE/FALSE)

Value

Vector, length of nrow(x), of function results

Note

This function iterates through a polygon object, masks the raster to each subset polygon and then coerces the subset raster to a vector object. The resulting vector is then passed to the specified statistic/function. This is much slower than zonal functions available in GIS software but has the notable advantage of being able to accept any custom function, passed to the 'stat' argument, appropriate for a vector object (see example).

Depends: sp, raster

Author(s)

Jeffrey S. Evans <[email protected]>

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
                                                                      
# skewness function
skew <- function(x, na.rm=FALSE) { 
   if (na.rm) 
       x <- x[!is.na(x)]
  sum( (x - mean(x)) ^ 3) / ( length(x) * sd(x) ^ 3 )  
  }   

# percent x >= p function
pct <- function(x, p=0.30) {
  if ( length(x[x >= p]) < 1 )  return(0) 
    if ( length(x[x >= p]) == length(x) ) return(1) 
     else return( length(x[x >= p]) / length(x) ) 
}

# create some example data
library(raster)
library(sp)   
p <- raster(nrow=10, ncol=10)
  p[] <- runif(ncell(p)) * 10
    p <- rasterToPolygons(p, fun=function(x){x > 9})
      r <- raster(nrow=100, ncol=100)
        r[] <- runif(ncell(r)) 
plot(r)
  plot(p, add=TRUE, lwd=4) 

# run zonal statistics using skew and pct functions   
z.skew <- zonal.stats(x=p, y=r, stat=skew, trace=TRUE, plot=TRUE) 
z.pct <- zonal.stats(x=p, y=r, stat=pct, trace=TRUE, plot=TRUE)
  ( z <- data.frame(ID=as.numeric(as.character(row.names(p@data))), 
                    SKEW=z.skew, PCT=z.pct) )  

jeffreyevans/spatialEco documentation built on Aug. 11, 2018, 1:08 p.m.