P13 Point Pattern.R

From gnewarchaeology wiki
Jump to: navigation, search


  • point_pattern.R

/!\ TODO: Syntax highlight + Literate programming with Org-mode

### p13 point pattern analysis ###
### workspace associato: point_pattern.RData ###

## load library 
library(sp)
library(lattice)
library(rgeos)
library(maptools)
library(spatial)
library(spatstat)
library(splancs)
library(spatialkernel)
library(scatterplot3d)
library(ks)
library(misc3d)

# SPATIAL OBJECTS

getClass("Spatial")
getClass("CRS")
getClass("SpatialPoints")

## read data.frame
points <- read.csv("/absolute/path/to/points_all.csv")
records <- read.csv("/absolute/path/to/records_all_var.csv")

## s_index (indice di sfericità=l_max/l_med)
index <- records$l_max/records$l_med
records_index <- cbind(records,index)
class(records_index)

## matrix xyz
xyz <- cbind(points$x, points$y, points$z)
row.names(xyz) <- 1:nrow(xyz)
str(xyz)

## spatial points
#p13CRS <- CRS("+proj=longlat +ellps=WGS84")
#p13crs <- CRS(as.character(NA))
xyz_sp <- SpatialPoints(xyz)
summary(xyz_sp)
bbox(xyz_sp)
A09B5S01 <- which(points$codice == "A09B5S01")
A09B5S01
coordinates(xyz_sp)[A09B5S01, ]
summary(xyz_sp[A09B5S01, ])
print(xyz_sp)

## what's class?
class(xyz_sp)
class(xyz)
class(points)
class(records)
names(points)

## spatial points dataframe
getClass("SpatialPointsDataFrame")
xyz_spdt <- SpatialPointsDataFrame(xyz, points, match.ID =TRUE)
xyz_spdt
str(xyz_spdt)
names(xyz_spdt)
xyz_spdt1 <- SpatialPointsDataFrame(xyz_sp, points, match.ID =TRUE)
all.equal(xyz_spdt, xyz_spdt1)

## spatial points dataframe + records!
## method 1?
records_sp1 <- SpatialPointsDataFrame(xyz_sp, records_index, match.ID =TRUE)
## method 2?
records_sp2 <- records_index
coordinates(records_sp2) <- xyz
all.equal(records_sp1, records_sp2)
str(records_sp1)
## coords.nrs?
## method 3 spCbind !!!
o <- match(xyz_spdt$codice, records_index$codice)
records1 <- records_index[o,]
row.names(records1) <- xyz_spdt$codice
records_sp3 <- spCbind(xyz_spdt, records1)
names(records_sp3)
identical(records$codice, records_sp3$codice)
str(records_sp3)
summary(records_sp3)

# RECORDS_SP3 and subsets
bone_sp <- records_sp3[which(records_sp3$descrizione == "OSSO"), ]
#bone_index_sp <- records_ppp[which(records_ppp$descrizione == "OSSO" && records_ppp$index >= 1.5), ]
summary(bone_sp)
stone_sp <- records_sp3[which(records_sp3$descrizione == "PIETRA"), ]
silex_sp <- records_sp3[which(records_sp3$descrizione == "SELCE"), ]
  
## convert between ppp objects and SpatialPoints
records_ppp <- as(records_sp3, "ppp")
class(records_ppp)
summary(records_ppp)
print(records_ppp)
## alternativa:
records_ppp1 <- as.ppp(records_sp3)
summary(records_ppp1)
print(records_ppp)
bone_ppp <- as.ppp(bone_sp)
stone_ppp <- as.ppp(stone_sp)
silex_ppp <- as.ppp(silex_sp)

## pp3 3D object!!!
records_pp3 <- pp3(records_sp3$x, records_sp3$y, records_sp3$z, box3(c(3,9),c(1,4),c(-2,0)))
summary(records_pp3)
print(records_pp3)
plot(records_pp3)
bone3 <- records_sp3[which(records_sp3$descrizione == "OSSO"), ]
bone_pp3 <- pp3(bone3$x, bone3$y, bone3$z, box3(c(3,9),c(1,4),c(-2,0)))
stone3 <- records_sp3[which(records_sp3$descrizione == "PIETRA"), ]
stone_pp3 <- pp3(stone3$x, stone3$y, stone3$z, box3(c(3,9),c(1,4),c(-2,0)))
silex3 <- records_sp3[which(records_sp3$descrizione == "SELCE"), ]
silex_pp3 <- pp3(silex3$x, silex3$y, silex3$z, box3(c(3,9),c(1,4),c(-2,0)))

# SIMPLE PLOT SYSTEM

spplot(records_sp2, "descrizione", do.log=TRUE)
plot(records_sp3)
plot(records_sp3, pch=16, cex=.5, axes=T)
plot(unmark(records_ppp))
identify.ppp(unmark(records_ppp))

# POINT PATTERN (DENSITY) ANALYSIS
# preliminary analysis of a point pattern

## G function: distance to the nearest event
g <- G3est(records_pp3)
if(interactive()) plot(g)

plot(records_pp3, main="records")
plot(G3est(records_pp3), main="G function")
plot(bone_pp3)
plot(Gbone <- G3est(bone_pp3), main="G function: bones")
plot(stone_pp3)
plot(Gstone <- G3est(stone_pp3), main="G function: stones")
plot(silex_pp3)
plot(Gsilex <- G3est(silex_pp3), main="G function: silex")

## F function: distance from a point to the nearest event
x11()
layout(matrix(c(1,2),1,2))
plot(records_pp3)
plot(F3est(records_pp3), main="F function")

plot(F3est(bone_pp3), main="F function: bones")
plot(F3est(stone_pp3), main="F function: stones")
plot(F3est(silex_pp3), main="F function: silex")

# statistical analysis of spatial point processes
# first order properties11111
# second order properties

## K function (homogeneous point processes)
#par(mfrow =c(1,2))
layout(matrix(c(1,2),1,2))
plot(records_pp3, main ="records")
plot(K3est(records_pp3), main ="K function")

plot(bone_pp3)
plot(K3est(bone_pp3))
plot(stone_pp3)
plot(K3est(stone_pp3))
plot(silex_pp3)
plot(K3est(silex_pp3))

## summary statistics (K, F, G, J)
## estimate of the K function
plot(Kest(records_ppp))
## k3est
plot(K3est(records_pp3))

## 4 main summary functions (K, F, G, J)
#plot(allstats(records_ppp))

x11()
layout(matrix(c(1,2),1,2))
plot(K3est(bone_pp3), main="K function: bones")
plot(F3est(bone_pp3), main="F function: bones")


## Compute a kernel smoothed intensity function from a point pattern
density.ppp(records_ppp)
plot(density.ppp(records_ppp))
plot(density.ppp(bone_ppp))
plot(density.ppp(stone_ppp))
plot(density.ppp(silex_ppp))
persp(density.ppp(silex_ppp))

plot(density(bone_ppp, sigma=bw.diggle(bone_ppp)), main="kernel density estimation")
persp(density.ppp(bone_ppp), main="bones", xlab="x", ylab="y", zlab="z")
plot(density(stone_ppp, sigma=bw.diggle(stone_ppp)), main="kernel density estimation")
persp(density.ppp(stone_ppp), main="stones", xlab="x", ylab="y", zlab="z")
plot(density(silex_ppp, sigma=bw.diggle(silex_ppp)), main="kernel density estimation")
persp(density.ppp(silex_ppp), main="silex", xlab="x", ylab="y", zlab="z")

## kde3d
#contour3d(kde3d(records_sp3$x, records_sp3$y, records_sp3$z))#,# lims = c(6,4,3))

# TODO...

# CLUSTER ANALYSIS
# SPATIAL AUTOCORRELATION
# GEOSTATISTICS