WS 2016/2017
Prerequisites
setwd, library, data, cleaning
setwd("YOUR-PATH")
raw_data <- read.table(file = "data/FS_UTM.csv",
header = TRUE,
sep = ";",
stringsAsFactors = FALSE,
quote = "")
## check for issues/problems
## str(raw_data)
## apply(X=raw_data[,4:20], MARGIN = 2, FUN=function(x){x[x>1]})
raw_data$settlement[raw_data$settlement > 1] <- 1
library(sp)
sites <- raw_data[(raw_data$settlement == 1 & raw_data$iron.age == 1),]
coordinates(sites) <- ~xUTM+yUTM
proj4string(sites) <- CRS("+init=epsg:32634")
Create Point Pattern dataset
library(spatstat)
sites_pp <- ppp(x = sites@coords[,1],
y = sites@coords[,2],
window = owin(xrange = c(min(sites@coords[,1]),
max(sites@coords[,1])
),
yrange = c(min(sites@coords[,2]),
max(sites@coords[,2])
)
)
)
unitname(sites_pp) <- c("meter", "meters")
anyDuplicated(sites_pp)
## [1] 0
sites_pp <- sites_pp[duplicated(sites_pp)==FALSE] anyDuplicated(sites_pp)
## [1] 0
#plot(sites_pp)
Caluclate density using graph methods
Distance based density using "largest empty circle"
Idea: calculate the centres of the largest empty circles – corresponds to the edges of the Voronoi graph
library(tripack)
sites_vor <- voronoi.mosaic(x = sites_pp$x,
y = sites_pp$y,
duplicate = "remove"
)
rad <- sites_vor$radius
lec <- SpatialPointsDataFrame(coords = cbind(sites_vor$x, sites_vor$y),
data = as.data.frame(rad), # length of edge
proj4string=CRS("+init=epsg:32634")
)
Caluclate density using graph methods
plot(lec)
## remove outlier lec <- lec[lec@coords[,2]>min(lec@coords[,2]),] plot(lec)
Caluclate density using graph methods
plot(sites_pp) plot(sites_vor, col = "red", add=TRUE)
Caluclate density using graph methods
plot(sites_pp, pch = 19)
points(lec, col="red")
library(plotrix)
for (i in 1:length(lec@coords[,1])) {
draw.circle(x = lec@coords[i,1],
y = lec@coords[i,2],
radius = lec@data$rad[i],
border = "gray")
}
Interpolation
We're lazy…instead of creating a raster we use one that is "already there":
library(raster)
srtm <- raster("results/srtm.tif")
srtm <- aggregate(x = srtm, fact = 9)
library(maptools)
srtm <- as(object = srtm,
Class = "SpatialGridDataFrame"
)
str(srtm)
## Formal class 'SpatialGridDataFrame' [package "sp"] with 4 slots ## ..@ data :'data.frame': 64220 obs. of 1 variable: ## .. ..$ srtm: num [1:64220] 74.7 75 76.6 78.7 75.1 ... ## ..@ grid :Formal class 'GridTopology' [package "sp"] with 3 slots ## .. .. ..@ cellcentre.offset: Named num [1:2] 433526 5001984 ## .. .. .. ..- attr(*, "names")= chr [1:2] "s1" "s2" ## .. .. ..@ cellsize : num [1:2] 810 810 ## .. .. ..@ cells.dim : int [1:2] 338 190 ## ..@ bbox : num [1:2, 1:2] 433121 5001579 706901 5155479 ## .. ..- attr(*, "dimnames")=List of 2 ## .. .. ..$ : chr [1:2] "s1" "s2" ## .. .. ..$ : chr [1:2] "min" "max" ## ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot ## .. .. ..@ projargs: chr "+proj=utm +zone=34 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"
##image(srtm) ##points(sites, pch=19, cex=.5)
Interpolation - Inverse Distance Weighting
library(gstat)
crs(srtm) <- lec@proj4string ## BE CAREFUL!
lec_gstat_p2 <- gstat(formula = lec@data$rad ~ 1,
locations = lec,
set = list(idp = 2)
)
lec_gstat_p25 <- gstat(formula = lec@data$rad ~ 1,
locations = lec,
set = list(idp = 2.5)
)
lec_idw2 <- predict(object = lec_gstat_p2,
newdata = srtm)
## [inverse distance weighted interpolation]
head(as.data.frame(lec_idw2))
## var1.pred var1.var s1 s2 ## 1 18090.80 NA 433526.3 5155074 ## 2 18124.37 NA 434336.3 5155074 ## 3 18158.62 NA 435146.3 5155074 ## 4 18193.53 NA 435956.3 5155074 ## 5 18229.13 NA 436766.3 5155074 ## 6 18265.40 NA 437576.3 5155074
lec_idw25 <- predict(object = lec_gstat_p25,
newdata = srtm)
## [inverse distance weighted interpolation]
Interpolation - Inverse Distance Weighting
Interpolation - Trend surfaces
lecp2 <- lm(
lec@data$rad ~ poly(x = lec@coords[,1],
y = lec@coords[,2],
degree = 2
)
)
lecp3 <- lm(
lec@data$rad ~ poly(x = lec@coords[,1],
y = lec@coords[,2],
degree = 3
)
)
plot(lec@data$rad ~ 1)
lines(lecp2$fitted.values, col = "red")
lines(lecp3$fitted.values, col = "blue")
legend("top",
lty = c(1,1,1),
col = c("red","blue","green"),
legend = c("2nd degree polynom",
"3rd degree polynom"
)
)
Interpolation - Trend surfaces
lec_poly <- krige(
formula = lec@data$rad ~ 1,
locations = lec,
newdata = srtm,
degree = 3
)
## [ordinary or weighted least squares prediction]
plot(raster(lec_poly))
Interpolation
Correlation between point pairs at different separation distances.
library(lattice)
hscat(formula = lec@data$rad ~ 1,
data = lec,
breaks = seq(0, 10, 1)*1000
)
Interpolation – Variogram
library(gstat)
plot(variogram(object = lec@data$rad ~ 1,
locations = lec,
cloud = TRUE
)
)
Interpolation - Kriging
library(gstat) plot(variogram(lec$rad ~ 1, lec, cloud = TRUE))
library(gstat) plot(variogram(lec$rad ~ 1, lec))
Interpolation – Variogram
library(lattice); library(gstat)
v <- variogram(lec@data$rad ~ 1, lec)
xyplot(
x = gamma/1e8 ~ dist,
data = v,
pch = 3,
type = 'b',
lwd = 2,
panel = function(x, y, ...) {
for (i in 1:500) {
lec$random = sample(lec$rad)
v = variogram(random ~ 1, lec)
llines(x = v$dist,
y = v$gamma/1e8,
col = rgb(.5,.5,.5,.2)
)
}
panel.xyplot(x, y, ...)
},
xlab = 'distance',
ylab = 'semivariance/1e8'
)
Interpolation
show.vgms()
Interpolation
plot(variogram(
object = lec@data$rad ~ 1,
locations = lec,
boundaries = c(seq(0, 1e+3, 1e+2),
seq(2e+3, 1e+4, 1e+3)
),
cutoff = 40000
)
)
Interpolation
plot(variogram(object = lec@data$rad ~ 1,
locations = lec,
cloud = FALSE,
width = 1000,
cutoff = 40000
)
)
Interpolation - Kriging
library(gstat) plot(variogram(lec$rad ~ 1, lec, cutoff = 40000, alpha = c(0,45,90,135)))
Interpolation - Kriging
When fitting goes wrong…
vt <- variogram(object = lec@data$rad ~ 1,
locations = lec,
width = 5000,
cutoff = 40000
)
vt_vgm <- vgm(psill = 8.0e+07,
model = "Gau",
range = 40000,
nugget = .2e+07)
v_fit <- fit.variogram(vt,vt_vgm)
## Warning in fit.variogram(vt, vt_vgm): No convergence after 200 iterations: ## try different initial values?
plot(vt,v_fit)
Interpolation – tool: eyefit
library(geoR) v_eye <- eyefit(variog(as.geodata(lec["rad"]), max.dist = 40000)) ve_fit <- as.vgm.variomodel(v_eye[[1]]) v_fit <- fit.variogram(vt,ve_fit) plot(vt,v_fit)
Interpolation - Kriging
As it gets obvious: the chosen method (and/or) its application leads to completely useless results.
lec_kri <- krige(formula = lec$rad ~ 1,
location = lec,
newdata = srtm,
model = v_fit
)
## [using ordinary kriging]
plot(raster(lec_kri),
main="Kriging using Gaussian fit"
)
points(sites,pch=20,cex=.4)
Inverse Distance Weighting vs. Kriging
Interpolation
So, did any of the things calculated so far make any sense?
Let's assume there is some value in it, what does it mean (landscape) archaeologically?
Why is the analysis at the current stage problematic?
What might be the reasons for the failure?