Monthly Archives: November 2013

VpSpPhWS13-L04-Enten_UmR

The Following codes all together are parts of the calculation of the UTM-Coordinates of trees.

Setting the path and working directory

path <- "C:/Users/Al Bern/Desktop/Dropbox/Uni/Vertiefung/Baumaufnahme/CSV"
setwd(path)

Reading in the table with the base data

data_complete <- read.csv("C:/Users/Al Bern/Desktop/Dropbox/Uni/Vertiefung/Baumaufnahme/CSV/complete.data.csv")

Changing the angles to radiants

rad <- function(dir) {
    radiants <- dir * (pi/180)
    return(radiants)
}

Changing the angles of dir and inc to radiants; calculating the corrected direction; writing it into data_complete as new columns


dir_rads <- rad(data_complete$dir)
data_complete$dir_rads <- dir_rads

inc_rads <- rad(data_complete$inc)
data_complete$inc_rads <- inc_rads

dir_rads_correct <- rad((data_complete$dir + 1.68))
data_complete$dir_rads_correct <- dir_rads_correct

Calculating the complete distance (With cbh in m)

distance_comp <- function(inc_rads, dis, cbh) {
    distance_comp <- cos(inc_rads) * (dis + ((cbh/100)/(2 * pi)))
    return(distance_comp)
}

Writing the complete distance

dis_comp <- distance_comp(data_complete$inc_rads, data_complete$dis, data_complete$cbh)
data_complete$dis_comp <- dis_comp

Calculating the x and y deviation

x_dev <- function(dir_rads_correct, dis_comp) {
    x_dev <- (sin(dir_rads_correct)) * (dis_comp)
    return(x_dev)
}

y_dev <- function(dir_rads_correct, dis_comp) {
    y_dev <- (cos(dir_rads_correct)) * (dis_comp)
    return(y_dev)
}

Calculating x and y deviation and writing it into data_complete as new column

x_dev <- x_dev(data_complete$dir_rads_correct, data_complete$dis_comp)
data_complete$x_dev <- x_dev

y_dev <- y_dev(data_complete$dir_rads_correct, data_complete$dis_comp)
data_complete$y_dev <- y_dev

Defining the trees that operated as origins for other trees

B14_x <- 481036.5 + sin(2.8972466) * 9.474704
B13_x <- 481036.5 + sin(3.1764992) * 11.714926
B21_x <- 481036.5 + sin(1.8325957) * 14.004796
B18_x <- B21_x + sin(3.9095375) * 7.928809

B14_y <- 5645561 + cos(2.8972466) * 9.474704
B13_y <- 5645561 + cos(3.1764992) * 11.714926
B21_y <- 5645561 + cos(1.8325957) * 14.004796
B18_y <- B21_y + cos(3.9095375) * 7.928809

Defining the Origins

data_complete$x_origin <- ifelse(data_complete$origin == "P1", 481073, ifelse(data_complete$origin == 
    "P2", 481073 + sin(rad(275.68)) * 40, ifelse(data_complete$origin == "P4", 
    481073 + sin(rad(5 + 1.68)) * 40, ifelse(data_complete$origin == "P3", 481073 + 
        sin(rad(5 + 1.68)) * 40 - 40, ifelse(data_complete$origin == "B14", 
        B14_x, ifelse(data_complete$origin == "B13", B13_x, ifelse(data_complete$origin == 
            "B21", B21_x, ifelse(data_complete$origin == "B18", B18_x, NA))))))))
data_complete$y_origin <- ifelse(data_complete$origin == "P1", 5645524, ifelse(data_complete$origin == 
    "P2", 5645524 + cos(rad(275 + 1.68)) * 40, ifelse(data_complete$origin == 
    "P4", 5645524 + cos(rad(5 + 1.68)) * 40, ifelse(data_complete$origin == 
    "P3", 5645524 + cos(rad(275 + 1.68)) * 40 + 40, ifelse(data_complete$origin == 
    "B14", B14_y, ifelse(data_complete$origin == "B13", B13_y, ifelse(data_complete$origin == 
    "B21", B21_y, ifelse(data_complete$origin == "B18", B18_y, NA))))))))

Calculating and writing the coordinates into the table data_complete as new column

x_coor <- function(x_origin, x_dev) {
    x_coor <- x_origin + x_dev
    return(x_coor)
}

x_coor <- x_coor(data_complete$x_origin, data_complete$x_dev)
data_complete$x_coor <- x_coor

y_coor <- function(y_origin, y_dev) {
    y_coor <- y_origin + y_dev
    return(y_coor)
}

y_coor <- y_coor(data_complete$y_origin, data_complete$y_dev)
data_complete$y_coor <- y_coor

height calculation (origins)

p2_height <- ((tan(data_complete[c(37), 9]) * data_complete[c(37), 11]) + 296) - 
    ((tan(data_complete[c(38), 9]) * data_complete[c(38), 11]))

p4_height <- ((tan(data_complete[c(46), 9]) * data_complete[c(45), 11]) + 296) - 
    ((tan(data_complete[c(45), 9]) * data_complete[c(46), 11]))

p3_height <- ((tan(data_complete[c(66), 9]) * data_complete[c(66), 11]) + 309.2018) - 
    ((tan(data_complete[c(64), 9]) * data_complete[c(64), 11])) - ((tan(data_complete[c(65), 
    9]) * data_complete[c(65), 11]))

b21_height <- ((tan(data_complete[c(64), 9]) * data_complete[c(64), 11])) + 
    p3_height

b14_height <- ((tan(data_complete[c(56), 9]) * data_complete[c(56), 11])) + 
    p3_height

b13_height <- ((tan(data_complete[c(55), 9]) * data_complete[c(55), 11])) + 
    p3_height

b18_height <- ((tan(data_complete[c(60), 9]) * data_complete[c(60), 11])) + 
    b21_height

writing height calculation in new column

data_complete$z_origin <- ifelse(data_complete$origin == "P1", 296, ifelse(data_complete$origin == 
    "P2", p2_height, ifelse(data_complete$origin == "P4", p4_height, ifelse(data_complete$origin == 
    "P3", p3_height, ifelse(data_complete$origin == "B14", b14_height, ifelse(data_complete$origin == 
    "B13", b13_height, ifelse(data_complete$origin == "B21", b21_height, ifelse(data_complete$origin == 
    "B18", b18_height, NA))))))))

calculating height (trees)

data_complete$z_value <- (tan(data_complete[c(1:113), 9]) * data_complete[c(1:113), 
    11] + data_complete$z_origin)

load libraries

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

create spatial data object

dat.sp <- data_complete[complete.cases(data_complete), ]
coordinates(dat.sp) <- c("x_coor", "y_coor")
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(z_value ~ 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-17

VpSpPhWS13-L04-Shades_of_Geography

The following script calculates the x/y/z – coordinates and a DEM.

### load packages ###
library(gstat)
library(knitr)
library(latticeExtra)
library(spatstat)

### set path and read data ###

path <- "D:/SpPhyGeo_WS13/tree_data_burgwald_csv/"
setwd("D:/SpPhyGeo_WS13/tree_data_burgwald_csv/")

data1 <- read.csv2("complete.data1.csv")

species <- read.csv("species_inventory.csv")
names(species) <- c("treeID", "cbh", "spec")
data <- merge(data1, species)


################################################ FUNCTIONS ####################


### convert cbh from m to cm ###

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

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


#################################### change degrees into radiants ###

# inc

radinc <- function(inc) {
    inclii <- inc * (pi/180)
    return(inclii)
}
incli <- radinc(data$inc)
data$radinc <- incli

# dir (add the deviation from magnetic north)

raddir <- function(dir) {
    rad <- (dir + 1.68) * (pi/180)
    return(rad)
}
rads <- raddir(data$dir)
data$raddir <- rads


####################################### calculate the complete distance ###

# add the tree radius to the calculated air line


fulldist <- function(cbhm, dis, radinc) {
    comp_distel <- (cbhm/(2 * pi)) + (cos(radinc) * dis)
    return(comp_distel)
}

comp_dist <- fulldist(data$cbhm, data$dis, data$radinc)
data$fulldist <- comp_dist


################### set origins ###

P1_x <- 481073
P1_y <- 5645524
P1_z <- 296


P4_x <- P1_x + ((sin(6.68 * (pi/180))) * 40)
P4_y <- P1_y + ((cos(6.68 * (pi/180))) * 40)
P4_z <- 309.8


P3_x <- P4_x + ((sin(276.68 * (pi/180))) * 40)
P3_y <- P4_y + ((cos(276.68 * (pi/180))) * 40)
P3_z <- 309.5


P2_x <- P3_x + ((sin(186.68 * (pi/180))) * 40)
P2_y <- P3_y + ((cos(186.68 * (pi/180))) * 40)
P2_z <- 295.5

############################# calculate treeorigins ###

B13_x <- 481037.9 + sin(3.20582077) * (11.715509)
B14_x <- 481037.9 + sin(2.92656809) * (9.47474)
B21_x <- 481037.9 + sin(1.86191725) * (14.004838)
B18_x <- B21_x + sin(3.93885906) * (7.930687)


B13_y <- 5645568 + cos(3.20582077) * (11.715509)
B14_y <- 5645568 + cos(2.92656809) * (9.47474)
B21_y <- 5645568 + cos(1.86191725) * (14.004838)
B18_y <- B21_y + cos(3.93885906) * (7.930687)

B13_z <- P3_z + sin(-0.1134464) * 11.7
B14_z <- P3_z + sin(-0.03490659) * 9.42
B21_z <- P3_z + sin(0.01745329) * 13.73
B18_z <- B21_z + sin(-0.20943951) * 8.02


################################# calculate x - coordinates ###


data$x <- ifelse(data$origin == "P1", P1_x + sin(data$raddir) * data$fulldist, 
    ifelse(data$origin == "P2", P2_x + sin(data$raddir) * data$fulldist, ifelse(data$origin == 
        "P3", P3_x + sin(data$raddir) * data$fulldist, ifelse(data$origin == 
        "P4", P4_x + sin(data$raddir) * data$fulldist, ifelse(data$origin == 
        "B13", B13_x + sin(data$raddir) * data$fulldist, ifelse(data$origin == 
        "B14", B14_x + sin(data$raddir) * data$fulldist, ifelse(data$origin == 
        "B21", B21_x + sin(data$raddir) * data$fulldist, ifelse(data$origin == 
        "B18", B18_x + sin(data$raddir) * data$fulldist, NA))))))))

################################# calculate y - coordinates ###

data$y <- ifelse(data$origin == "P1", P1_y + cos(data$raddir) * data$fulldist, 
    ifelse(data$origin == "P2", P2_y + cos(data$raddir) * data$fulldist, ifelse(data$origin == 
        "P3", P3_y + cos(data$raddir) * data$fulldist, ifelse(data$origin == 
        "P4", P4_y + cos(data$raddir) * data$fulldist, ifelse(data$origin == 
        "B13", B13_y + cos(data$raddir) * data$fulldist, ifelse(data$origin == 
        "B14", B14_y + cos(data$raddir) * data$fulldist, ifelse(data$origin == 
        "B21", B21_y + cos(data$raddir) * data$fulldist, ifelse(data$origin == 
        "B18", B18_y + cos(data$raddir) * data$fulldist, NA))))))))




################################# calculate z - coordinates ###

data$z <- ifelse(data$origin == "P1", P1_z + sin(data$inc) * data$dis, ifelse(data$origin == 
    "P2", P2_z + sin(data$inc) * data$dis, ifelse(data$origin == "P3", P3_z + 
    sin(data$inc) * data$dis, ifelse(data$origin == "P4", P4_z + sin(data$inc) * 
    data$dis, ifelse(data$origin == "B13", B13_z + sin(data$inc) * data$dis, 
    ifelse(data$origin == "B14", B14_z + sin(data$inc) * data$dis, ifelse(data$origin == 
        "B21", B21_z + sin(data$inc) * data$dis, ifelse(data$origin == "B18", 
        B18_z + sin(data$inc) * data$dis, NA))))))))


### load libraries
library(sp)
library(latticeExtra)
library(automap)
library(raster)

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

## [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