VpSpPhWS13-L04-GISnost

This script adds the z-values to the existing x- and y-values. For this the measured distance (dis) and the radiation inclination (incrad) is needed. After that a raster DEM ist calculated.

library(latticeExtra)
library(spatstat)
library(sp)
library(automap)
library(raster)

# Set path
path <- "H:/.ntdesktop/VpSpPh2013WS-L04"
setwd(path)

# Read data
data <- read.csv("complete.data.csv")

# Set value of P1
P1x <- 481073
P1y <- 5645524
P1z <- 296

# Calculate values of P2, P3 and P4
P4x <- P1x + ((sin(5 * (pi/180))) * 40)
P4y <- P1y + ((cos(5 * (pi/180))) * 40)
P3x <- P4x + ((sin(275 * (pi/180))) * 40)
P3y <- P4y + ((cos(275 * (pi/180))) * 40)
P2x <- P3x + ((sin(185 * (pi/180))) * 40)
P2y <- P3y + ((cos(185 * (pi/180))) * 40)

# Set z value
P2z <- 295.5
P3z <- 309.8
P4z <- 309.9



# Calculate inclination degree to radiation
incrad <- function(inc) {
    incradiation <- inc * (pi/180)
    return(incradiation)
}

# Set column for inclination radiants
incrads <- incrad(data$inc)
data$incrads <- incrads


# Calculate 2D-beeline
bl <- function(incrad, dis) {
    beeline <- cos(incrad) * dis
    return(beeline)
}

# Set column for 2D-beeline
bl <- bl(data$incrad, data$dis)
data$bl <- bl


# Convert cbh from centimetres to metres
convcbh <- function(cbh) {
    cbhm <- cbh/100
    return(cbhm)
}

# Set column for converted cbh, now cbhm
cbhm <- convcbh(data$cbh)
data$cbhm <- cbhm


# Calculate full distance (2D-beeline + CBH)
fulldist <- function(bl, cbhm) {
    fulldistance <- bl + (cbhm/(2 * pi))
    return(fulldistance)
}

# Set column for full distance
fulldist <- fulldist(data$bl, data$cbhm)
data$fulldist <- fulldist


# Convert degree to radiation
rad <- function(dir) {
    radiation <- dir * (pi/180)
    return(radiation)
}

# Set column for radiation
rads <- rad(data$dir)
data$rads <- rads


# Calculate x, y and z values of the start point trees
B13x <- P3x + sin(3.1764992 * (pi/180)) * 12.520145
B14x <- P3x + sin(2.8972466 * (pi/180)) * 10.011164
B21x <- P3x + sin(1.8325957 * (pi/180)) * 16.461094
B18x <- B21x + sin(3.9095375 * (pi/180)) * 8.692974

B13y <- P3y + cos(3.1764992 * (pi/180)) * 12.520145
B14y <- P3y + cos(2.8972466 * (pi/180)) * 10.011164
B21y <- P3y + cos(1.8325957 * (pi/180)) * 16.461094
B18y <- B21y + cos(3.9095375 * (pi/180)) * 8.692974

B13z <- P3z + sin(-6.5 * (pi/180)) * 11.7
B14z <- P3z + sin(-2 * (pi/180)) * 9.42
B21z <- P3z + sin(1 * (pi/180)) * 13.73
B18z <- B21z + sin(-12 * (pi/180)) * 8.02

# Calculate final value x
data$finalvaluex <- ifelse(data$origin == "P1", P1x + (sin(data$rads) * data$fulldist), 
    ifelse(data$origin == "P2", P2x + (sin(data$rads) * data$fulldist), ifelse(data$origin == 
        "P3", P3x + (sin(data$rads) * data$fulldist), ifelse(data$origin == 
        "P4", P4x + (sin(data$rads) * data$fulldist), ifelse(data$origin == 
        "B13", B13x + (sin(data$rads) * data$fulldist), ifelse(data$origin == 
        "B14", B14x + (sin(data$rads) * data$fulldist), ifelse(data$origin == 
        "B21", B21x + (sin(data$rads) * data$fulldist), ifelse(data$origin == 
        "B18", B18x + (sin(data$rads) * data$fulldist), NA))))))))

# Calculate final value y
data$finalvaluey <- ifelse(data$origin == "P1", P1y + (cos(data$rads) * data$fulldist), 
    ifelse(data$origin == "P2", P2y + (cos(data$rads) * data$fulldist), ifelse(data$origin == 
        "P3", P3y + (cos(data$rads) * data$fulldist), ifelse(data$origin == 
        "P4", P4y + (cos(data$rads) * data$fulldist), ifelse(data$origin == 
        "B13", B13y + (cos(data$rads) * data$fulldist), ifelse(data$origin == 
        "B14", B14y + (cos(data$rads) * data$fulldist), ifelse(data$origin == 
        "B21", B21y + (cos(data$rads) * data$fulldist), ifelse(data$origin == 
        "B18", B18y + (cos(data$rads) * data$fulldist), NA))))))))

# Calculate final value z
data$finalvaluez <- ifelse(data$origin == "P1", P1z + (sin(data$incrads) * data$dis), 
    ifelse(data$origin == "P2", P2z + (sin(data$incrads) * data$dis), ifelse(data$origin == 
        "P3", P3z + (sin(data$incrads) * data$dis), ifelse(data$origin == "P4", 
        P4z + (sin(data$incrads) * data$dis), ifelse(data$origin == "B13", B13z + 
            (sin(data$incrads) * data$dis), ifelse(data$origin == "B14", B14z + 
            (sin(data$incrads) * data$dis), ifelse(data$origin == "B21", B21z + 
            (sin(data$incrads) * data$dis), ifelse(data$origin == "B18", B18z + 
            (sin(data$incrads) * data$dis), NA))))))))


### create spatial data object
dat.sp <- data[complete.cases(data), ]
coordinates(dat.sp) <- c("finalvaluex", "finalvaluey")
proj4string(dat.sp) <- "+proj=utm +ellps=WGS84 +zone=32 +units=m +north"

###### CREATE DEM kriging interpolation of z values
bbx <- dat.sp@bbox
ext <- extent(bbx)
grd <- raster(ext, 40, 40)

grd <- setValues(grd, rnorm(ncell(grd)))
grdsp <- as(grd, "SpatialPixelsDataFrame")

dem <- autoKrige(data$finalvaluez ~ 1, dat.sp, grdsp)
## [using ordinary kriging]


clrs <- colorRampPalette(rev(brewer.pal(9, "YlGnBu")))  #Spectral

dat.xyz <- data.frame(x = dem$krige_output@coords[, 1], y = dem$krige_output@coords[, 
    2], z = dem$krige_output@data$var1.pred)

zx.ratio <- diff(range(dat.xyz$z))/diff(range(dat.xyz$x))

wireframe(z ~ x + y, data = dat.xyz, col.regions = clrs(1000), asp = c(1, zx.ratio), 
    par.box = list(col = NA), ylab = "", xlab = "", zlab = "", scales = list(draw = FALSE), 
    shade = F, screen = list(z = 0, x = -65, y = 30), drape = T, col = "black", 
    zoom = 1.2, colorkey = FALSE)

plot of chunk unnamed-chunk-1

Leave a Reply