From 44c3dbbb41d91213626e9e6c83f5cc5d7545dbe3 Mon Sep 17 00:00:00 2001 From: nicsoattc <nicolas.cattaneo@nibio.no> Date: Fri, 20 Oct 2023 14:40:20 +0200 Subject: [PATCH] update repo --- CodeR/03_Compute_angles_cone.r | 144 +++++++++++++++++++++++++++++++-- 1 file changed, 138 insertions(+), 6 deletions(-) diff --git a/CodeR/03_Compute_angles_cone.r b/CodeR/03_Compute_angles_cone.r index e408c75..eb0af4d 100644 --- a/CodeR/03_Compute_angles_cone.r +++ b/CodeR/03_Compute_angles_cone.r @@ -29,27 +29,159 @@ source(paste0(CodeR, "/Functions/Functions.r")) ## plot data - full path DataPaths <- list.files(DataR, ## pattern = "plot_", - pattern = "Plot_", - ## pattern = "merg", + ## pattern = "Plot_", + pattern = "merg", ## pattern = ".rds", ## pattern = "Tower", full.names = T) -Data <- read.las(DataPaths) -Data[, 1:3] +Data <- read.las(DataPaths[1]) +Data <- Data[, 1:3] + + + + +# Create a simple 3D arrow function +drawArrow <- function(start, end, col = "red", lwd = 5, arrowSize = 0.1) { + segments3d(rbind(start, end), col = col, lwd = lwd) + + # Compute the direction vector + dirVector <- end - start + + # Normalize the direction vector + dirVector <- dirVector / sqrt(sum(dirVector^2)) + + # Create the arrow head (as three segments emanating from the arrow tip) + segments3d(rbind(end, end - arrowSize * dirVector + c( arrowSize, 0, 0)), col = col, lwd = lwd) + segments3d(rbind(end, end - arrowSize * dirVector + c(-arrowSize, 0, 0)), col = col, lwd = lwd) + segments3d(rbind(end, end - arrowSize * dirVector + c(0, -arrowSize, 0)), col = col, lwd = lwd) +} + +# Example data +set.seed(123) +n <- 1000 +x <- runif(n, 500000, 505000) # UTM x-coordinates +y <- runif(n, 5000000, 5005000) # UTM y-coordinates +z <- runif(n, 0, 100) # just some random z values + +# Plot the point cloud +plot3d(x, y, z, type = "p", col = "blue") + +# Define the "Devise" point +devise <- c(502500, 5002500, 50) # For instance +spheres3d(devise, radius = 100) + + +# Compute the arrow length based on the data range +arrow_length <- max(c(diff(range(x)), diff(range(y)), diff(range(z)))) / 10 + +# South (S) +drawArrow(devise, devise + c(0, -arrow_length, 0), col="red") + +# West (W) +drawArrow(devise, devise + c(-arrow_length, 0, 0), col="green") + +# North (N) +drawArrow(devise, devise + c(0, arrow_length, 0), col="blue") + +# East (E) +drawArrow(devise, devise + c(arrow_length, 0, 0), col="yellow") + +# Southwest (SW) +drawArrow(devise, devise + c(-arrow_length, -arrow_length, 0), col="purple") + +# Southeast (SE) +drawArrow(devise, devise + c(arrow_length, -arrow_length, 0), col="cyan") + +# Northwest (NW) +drawArrow(devise, devise + c(-arrow_length, arrow_length, 0), col="magenta") + +# Northeast (NE) +drawArrow(devise, devise + c(arrow_length, arrow_length, 0), col="orange") + + + +getArrowTipCoords <- function(base, length, direction) { + if (direction == "S") { + return(base + c(0, -length, 0)) + } else if (direction == "W") { + return(base + c(-length, 0, 0)) + } else if (direction == "N") { + return(base + c(0, length, 0)) + } else if (direction == "E") { + return(base + c(length, 0, 0)) + } else if (direction == "SW") { + return(base + c(-length, -length, 0)) + } else if (direction == "SE") { + return(base + c(length, -length, 0)) + } else if (direction == "NW") { + return(base + c(-length, length, 0)) + } else if (direction == "NE") { + return(base + c(length, length, 0)) + } else { + stop("Invalid direction provided!") + } +} + +# Define the "Devise" point and arrow length as before +devise <- c(502500, 5002500, 50) # For instance +arrow_length <- max(c(diff(range(x)), diff(range(y)), diff(range(z)))) / 10 + +# Get the coordinates for each direction +s_tip <- getArrowTipCoords(devise, arrow_length, "S") +w_tip <- getArrowTipCoords(devise, arrow_length, "W") +n_tip <- getArrowTipCoords(devise, arrow_length, "N") +e_tip <- getArrowTipCoords(devise, arrow_length, "E") +sw_tip <- getArrowTipCoords(devise, arrow_length, "SW") +se_tip <- getArrowTipCoords(devise, arrow_length, "SE") +nw_tip <- getArrowTipCoords(devise, arrow_length, "NW") +ne_tip <- getArrowTipCoords(devise, arrow_length, "NE") + +# Print the coordinates +print(list(S = s_tip)) + + + + + + + ## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ +## Hi Nico, +## Assuming the coordinates in the PointCloud are UTM, +## the preferred location of the SIF sensor would be at + +## N 6694611.571 E 614675.0338 and height 315.8627 masl + +## This is a point in the middle of the SW-SE edge of the tower +## top, and the preferred viewing direction is more or less south +## since there is a dense group of trees there; the crucial +## question is the viewing angle. + +## Cheers, Holger + +Devise <- data.table(X = 614675.0338, + Y = 6694611.571, + Z = 315.8627) + +## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ Data <- Data[sample(1:nrow(Data), nrow(Data)*0.10), ] + +Data <- rbind(Data, Devise) + Data[, X := Scale(Data$X)] Data[, Y := Scale(Data$Y)] -Data[, Z := Scale(Data$Z)] Data[, X := X-mean(X)] Data[, Y := Y-mean(Y)] +Devise <- tail(Data, 1) -plot3d(Data, aspect = "iso", size = 0.5) +plot3d(Data, aspect = "iso", size = 0.2) pan3d(3) +spheres3d(Devise, radius = 0.5, col = "red") + ## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ## Create an artificial tower. Just the corners -- GitLab