In dakni/TestPackage: Just a test package

knitr::opts_chunk$set(echo = TRUE,
                      cache = FALSE)

aha

Look at Figure \@ref(fig:testing), it shows the output of a magic function.

library(TestPackage)
a <- add_2(seq(1:10))
plot(a)

Adding 2 to 4 we get r add_2(4).

knitr::kable(head(cars),
             caption = "A table of the cars package"
             )

Look at the table \@ref(tab:table).

Thanks to @bookdown, @rmarkdown, @knitr

Xie knows what he (?) is doing [-@knitr].

add_3(5)

Caching Tests

library(tidyverse)
x <- rnorm(1e6) - 5e3
knitr::kable(head(x))

knitr::kable(head(x + 2))

second column

selsecond(mtcars)

Packrat test

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

harran <- read.table("../data/Sites_HarranPlain.csv", sep = ",", header=TRUE)
str(harran)

library(sp)
coordinates(harran) <- ~X+Y
proj4string(harran) <- CRS("+init=epsg:4326")
str(harran)

library(raster)
##srtm <- getData("SRTM", lon=38, lat=37)
srtm <- raster("srtm_44_05.tif")
plot(srtm)
points(harran)
srtm <- crop(srtm, extent(harran)+1)
plot(srtm)
srtm <- projectRaster(srtm, crs = CRS("+init=epsg:32637"))
srtm2 <- aggregate(srtm, fact = 2)
writeRaster(srtm2, "data/dem.tif", overwrite = TRUE)

create point pattern object

harran <- spTransform(harran, CRSobj = CRS("+init=epsg:32637"))
library(spatstat)
harran_ppp <- ppp(x = harran@coords[,1],
              y = harran@coords[,2],
              window = owin(xrange = harran@bbox[1,],
                            yrange = c(min(harran@coords[,2]), min(harran@coords[,2]+52000))))
plot(harran)

challenge: delete duplicated points

anyDuplicated(harran_ppp)
harran <- unique(harran_ppp)
plot(harran_ppp)

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

challenge create kernel density estimation

harran_kde <- density.ppp(x = harran_ppp, sigma = mean(harran_ppp_nn))
plot(harran_kde)

raster

library(raster)
dem <- raster("../data/dem.tif")

im_dem <- as.im(as.image.SpatialGridDataFrame(as(dem, "SpatialGridDataFrame")))

harran_rhohat <- rhohat(object = harran_ppp, 
                        covariate = im_dem,
                        bw = 200
                        )
plot(harran_rhohat)
str(harran_rhohat)
rho_dem <- predict(harran_rhohat)
plot(rho_dem)
diff_rho <- harran_kde - rho_dem

create random points with rpoispp function that have the same intensity like our point pattern.

set.seed(123)
harran_poispp2 <- rpoispp(lambda = harran_ppp$n/area.owin(harran_ppp$window), win = harran_ppp$window)
set.seed(123)
harran_poispp3 <- rpoispp(intensity(harran_ppp), win=Window(harran_ppp))
set.seed(123)
harran_poispp4 <- rpoispp(ex = harran_ppp)

plot(harran_ppp)
points(harran_poispp2, col = "red")
points(harran_poispp3, col = "blue")
points(harran_poispp4, col = "green")

Second order effects

harran_g <- Gest(harran_ppp)
str(harran_g)
plot(harran_g)
harran_ge <- envelope(harran_ppp, fun = "Gest", nsim = 999)
plot(harran_ge)

harran_f <- Fest(harran_ppp)
harran_fe <- envelope(harran_ppp, fun = "Fest", nsim = 99)

harran_k <- Kest(harran_ppp)
harran_ke <- envelope(harran_ppp, fun = "Kest", nsim = 99)

plot(harran_fe)
plot(harran_ke)

Inhomogeneous G/F/K

harran_gi <- Ginhom(harran_ppp, lambda = predict(harran_rhohat))
#par(mfrow = c(1,2))
#plot(harran_gi, xlim = c(0,6000))
#plot(harran_g, xlim = c(0,6000))

Interpolation

load("../data/Precipitation.RData")
test <- data.frame(test)
library(sp)
coordinates(test) <- ~lon+lat
proj4string(test) <- CRS("+init=epsg:4326")
test2 <- spTransform(test, CRS("+init=epsg:32634"))
plot(test2)
library(raster)
##srtm <- getData("SRTM", lon=mean(coordinates(test)[,1]), lat=mean(coordinates(test)[,2]))
srtm <- raster("../data/srtm_41_05.tif")
srtm <- crop(srtm, extent(test)+1)
plot(srtm)
srtm3 <- projectRaster(srtm, crs=CRS("+init=epsg:32634"))
srtm3 <- aggregate(srtm3, fact = 3)
library(gstat)

IDW

rain_idw <- idw(mean_r~1, test2, as(srtm3, "SpatialGridDataFrame"))
rain_idw_cv <- krige.cv(mean_r~1, test2)
head(rain_idw_cv)
plot(rain_idw)

Kriging

plot(variogram(mean_r~1, loc = test2, cloud = TRUE))
va <- variogram(mean_r~1, loc = test2, cutoff = 80000, width = 5000)
plot(va)
plot(va, vgm(8e+04, "Exp", 5e+04, 5000))
fva <- fit.variogram(va, vgm(8e+04, "Exp", 5e+04, 5000),
fit.method = 7) ## ordinary least squares; default
plot(va, fva)
rain_krige <- krige(mean_r~1, test2, as(srtm3, "SpatialGridDataFrame"), fva)
rain_krige_cv <- krige.cv(mean_r~1, test2, fva)
head(rain_krige_cv)
plot(rain_krige)
rain_krige2 <- brick(rain_krige)
plot(rain_krige2)

Automatic Kriging

library(automap)
rainautkrige <- automap::autoKrige(mean_r~1, test2, as(srtm3, "SpatialGridDataFrame"))
str(rainautkrige)
plot(rainautkrige)
tmp <- rainautkrige[[1]]
tmp <- raster(tmp)
tmp2 <- raster(rain_krige)
tmp3 <- tmp2-tmp
plot(tmp3)
## or using Reduce function (that needs a list)
plot(Reduce("-", list(tmp2,tmp)))

Kriging with external drift

names(srtm3) <- "altitude"
test3 <- test2[!is.na(test2$altitude),]
rainkwed <- automap::autoKrige(mean_r~altitude, test3, as(srtm3, "SpatialGridDataFrame"))
plot(rainkwed)

Kriging with external drift (manual)

see Wolfgang's script.

linmod <- lm(mean_r~altitude, test3)
linmod

plot(variogram(linmod$residuals~1, loc = test3, cloud = TRUE))

the easy way: altitude + x

see Wolfgang's script.

test3$x <- coordinates(test3)[,1]

library(lattice); library(gstat)

v <- variogram(mean_r ~ 1, test3)
xyplot(
    x = gamma/1e8 ~ dist,
    data = v,
    pch = 3,
    type = 'b',
    lwd = 2,
    panel = function(x, y, ...) {
        for (i in 1:500) {
            test3$random = sample(test3$mean_r)
            v = variogram(random ~ 1, test3)
            llines(x = v$dist,
                   y = v$gamma/1e8,
                   col = rgb(.5,.5,.5,.2)
                   )
        }
        panel.xyplot(x, y, ...)
    },
    xlab = 'distance',
    ylab = 'semivariance/1e8'
)

References

dakni/TestPackage documentation built on May 30, 2019, 8:02 a.m.