Category Archives: R

VpSpPh2013WS-W11-1-CS-tappelhans

Code snippet: extraction of points in polygons for point pattern analysis


### shape data
library(rgdal)
library(maptools)
library(raster)
library(spatstat)

setwd("/media/windows/tappelhans/uni/marburg/lehre/2013/ws/spezielle")

lyr <- ogrListLayers("lidar_data/trees.shp")
trees <- readOGR("lidar_data/trees.shp", layer = lyr)
rst <- raster("lidar_data/rst/grid_raster1.tif")

grid <- rasterToPolygons(rst, na.rm = FALSE)
grid.pol <- as(grid, "SpatialPolygons")

trees.polID <- over(x = trees[, 3], y = grid.pol)
trees$polID <- trees.polID

trees.polID.compl <- trees.polID[!is.na(trees.polID)]
trees.compl <- trees[!is.na(trees.polID), ]

dat.ppp <- lapply(unique(trees.polID.compl), function(i) {
  if (any(unique(trees.polID.compl == i))) {
    tmp <- trees.compl[trees.polID.compl == i, ]
    print(i)
    p <- as(tmp[, 4], "ppp")
    w <- as(grid.pol[i], "owin") 
    p$window <- w
  }
  return(p)
})

dat.ppp <- dat.ppp[!sapply(dat.ppp, is.null)] 

### visual interpretation of point pattern (exploratory)
plot(dat.ppp[[67]])
plot(density(dat.ppp[[100]]))
plot(Kest(dat.ppp[[9]]))
plot(Lest(dat.ppp[[100]]))

###

kest.ls <- lapply(seq(dat.ppp), function(i) {
  if (dat.ppp[[i]]$n > 10) {
    print(i)
    tmp <- Kest(dat.ppp[[i]])
    mn <- mean(tmp$iso - tmp$theo) / max(tmp$theo)
    return(mn)
  }  
})

kest.ls[sapply(kest.ls, is.null)] <- NA

df.pol.kest <- do.call("rbind", lapply(seq(dat.ppp), function(i) {
  data.frame(pol = unique(dat.ppp[[i]]$marks),
             mkest = kest.ls[[i]])
}))

pols <- unlist(lapply(seq(nrow(grid)), function(i) {
  as.numeric(grid@polygons[[i]]@ID)
}))

rst2 <- rst 
rst2.vals <- unique(merge(data.frame(pol = pols), 
                          df.pol.kest, all.x = TRUE))
rcl <- as.matrix(rst2.vals)

rst2[] <- pols
rst2 <- reclassify(rst2, rcl)

plot(rst2)