Monthly Archives: November 2013

VpSpPhWS13-L04-Geo2013

This is our CODE and it does the following:

After calculating all important values (rads, final_distance.) we set plotorigins P1-P4 in order to calculate the final treepositions.
Therefore we define the x and y-coordinates for the further origins which were used for the treemeasurement at P3.
Furthermore we use an ifelse-statement and plot the results.

Homework L04:

Calculating z-coordinates and visualizing with a 3D-plot

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

dat <- read.csv("H:/.ntdesktop/Felddaten_csv/complete.data.csv")

### set origin

P1X <- 481073
P1Y <- 5645524
zP1 <- 296

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)

### calculating z p1 to a44
zp1a44 = (sin(0.03490659)) * 26.53

# a44 to p2
za44p2 = (sin(0.1047197)) * 13.6
#### HEIGHT P2
zP2 = zP1 + (zp1a44 - za44p2)
# p1 to a31
zp31 = (sin(0.33161256)) * 26.35

# a31 to p4

z31p4 = (sin(-0.2268928)) * 23
#### HEIGHT P4
zP4 = zP1 + (zp31 - z31p4)

# p4 to b22

zp4b22 = 0

# b22 to b21

zb22b21 = (sin(-0.06981317)) * 5.03

# b21 to p3

zb21p3 = (sin(0.01745329)) * 13.37
#### HEIGHT P3
zP3 = zP4 + ((zp4b22) - (zb22b21) - (zb21p3))
### set column for radiants and inclination

### zb13,b14,b18,b21

zB13P3 = sin(-0.1134464) * 11.7
zB13 = zP3 + zB13P3

zB14P3 = sin(-0.03490659) * 9.42
zB14 = zP3 + zB14P3

zB18P3 = sin(-0.20943951) * 8.02
zB18 = zP3 + zB18P3

zB21P3 = sin(0.01745329) * 13.73
zB21 = zP3 - (zB21P3)

rad <- function(dir) {
    radiants <- (dir + 1.68) * (pi/180)  # radianten=direction*(pi/180)
    return(radiants)
}
rads <- rad(dat$dir)
dat$rads <- rads

incli <- function(inc) {
    incl <- inc * (pi/180)  # radinc=inclination*(pi/180)
    return(incl)
}
radinc <- incli(dat$inc)
dat$radinc <- radinc

### set column for cbhm

cbhm <- function(cbh) {
    cbhmeter <- cbh/100
    return(cbhmeter)
}

cbhmet <- cbhm(dat$cbh)
dat$cbhmet <- cbhmet

### set column for final_distance

fdist <- function(cbhmet, dis, radinc) {
    fudi <- (cbhmet/(2 * pi)) + (cos(radinc) * dis)  # Finale Distanz=Baumradius+(cos(radinc)*distanz)
    return(fudi)
}

final_distance <- fdist(dat$cbhmet, dat$dis, dat$radinc)
dat$final_distance <- final_distance

### calculate further origins
XB13 <- P3X + sin(3.1764992) * (12.520145)
XB14 <- P3X + sin(2.8972466) * (10.01)
XB18 <- P3X + sin(3.9095375) * (8.692)
XB21 <- P3X + sin(1.8325957) * (16.46)
YB13 <- P3Y + cos(3.1764992) * (12.520145)
YB14 <- P3Y + cos(2.8972466) * (10.01)
YB18 <- P3Y + cos(3.9095375) * (8.692)
YB21 <- P3Y + cos(1.8325957) * (16.46)

### calculate x-coordinates

dat$xfinal <- ifelse(dat$origin == "P1", P1X + sin(dat$rads) * dat$final_distance, 
    ifelse(dat$origin == "P2", P2X + sin(dat$rads) * dat$final_distance, ifelse(dat$origin == 
        "P3", P3X + sin(dat$rads) * dat$final_distance, ifelse(dat$origin == 
        "P4", P4X + sin(dat$rads) * dat$final_distance, ifelse(dat$origin == 
        "B13", XB13 + sin(dat$rads) * dat$final_distance, ifelse(dat$origin == 
        "B14", XB14 + sin(dat$rads) * dat$final_distance, ifelse(dat$origin == 
        "B18", XB18 + sin(dat$rads) * dat$final_distance, ifelse(dat$origin == 
        "B21", XB21 + sin(dat$rads) * dat$final_distance, NA))))))))

### calculate y-coordinates

dat$yfinal <- ifelse(dat$origin == "P1", P1Y + cos(dat$rads) * dat$final_distance, 
    ifelse(dat$origin == "P2", P2Y + cos(dat$rads) * dat$final_distance, ifelse(dat$origin == 
        "P3", P3Y + cos(dat$rads) * dat$final_distance, ifelse(dat$origin == 
        "P4", P4Y + cos(dat$rads) * dat$final_distance, ifelse(dat$origin == 
        "B13", YB13 + cos(dat$rads) * dat$final_distance, ifelse(dat$origin == 
        "B14", YB14 + cos(dat$rads) * dat$final_distance, ifelse(dat$origin == 
        "B18", YB18 + cos(dat$rads) * dat$final_distance, ifelse(dat$origin == 
        "B21", YB21 + cos(dat$rads) * dat$final_distance, NA))))))))

### calculate z
dat$zfinal <- ifelse(dat$origin == "P1", zP1 + sin(dat$radinc) * dat$dis, ifelse(dat$origin == 
    "P2", zP2 + sin(dat$radinc) * dat$dis, ifelse(dat$origin == "P3", zP3 + 
    sin(dat$radinc) * dat$dis, ifelse(dat$origin == "P4", zP4 + sin(dat$radinc) * 
    dat$dis, ifelse(dat$origin == "B13", zB13 - sin(dat$radinc) * dat$dis, ifelse(dat$origin == 
    "B14", zB14 + sin(dat$radinc) * dat$dis, ifelse(dat$origin == "B18", zB18 - 
    sin(dat$radinc) * dat$dis, ifelse(dat$origin == "B21", zB21 + sin(dat$radinc) * 
    dat$dis, NA))))))))
wndw <- owin(poly = list(x = c(P1X, P4X, P3X, P2X), y = c(P1Y, P4Y, P3Y, P2Y)))

datpic <- subset(dat, species == "Picea")
datpin <- subset(dat, species == "Pinus")
datpicppp <- ppp(x = datpic$xfinal, y = datpic$yfinal, window = wndw)
## Warning: 5 points were rejected as lying outside the specified window

datpinppp <- ppp(x = datpin$xfinal, y = datpin$yfinal, window = wndw)
## Warning: 1 point was rejected as lying outside the specified window
set.seed(234)
# env.pin <- envelope(datpinppp, fun = Lest, nsim = 99)
# 
# env.pic <- envelope(datpicppp, fun = Lest, nsim = 99)
# 
# plot(env.pin) ### The scatter of pinus is rather common, since it's
# covering with the envelope plot(env.pic) The scatter of picea is clearly
# clustered, since its y-value is way higher compared to the scatter of
# all pinus and the drawn envelope

# dat.ppp <- ppp(x = dat$xfinal, y = dat$yfinal, window = wndw)

# dat.ha2711 <- ha2711(x = ) plot(density(dat.ppp)) plot(Kest(dat.ppp))
# plot(Lest(dat.ppp)) enp<-envelope(dat.ppp,fun=Lest,nsim=99) plot(enp)
# create plot object for xfinal and yfinal

# p <- xyplot(dat$yfinal ~ dat$xfinal, dat = longleaf)

### print print(p) load libraries

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

### create spatial data object
dat.sp <- dat[complete.cases(dat), ]
coordinates(dat.sp) <- c("xfinal", "yfinal")
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(zfinal ~ 1, dat.sp, grdsp)
## [using ordinary kriging]


##

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

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

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