Skip to content
Snippets Groups Projects
Commit cae1c3f8 authored by fvafrcu's avatar fvafrcu
Browse files

Rename the package, split code files, spelling

parent fea69096
No related branches found
No related tags found
No related merge requests found
Showing with 226 additions and 378 deletions
Package: treePlotAreas
Package: treePlotArea
Title: Correction Factors for Tree Plot Areas Intersected by Stand
Boundaries
Version: 0.3.0.9000
......
# treePlotAreas 0.3.0.9000
# treePlotArea 0.3.0.9000
* FIXME
......
boundaries2coords <- function(x, impute_nook = FALSE) {
res <- x
res[["phia"]] <- geodatic2math(gon2degree(res[["spa_gon"]]))
res[["phie"]] <- geodatic2math(gon2degree(res[["spe_gon"]]))
res[["phik"]] <- geodatic2math(gon2degree(res[["spk_gon"]]))
res <- as.data.frame(cbind(polar2cartesian(res[["phia"]], res[["spa_m"]]),
polar2cartesian(res[["phik"]], res[["spk_m"]]),
polar2cartesian(res[["phie"]], res[["spe_m"]])))
names(res) <- c("x1", "y1", "x0", "y0", "x2", "y2")
if (isTRUE(impute_nook)) res <- impute_missing_nook(res)
return(res)
}
boundaries2polygons <- function(x, use_only_two_boundaries = TRUE, ...) {
current_boundaries <- get_current_boundaries(x)
if (isTRUE(use_only_two_boundaries) && nrow(current_boundaries) > 2)
current_boundaries <- delete_extraneous_boundaries(current_boundaries,
...)
coords <- boundaries2coords(current_boundaries)
coords_split <- split_coords(coords)
res <- apply(coords_split, 1, get_polygon, simplify = FALSE)
return(res)
}
#' Convert BWI Coordinates to Cartesian Coordinates
#'
#' BWI coordinates are measured in Gon eastward from north an cm distance.
#' We need cartesian coordinates for relational computations.
#' @param azimuth The azimuths, from north, eastern side, in gon.
#' @param distance The distances from the origin, typically measured in
#' centimeter.
#' @return Matrix of cartesian coordinates in the unit of \code{distance}.
#' @export
#' @keywords internal
#' @examples
#' a1 <- c(0, 100)
#' d1 <- c(100, 200)
#' print(coords <- bwi2cartesian(a1, d1))
#' all.equal(coords, matrix(c(0, 100, 200, 0), nrow = 2, byrow = TRUE),
#' check.attributes = FALSE)
bwi2cartesian <- function(azimuth, distance) {
res <- polar2cartesian(geodatic2math(gon2degree(azimuth)), distance)
return(res)
}
circle2polygon <- function(origin = c(x = 0, y = 0), r = 1, n = 400) {
circle <- NULL
for (i in seq(-r, r, length.out = n / 2)) {
si <- secant_intersections(Inf, i, r)
circle <- rbind(circle, c(si[1, TRUE], si[2, TRUE]))
}
circle <- rbind(circle[TRUE, c(1, 2)],
circle[order(circle[TRUE, 3], decreasing = TRUE), c(3, 4)])
circle[TRUE, "x"] <- circle[TRUE, "x"] + origin["x"]
circle[TRUE, "y"] <- circle[TRUE, "y"] + origin["y"]
circle <- rbind(circle, circle[1, TRUE])
return(circle)
}
#' Convert a Distance Into a Minimal Radius
#'
#' For a measure Distance, the tree will have to have at least that radius to be
#' part of the sample.
#' This is just a little helper for the field measurements.
#' @param r Boundary radius in meter.
#' @return Minimum dbh in cm.
#' @export
#' @keywords internal
#' @examples
#' get_min_dbh(12.5)
get_min_dbh <- function(r) {
res <- r * 4
attr(res, "unit") <- "cm"
return(res)
}
#' Get a Theoretical Maximum Distance for a Tree
#'
#' Maximum distance is of interest as boundaries that are more than double that
#' distance away are of no interest.
#' @return A theoretical maximum distance in centimeter. Based on the assumption
#' that trees have a maximum dbh of 200 cm.
#' @export
#' @keywords internal
#' @examples
#' get_r_max()
get_r_max <- function() {
# Unit needs to be centimeter as this is the way the horizontal distance is
# recorded.
# A tree with a dbh of 2000 millimeter. We don't expect any more.
res <- get_boundary_radius(2000, unit = "cm")
return(res)
}
gon2degree <- function(x) return(x * .9)
geodatic2math <- function(x, gon = FALSE) {
if (isTRUE(gon)) {
max <- 400
} else {
max <- 360
}
res <- max - x
res <- ifelse(res < 3/4 * max, res + 90, res - 3/4 * max)
return(res)
}
deg2rad <- function(x) {
res <- x * pi / 180
return(res)
}
polar2cartesian <- function(phi, distance) {
x <- cos(deg2rad(phi)) * distance
y <- sin(deg2rad(phi)) * distance
res <- cbind(x = x, y = y)
return(res)
}
#' Convert BWI Coordinates to Cartesian Coordinates
#'
#' BWI coordinates are measured in Gon eastward from north an cm distance.
#' We need cartesian coordinates for relational computations.
#' @param azimuth The azimuths, from north, eastern side, in gon.
#' @param distance The distances from the origin, typically measured in
#' centimeter.
#' @return Matrix of cartesian coordinates in the unit of \code{distance}.
#' @export
#' @keywords internal
#' @examples
#' a1 <- c(0, 100)
#' d1 <- c(100, 200)
#' print(coords <- bwi2cartesian(a1, d1))
#' all.equal(coords, matrix(c(0, 100, 200, 0), nrow = 2, byrow = TRUE),
#' check.attributes = FALSE)
bwi2cartesian <- function(azimuth, distance) {
res <- polar2cartesian(geodatic2math(gon2degree(azimuth)), distance)
return(res)
}
impute_missing_nook <- function(x) {
index <- is.na(x[["x0"]])
x[index, "x0"] <- (x[["x1"]][index] + x[["x2"]][index]) / 2
x[index, "y0"] <- (x[["y1"]][index] + x[["y2"]][index]) / 2
return(x)
}
boundaries2coords <- function(x, impute_nook = FALSE) {
res <- x
res[["phia"]] <- geodatic2math(gon2degree(res[["spa_gon"]]))
res[["phie"]] <- geodatic2math(gon2degree(res[["spe_gon"]]))
res[["phik"]] <- geodatic2math(gon2degree(res[["spk_gon"]]))
res <- as.data.frame(cbind(polar2cartesian(res[["phia"]], res[["spa_m"]]),
polar2cartesian(res[["phik"]], res[["spk_m"]]),
polar2cartesian(res[["phie"]], res[["spe_m"]])))
names(res) <- c("x1", "y1", "x0", "y0", "x2", "y2")
if (isTRUE(impute_nook)) res <- impute_missing_nook(res)
return(res)
}
project_along_r <- function(x, y, r) {
z <- create_xy()
if (x == 0) {
z["x"] <- 0
z["y"] <- r
} else {
b <- as.numeric(y / x)
z["x"] <- sqrt(1 / (b^2 + 1) * r^2)
z["y"] <- b * z["x"]
}
return(z)
}
vector_length <- function(v) {
z <- sqrt(sum(as.numeric(v)^2))
return(z)
}
#' Get a Slope–intercept Form from a Two-point Form of an Equation
#'
#' Two-point from is often seen in BWI data, we want to get an equation of form
#' \emph{y = a + bx}.
#' @param p1 The first point (x, y).
#' @param p2 The second point (x, y).
#' @return A named vector with intercept ["a"] and slope ["b"].
#' If both points have the same value for x, no function exists. Then
#' the intercept is \code{\link{NA}} and the slope gives the value of x.
#' @export
#' @keywords internal
#' @examples
#' points2equation(c(0, 4), c(1, 5))
points2equation <- function(p1, p2 = c(0, 0)) {
x1 <- as.numeric(p1[1])
y1 <- as.numeric(p1[2])
x2 <- as.numeric(p2[1])
y2 <- as.numeric(p2[2])
if (isTRUE(all.equal(x2, x1))) {
b <- x2
a <- Inf
} else {
b <- (y2 - y1) / (x2 - x1)
a <- y1 - x1 * b
}
res <- c(a = a, b = b)
return(res)
}
create_xy <- function() {
z <- vector("numeric", length = 2)
names(z) <- c("x", "y")
return(z)
}
has_nook <- function(boundary) {
res <- unlist(!is.na(unlist(boundary["x0"])))
return(res)
}
is_boundary_encloses_origin <- function(boundary) {
if (has_nook(boundary)) {
tri <- create_triangle(boundary, r = 1000 * get_r_max())
# plot(sf::st_polygon(list(tri)), add = TRUE)
res <- sf::st_contains(sf::st_polygon(list(tri)),
sf::st_point(c(0, 0)),
sparse = FALSE)
} else {
res <- FALSE
}
return(as.vector(res))
}
do_split <- function(boundary) {
res <- is_boundary_encloses_origin(boundary)
return(res)
}
split_coord <- function(coord) {
# split a triple of coords into to lines
if (isTRUE(do_split(coord))) {
res <- data.frame(x1 = c(coord[["x1"]], coord[["x0"]]),
y1 = c(coord[["y1"]], coord[["y0"]]),
x0 = c(NA, NA),
y0 = c(NA, NA),
x2 = c(coord[["x0"]], coord[["x2"]]),
y2 = c(coord[["y0"]], coord[["y2"]])
)
} else {
res <- coord
}
return(res)
}
split_coords <- function(coords) {
res <- do.call("rbind", apply(coords, 1, split_coord, simplify = FALSE))
return(res)
}
#' FIXME:
#'
#' @param a FIXME:
#' @param b FIXME:
#' @param r FIXME:
#' @return FIXME:
#' @export
#' @keywords internal
#' @examples
#' plot(0, 0, col = "red", pch = "+",
#' xlim = c(-2, 2),
#' ylim = c(-2, 2))
#' for (i in seq(-1, 1, by = 0.01)) {
#' points(secant_intersections(Inf, i, 1), pch = "+")
#' }
secant_intersections <- function(a, b, r) {
if (is.infinite(a)) {
x1 <- x2 <- as.numeric(b)
y1 <- sqrt(abs(r^2 - b^2))
y2 <- -y1
} else {
p <- (2 * a * b) / (b^2 + 1)
q <- (a^2 - r^2) / (b^2 + 1)
x1 <- - p/2 + sqrt((p/2)^2 - q)
y1 <- a + b * x1
x2 <- - p/2 - sqrt((p/2)^2 - q)
y2 <- a + b * x2
}
res <- matrix(c(x1, y1, x2, y2), nrow = 2, byrow = TRUE,
dimnames = list(c("1", "2"), c("x", "y")))
return(res)
}
#' FIXME:
#'
#' @param b FIXME:
#' @param xy FIXME:
#' @return FIXME:
#' @export
#' @keywords internal
#' @examples
#' orthogonal(1, c(x = 0, y = 0))
#' orthogonal(0, c(x = 4, y = 0))
orthogonal <- function(b, xy) {
if (identical(unname(b), 0)) {
res <- c(a = Inf, b = unname(xy["x"]))
} else {
b1 <- - 1 / unname(b)
a <- b1 * unname(xy["x"]) + unname(xy["y"])
res <- c(a = a, b = b1)
}
return(res)
}
create_tetragon <- function(xy1xy2, verbose = FALSE, r = 2 * get_r_max()) {
if(isTRUE(verbose)) message("# is one line -> need tetragon")
straight_line <- points2equation(c(xy1xy2["x1"], xy1xy2["y1"]),
c(xy1xy2["x2"], xy1xy2["y2"]))
si <- secant_intersections(straight_line["a"], straight_line["b"], r)
length_of_chord <- vector_length(c(si[1, "x"] - si[2, "x"],
si[1, "y"] - si[2, "y"]))
o <- orthogonal(straight_line["b"], c(x = si[1, "x"], y = si[1, "y"]))
b <- o["b"]
dx <- sqrt(length_of_chord^2 / (b^2 + 1))
dy <- b * dx
xyplus <- si[1, TRUE] + c(dx, dy)
xyminus <- si[1, TRUE] + c(-dx, -dy)
corners_away_from_origin <- NULL
if (vector_length(xyplus) > vector_length(xyminus)) {
corners_away_from_origin <- matrix(c(si[2, TRUE] + c(dx, dy), xyplus),
byrow = TRUE, ncol = 2,
dimnames = list(c("3", "4"),
c("x", "y")))
} else {
corners_away_from_origin <- matrix(c(si[2, TRUE] - c(dx, dy), xyminus),
byrow = TRUE, ncol = 2,
dimnames = list(c("3", "4"),
c("x", "y")))
}
tetragon <- rbind(si, corners_away_from_origin, si[1, TRUE])
res <- tetragon
return(res)
}
create_triangle <- function(xy1xy0xy2, verbose = FALSE, r = 2 * get_r_max()) {
if(isTRUE(verbose)) message("# is acute angle -> triangle will do")
# first side
straight_line <- points2equation(c(xy1xy0xy2["x1"], xy1xy0xy2["y1"]),
c(xy1xy0xy2["x0"], xy1xy0xy2["y0"]))
si <- secant_intersections(straight_line["a"], straight_line["b"], r)
s10 <- si[1, TRUE] - c(xy1xy0xy2[["x0"]], xy1xy0xy2[["y0"]])
s11 <- si[1, TRUE] - c(xy1xy0xy2[["x1"]], xy1xy0xy2[["y1"]])
if(vector_length(s11) < vector_length(s10)) {
i1 <- si[1, TRUE]
} else {
i1 <- si[2, TRUE]
}
# second side
straight_line <- points2equation(c(xy1xy0xy2["x2"], xy1xy0xy2["y2"]),
c(xy1xy0xy2["x0"], xy1xy0xy2["y0"]))
si <- secant_intersections(straight_line["a"], straight_line["b"], r)
s10 <- si[1, TRUE] - c(xy1xy0xy2[["x0"]], xy1xy0xy2[["y0"]])
s11 <- si[1, TRUE] - c(xy1xy0xy2[["x2"]], xy1xy0xy2[["y2"]])
if(vector_length(s11) < vector_length(s10)) {
i2 <- si[1, TRUE]
} else {
i2 <- si[2, TRUE]
}
res <- rbind(
c(xy1xy0xy2[["x0"]], xy1xy0xy2[["y0"]]),
i1, i2,
c(xy1xy0xy2[["x0"]], xy1xy0xy2[["y0"]])
)
return(res)
}
create_pentagon <- function(xy1xy0xy2, verbose = FALSE, r = 2 * get_r_max()) {
triangle <- create_triangle(xy1xy0xy2, verbose = FALSE, r = r)
xy1xy2 <- c(triangle[2,], triangle[3,])
names(xy1xy2) <- c("x1", "y1", "x2", "y2")
tetragon <- create_tetragon(xy1xy2, verbose = FALSE, r = r)
# insert nook in between rows 1 and 2 of tetragon:
res <- do.call("rbind", list(tetragon[1, TRUE],
triangle[1, TRUE],
tetragon[2:5, TRUE])
)
return(res)
}
get_polygon <- function(coord, r = 2 * get_r_max()) {
if (is.matrix(coord)) {
cv <- as.vector(coord)
names(cv) <- colnames(coord)
} else {
cv <- coord
}
if (has_nook(cv)) {
res <- create_pentagon(xy1xy0xy2 = cv, verbose = FALSE, r = r)
} else {
res <- create_tetragon(xy1xy2 = cv, verbose = FALSE, r = r)
}
return(res)
}
circle2polygon <- function(origin = c(x = 0, y = 0), r = 1, n = 400) {
circle <- NULL
for (i in seq(-r, r, length.out = n / 2)) {
si <- secant_intersections(Inf, i, r)
circle <- rbind(circle, c(si[1, TRUE], si[2, TRUE]))
}
circle <- rbind(circle[TRUE, c(1, 2)],
circle[order(circle[TRUE, 3], decreasing = TRUE), c(3, 4)])
circle[TRUE, "x"] <- circle[TRUE, "x"] + origin["x"]
circle[TRUE, "y"] <- circle[TRUE, "y"] + origin["y"]
circle <- rbind(circle, circle[1, TRUE])
return(circle)
}
get_partial_circle <- function(circle, boundaries) {
st_boundaries <- sf::st_polygon(boundaries)
res <- sf::st_polygon(list(circle))
for (b in st_boundaries) {
res <- sf::st_difference(res, sf::st_polygon(list(b)))
}
return(res)
}
tree2polygon <- function(tree) {
origin <- polar2cartesian(geodatic2math(gon2degree(tree[["azi"]])),
tree[["hori"]])
boundary_radius <- get_boundary_radius(tree[["m_bhd"]], unit = "cm")
res <- circle2polygon(c(x = tree[["x"]], y = tree[["y"]]), boundary_radius)
return(res)
}
create_pentagon <- function(xy1xy0xy2, verbose = FALSE, r = 2 * get_r_max()) {
triangle <- create_triangle(xy1xy0xy2, verbose = FALSE, r = r)
xy1xy2 <- c(triangle[2,], triangle[3,])
names(xy1xy2) <- c("x1", "y1", "x2", "y2")
tetragon <- create_tetragon(xy1xy2, verbose = FALSE, r = r)
# insert nook in between rows 1 and 2 of tetragon:
res <- do.call("rbind", list(tetragon[1, TRUE],
triangle[1, TRUE],
tetragon[2:5, TRUE])
)
return(res)
}
create_tetragon <- function(xy1xy2, verbose = FALSE, r = 2 * get_r_max()) {
if(isTRUE(verbose)) message("# is one line -> need tetragon")
straight_line <- points2equation(c(xy1xy2["x1"], xy1xy2["y1"]),
c(xy1xy2["x2"], xy1xy2["y2"]))
si <- secant_intersections(straight_line["a"], straight_line["b"], r)
length_of_chord <- vector_length(c(si[1, "x"] - si[2, "x"],
si[1, "y"] - si[2, "y"]))
o <- orthogonal(straight_line["b"], c(x = si[1, "x"], y = si[1, "y"]))
b <- o["b"]
dx <- sqrt(length_of_chord^2 / (b^2 + 1))
dy <- b * dx
xyplus <- si[1, TRUE] + c(dx, dy)
xyminus <- si[1, TRUE] + c(-dx, -dy)
corners_away_from_origin <- NULL
if (vector_length(xyplus) > vector_length(xyminus)) {
corners_away_from_origin <- matrix(c(si[2, TRUE] + c(dx, dy), xyplus),
byrow = TRUE, ncol = 2,
dimnames = list(c("3", "4"),
c("x", "y")))
} else {
corners_away_from_origin <- matrix(c(si[2, TRUE] - c(dx, dy), xyminus),
byrow = TRUE, ncol = 2,
dimnames = list(c("3", "4"),
c("x", "y")))
}
tetragon <- rbind(si, corners_away_from_origin, si[1, TRUE])
res <- tetragon
return(res)
}
create_triangle <- function(xy1xy0xy2, verbose = FALSE, r = 2 * get_r_max()) {
if(isTRUE(verbose)) message("# is acute angle -> triangle will do")
# first side
straight_line <- points2equation(c(xy1xy0xy2["x1"], xy1xy0xy2["y1"]),
c(xy1xy0xy2["x0"], xy1xy0xy2["y0"]))
si <- secant_intersections(straight_line["a"], straight_line["b"], r)
s10 <- si[1, TRUE] - c(xy1xy0xy2[["x0"]], xy1xy0xy2[["y0"]])
s11 <- si[1, TRUE] - c(xy1xy0xy2[["x1"]], xy1xy0xy2[["y1"]])
if(vector_length(s11) < vector_length(s10)) {
i1 <- si[1, TRUE]
} else {
i1 <- si[2, TRUE]
}
# second side
straight_line <- points2equation(c(xy1xy0xy2["x2"], xy1xy0xy2["y2"]),
c(xy1xy0xy2["x0"], xy1xy0xy2["y0"]))
si <- secant_intersections(straight_line["a"], straight_line["b"], r)
s10 <- si[1, TRUE] - c(xy1xy0xy2[["x0"]], xy1xy0xy2[["y0"]])
s11 <- si[1, TRUE] - c(xy1xy0xy2[["x2"]], xy1xy0xy2[["y2"]])
if(vector_length(s11) < vector_length(s10)) {
i2 <- si[1, TRUE]
} else {
i2 <- si[2, TRUE]
}
res <- rbind(
c(xy1xy0xy2[["x0"]], xy1xy0xy2[["y0"]]),
i1, i2,
c(xy1xy0xy2[["x0"]], xy1xy0xy2[["y0"]])
)
return(res)
}
create_xy <- function() {
z <- vector("numeric", length = 2)
names(z) <- c("x", "y")
return(z)
}
deg2rad <- function(x) {
res <- x * pi / 180
return(res)
}
delete_extraneous_boundaries <- function(x, rart_var = "rart") {
if (nrow(x) > 2) {
if (1 %in% x[[rart_var]]) {
x <- x[x[[rart_var]] == 1, TRUE]
if (nrow(x) > 2) {
x <- x[1:2, TRUE]
}
} else {
x <- x[1:2, TRUE]
}
}
return(x)
}
do_split <- function(boundary) {
res <- is_boundary_encloses_origin(boundary)
return(res)
}
geodatic2math <- function(x, gon = FALSE) {
if (isTRUE(gon)) {
max <- 400
} else {
max <- 360
}
res <- max - x
res <- ifelse(res < 3/4 * max, res + 90, res - 3/4 * max)
return(res)
}
get_current_boundaries <- function(x, boundary_column = "rk") {
res <- x[x[[boundary_column]] < 9, TRUE]
return(res)
}
delete_extraneous_boundaries <- function(x, rart_var = "rart") {
if (nrow(x) > 2) {
if (1 %in% x[[rart_var]]) {
x <- x[x[[rart_var]] == 1, TRUE]
if (nrow(x) > 2) {
x <- x[1:2, TRUE]
}
} else {
x <- x[1:2, TRUE]
}
}
return(x)
}
boundaries2polygons <- function(x, use_only_two_boundaries = TRUE, ...) {
current_boundaries <- get_current_boundaries(x)
if (isTRUE(use_only_two_boundaries) && nrow(current_boundaries) > 2)
current_boundaries <- delete_extraneous_boundaries(current_boundaries,
...)
coords <- boundaries2coords(current_boundaries)
coords_split <- split_coords(coords)
res <- apply(coords_split, 1, get_polygon, simplify = FALSE)
return(res)
}
#' FIXME
#'
#' @param x FIXME
......@@ -46,4 +14,3 @@ get_boundaries_polygons <- function(x) {
return(res)
}
get_current_boundaries <- function(x, boundary_column = "rk") {
res <- x[x[[boundary_column]] < 9, TRUE]
return(res)
}
#' Convert a Distance Into a Minimal Radius
#'
#' For a measure Distance, the tree will have to have at least that radius to be
#' part of the sample.
#' This is just a little helper for the field measurements.
#' @param r Boundary radius in meter.
#' @return Minimum dbh in cm.
#' @export
#' @keywords internal
#' @examples
#' get_min_dbh(12.5)
get_min_dbh <- function(r) {
res <- r * 4
attr(res, "unit") <- "cm"
return(res)
}
get_partial_circle <- function(circle, boundaries) {
st_boundaries <- sf::st_polygon(boundaries)
res <- sf::st_polygon(list(circle))
for (b in st_boundaries) {
res <- sf::st_difference(res, sf::st_polygon(list(b)))
}
return(res)
}
get_polygon <- function(coord, r = 2 * get_r_max()) {
if (is.matrix(coord)) {
cv <- as.vector(coord)
names(cv) <- colnames(coord)
} else {
cv <- coord
}
if (has_nook(cv)) {
res <- create_pentagon(xy1xy0xy2 = cv, verbose = FALSE, r = r)
} else {
res <- create_tetragon(xy1xy2 = cv, verbose = FALSE, r = r)
}
return(res)
}
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment