loading packages
library(rgdal)
library(sp)
library(raster)
create raster
df <- data.frame( x = rep( 0:4, each=5 ),y = rep( 0:4, 5),l = 1:25)
coordinates(df) <- ~ x + y # coerce to spatial points data frame by assigning x and y values as coordinates
gridded(df) <- TRUE # coerce to SpatialPixelsDataFrame
r <- raster(df) # coerce to raster
create spatial polygon
coords = matrix(c(2,1.5, 2,3.5, 3.5,3.5, 3.5,1.5, 2,1.5 ), ncol = 2, byrow = TRUE) # coordinates
p1 = Polygon(coords)
pl = SpatialPolygons(list(Polygons(list(p1), ID = "a")))
plotting
plot(r)
text(r)
plot(pl, add = TRUE)
Simple operations:
taking spatial average of raster values inside polygon
r.l = extract(r,pl,weights = T) # taking area weights
r.l
## [[1]]
## value weight
## [1,] 14 0.1666667
## [2,] 19 0.3333333
## [3,] 13 0.1666667
## [4,] 18 0.3333333
v = sum(r.l[[1]][,1] * r.l[[1]][,2])
v
## [1] 16.83333
r.l2 = extract(r,pl)
r.l2
## [[1]]
## [1] 14 19 13 18
lapply(r.l2,FUN = mean) # without taking area weights
## [[1]]
## [1] 16
Comments
Post a Comment