In rsbraun/mytestpkg: demonstrates some concepts

knitr::opts_chunk$set(echo = TRUE)

library(mytestpkg)

Use new functions:

a <- c(5,3,9,5)
b <- c(4,6,7,3)
new_fun(a,b)

Create Table \@ref(fig:table) [after @schauer1968wirkungsspezifitat]:

tab <- matrix(c(3,5,8,3,5,3,6,4), ncol=2)
colnames(tab) <- c("a","b")
rownames(tab) <- c("test1", "test2", "test3", "test4")
knitr::kable(head(tab), 
             caption = "My new table")

Create Plot \@ref(fig:tableplot) :

plot(tab)

For help check @Hijmans2016

library(binford)
data(LRB)
knitr::kable(head(LRB))

Point Pattern

##------------------------
## First Order effects
##------------------------

harran <- read.table("../data/ReReLA/data/Sites_HarranPlain.csv", sep=",", header=TRUE)
head(harran)
## Create Spatial Object 
##--------------------------
library(sp)
coordinates(harran) <- ~X+Y
proj4string(harran) <- CRS("+init=epsg:4326")
harran <- spTransform(harran, CRSobj = CRS("+init=epsg:32637"))

## Create Point pattern object
##------------------------------
library(spatstat)
harran_ppp <- ppp(x = harran@coords[,1],
                  y = harran@coords[,2],
                  window = owin(xrange = harran@bbox[1,],
                                yrange = c(min(harran@bbox[2,]), min(harran@coords[,2]+52000)))
                  )
# remove duplicated points
str(harran_ppp)
harran_ppp2 <- harran_ppp[!duplicated(harran_ppp)]
str(harran_ppp2)
plot(harran_ppp2)

# or:
anyDuplicated(harran_ppp)
harran_ppp3 <- unique(harran_ppp)

harran_ppp_nn <- nndist(harran_ppp2)
str(harran_ppp_nn)
hist(harran_ppp_nn)

## Kernel Density Estimation
##-----------------------------
kde <- density.ppp(x=harran_ppp2, sigma = mean(harran_ppp_nn))   
plot(kde)

# use another bandwidth (sigma)
# for clustered data use diggle
bw.ppl(harran_ppp2)
bw.diggle(harran_ppp2)
plot(bw.ppl(harran_ppp2))
plot(bw.ppl(harran_ppp2), xlim=c(2000,5000))

## Add a Covariate
##---------------------------
# load raster
library(raster)
dem <- raster("../data/ReReLA/data/dem_harran.tif")

# convert raster to pixel image
t <- as(dem, "SpatialGridDataFrame")
im.dem <- as.im(as.image.SpatialGridDataFrame(as(dem, "SpatialGridDataFrame")))
harran_rhohat <-rhohat(object = harran_ppp2, covariate = im.dem, bw = mean(harran_ppp_nn))     # window of points matters
plot(harran_rhohat)   # y-axis: intensity of points
# check bandwidth (sigma):
str(harran_rhohat)
# Predict
rho_dem <- predict(harran_rhohat)
plot(rho_dem)
# Compare raster with real data with predicted raster
diff_rho <- kde - rho_dem
plot(raster(diff_rho))

## Test against random poisson process
## -------------------------------------
# create random points with the same density like the real points
# compute density - Points per area
dens1 <- harran_ppp2$n/area.owin(harran_ppp2$window)
#set.seed(123)
harran_poispp1 <- rpoispp(lambda = dens1, win = harran_ppp2$window)   # poisson - complete spatial randomness 

## or /error
dens2 <- intensity(harran_ppp2)
#set.seed(123)
harran_poispp2 <- rpoispp(lambda = dens2, win = harran_ppp2$window)

## or:
#set.seed(123)
harran_poispp3 <- (ex = harran_ppp)

plot(harran_ppp2)
points(harran_poispp1, col="red")
points(harran_poispp2, col="blue")
points(harran_poispp3, col="green")

##-------------------------
## Second order effects
##------------------------
# G-Function
harran_g <- Gest(harran_ppp2)
str(harran_g)
plot(harran_g)
 # generate 99 random points and run G-Function 
harran_ge <- envelope(harran_ppp2, fun="Gest")
plot(harran_ge)
# greyish curve - result of all 99 cases - becomes wider with more simulations
harran_ge1000 <- envelope(harran_ppp2, fun="Gest", nsim=1000)
plot(harran_ge1000)

# F-Function
harran_f <- Fest(harran_ppp2)
str(harran_f)
plot(harran_f)

harran_fe1000 <- envelope(harran_ppp2, fun="Fest", nsim=1000)
plot(harran_fe1000)

# K-Function
harran_k <- Kest(harran_ppp2)
str(harran_k)
plot(harran_k)

harran_ke1000 <- envelope(harran_ppp2, fun="Kest", nsim=1000)
plot(harran_ke1000)

# Inhomogeneous G/F/K:
harran_gi <- Ginhom(harran_ppp2, lambda=predict(harran_rhohat))
par(mfrow=c(1,2))
plot(harran_gi, xlim=c(0,500))
plot(harran_g, xlim=c(0,500))


#Finhom()
#Kinhom()