Skip to content
GitLab
Explore
Sign in
Register
Primary navigation
Search or go to…
Project
H
Holger_Tower_POV
Manage
Activity
Members
Labels
Plan
Issues
Issue boards
Milestones
Wiki
Code
Merge requests
Repository
Branches
Commits
Tags
Repository graph
Compare revisions
Snippets
Build
Pipelines
Jobs
Pipeline schedules
Artifacts
Deploy
Releases
Package registry
Container registry
Model registry
Operate
Environments
Terraform modules
Monitor
Incidents
Analyze
Value stream analytics
Contributor analytics
CI/CD analytics
Repository analytics
Model experiments
Help
Help
Support
GitLab documentation
Compare GitLab plans
Community forum
Contribute to GitLab
Provide feedback
Keyboard shortcuts
?
Snippets
Groups
Projects
Show more breadcrumbs
Nicolas Cattaneo
Holger_Tower_POV
Commits
22a52677
Commit
22a52677
authored
1 year ago
by
nicoscattaneo
Browse files
Options
Downloads
Patches
Plain Diff
update repo
parent
44c3dbbb
Branches
Branches containing commit
Tags
Tags containing commit
No related merge requests found
Changes
3
Hide whitespace changes
Inline
Side-by-side
Showing
3 changed files
CodeR/03_Compute_angles_cone.r
+108
-117
108 additions, 117 deletions
CodeR/03_Compute_angles_cone.r
CodeR/04_Compute_angles_cone_clear.r
+147
-0
147 additions, 0 deletions
CodeR/04_Compute_angles_cone_clear.r
CodeR/04_Compute_angles_cone_clear.r~
+224
-0
224 additions, 0 deletions
CodeR/04_Compute_angles_cone_clear.r~
with
479 additions
and
117 deletions
CodeR/03_Compute_angles_cone.r
+
108
−
117
View file @
22a52677
...
@@ -18,13 +18,18 @@ library("data.table")
...
@@ -18,13 +18,18 @@ library("data.table")
## DataR <- paste0("/home/nica/Documents/Holger_tower_pov/DataR")
## DataR <- paste0("/home/nica/Documents/Holger_tower_pov/DataR")
## CodeR <- paste0("/home/nica/Documents/Holger_tower_pov/CodeR")
## CodeR <- paste0("/home/nica/Documents/Holger_tower_pov/CodeR")
DataR
<-
paste0
(
"/home/nibio/REPOS/holger_tower_pov/DataR"
)
## DataR <- paste0("/home/nibio/REPOS/holger_tower_pov/DataR")
CodeR
<-
paste0
(
"/home/nibio/REPOS/holger_tower_pov/CodeR"
)
## CodeR <- paste0("/home/nibio/REPOS/holger_tower_pov/CodeR")
CodeR
<-
"/home/pepito/Documents/holger_tower_pov/CodeR/"
DataR
<-
"/home/pepito/Desktop/Temp_Holger_Data/"
## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
## Functions
## Functions
source
(
paste0
(
CodeR
,
"/Functions/Functions.r"
))
source
(
paste0
(
CodeR
,
"/Functions/Functions.r"
))
objects
()
dist3d
## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
## plot data - full path
## plot data - full path
DataPaths
<-
list.files
(
DataR
,
DataPaths
<-
list.files
(
DataR
,
...
@@ -35,159 +40,145 @@ DataPaths <- list.files(DataR,
...
@@ -35,159 +40,145 @@ DataPaths <- list.files(DataR,
## pattern = "Tower",
## pattern = "Tower",
full.names
=
T
)
full.names
=
T
)
Data
<-
read.las
(
DataPaths
[
1
])
## Data <- read.las(DataPaths[1])
Data
<-
fst
::
read_fst
(
DataPaths
[
1
],
as.data.table
=
T
)
Data
<-
Data
[,
1
:
3
]
Data
<-
Data
[,
1
:
3
]
## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
## Hi Nico,
## Assuming the coordinates in the PointCloud are UTM,
## the preferred location of the SIF sensor would be at
# Create a simple 3D arrow function
## N 6694611.571 E 614675.0338 and height 315.8627 masl
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)
## This is a point in the middle of the SW-SE edge of the tower
drawArrow
(
devise
,
devise
+
c
(
arrow_length
,
-
arrow_length
,
0
),
col
=
"cyan"
)
## 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.
# Northwest (NW)
## Cheers, Holger
drawArrow
(
devise
,
devise
+
c
(
-
arrow_length
,
arrow_length
,
0
),
col
=
"magenta"
)
## Devise <- data.table(X = 614675.0338,
## Y = 6694611.571,
## Z = 315.8627)
# Northeast (NE)
## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
drawArrow
(
devise
,
devise
+
c
(
arrow_length
,
arrow_length
,
0
),
col
=
"orange"
)
Data
<-
Data
[
sample
(
1
:
nrow
(
Data
),
nrow
(
Data
)
*
0.10
),
]
gc
()
## Data <- rbind(Data, Devise)
Data
[,
X
:=
Scale
(
Data
$
X
)]
Data
[,
Y
:=
Scale
(
Data
$
Y
)]
## Devise <- tail(Data, 1)
getArrowTipCoords
<-
function
(
base
,
length
,
direction
)
{
plot3d
(
Data
,
aspect
=
"iso"
,
size
=
0.2
)
if
(
direction
==
"S"
)
{
pan3d
(
3
)
return
(
base
+
c
(
0
,
-
length
,
0
))
## spheres3d(Devise, radius = 0.5, col = "red")
}
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 location is not in the area mentioned in the email.
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"
)
## Select the top of the tower
w_tip
<-
getArrowTipCoords
(
devise
,
arrow_length
,
"W"
)
## ids <- plot3d(Data, aspect = "iso", size = 4)
n_tip
<-
getArrowTipCoords
(
devise
,
arrow_length
,
"N"
)
## pan3d(3)
e_tip
<-
getArrowTipCoords
(
devise
,
arrow_length
,
"E"
)
## TopTowerIndex <- selectpoints3d(ids["data"], value = FALSE)
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
## TopTowerIndex[, 2]
print
(
list
(
S
=
s_tip
))
## plot3d(Data[-TopTowerIndex[, 2], ], aspect = "iso", size = 3)
## pan3d(3)
## plot3d(Data[TopTowerIndex[, 2], ], aspect = "iso",
## size = 3, col = "red", add = T)
## plot3d(Data[TopTowerIndex[, 2], ], aspect = "iso",
## size = 5, col = "red")
## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
## TopTowerP <- Data[TopTowerIndex[, 2], ]
## saveRDS(TopTowerP, paste0(DataR, "TopTowerP.rds"))
## TopTowerP <- readRDS(paste0(DataR, "TopTowerP.rds"))
## plot3d(TopTowerP, aspect = "iso",
## size = 10, col = "red")
## ToWerCenter <- TopTowerP[10, ]
## spheres3d(ToWerCenter, radius = 0.1, col = "black")
## TopTowerSW <- TopTowerP[306, ]
## spheres3d(TopTowerSW, radius = 0.1, col = "green")
## TopTowerSE <- TopTowerP[168, ]
## spheres3d(TopTowerSE, radius = 0.1, col = "blue")
## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
## 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
cloud
<-
data.frame
(
x
=
runif
(
10000
,
-20
,
20
),
y
=
runif
(
10000
,
-20
,
20
),
z
=
runif
(
10000
,
0
,
40
))
## 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
# Define the TopTowerSW point
TopTowerSW
<-
data.frame
(
x
=
0.5
,
y
=
0.5
,
x
=
40
)
Devise
<-
data.table
(
X
=
614675.0338
,
# Function to compute the distance between points and the line
Y
=
6694611.571
,
distance_to_line
<-
function
(
cloud
,
TopTowerSW
,
theta_deg
)
{
Z
=
315.8627
)
# Convert angle to radians
theta_rad
<-
theta_deg
*
pi
/
180
# Direction vector for the line
d
<-
c
(
0
,
-1
,
tan
(
theta_rad
))
# Filter the cloud to retain only points to the south of TopTowerSW
cloud_south
<-
cloud
[
cloud
$
y
<
TopTowerSW
$
y
,
]
# Convert both data.frames to matrices for easier computation
cloud_matrix
<-
as.matrix
(
cloud_south
)
TopTowerSW_matrix
<-
matrix
(
rep
(
as.matrix
(
TopTowerSW
),
nrow
(
cloud_south
)),
ncol
=
3
,
byrow
=
TRUE
)
# Compute the position vectors relative to TopTowerSW for all points
r
<-
cloud_matrix
-
TopTowerSW_matrix
# Compute the cross product between r and d for all points
cross_product
<-
cbind
(
r
[,
2
]
*
d
[
3
]
-
r
[,
3
]
*
d
[
2
],
r
[,
3
]
*
d
[
1
]
-
r
[,
1
]
*
d
[
3
],
r
[,
1
]
*
d
[
2
]
-
r
[,
2
]
*
d
[
1
])
# Compute the distance using the formula mentioned
distances
<-
sqrt
(
rowSums
(
cross_product
^
2
))
/
sqrt
(
sum
(
d
^
2
))
return
(
distances
)
}
#
#
@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
#
Example usage: Compute the distance with a vertical angle of 45 degrees
Data
<-
Data
[
sample
(
1
:
nrow
(
Data
),
nrow
(
Data
)
*
0.10
),
]
distances
<-
distance_to_line
(
cloud
,
TopTowerSW
,
-60
)
Data
<-
rbind
(
Data
,
Devise
)
Data
[,
X
:=
Scale
(
Data
$
X
)]
plot3d
(
cloud
,
size
=
5
)
Data
[,
Y
:=
Scale
(
Data
$
Y
)]
plot3d
(
TopTowerSW
,
size
=
10
,
add
=
T
,
col
=
"red"
)
Data
[,
X
:=
X
-
mean
(
X
)]
Puntitos
<-
as.numeric
(
names
(
distances
[
distances
<
4
]))
Data
[,
Y
:=
Y
-
mean
(
Y
)]
plot3d
(
cloud
[
Puntitos
,
],
size
=
10
,
col
=
"blue"
,
add
=
T
)
Devise
<-
tail
(
Data
,
1
)
compute_endpoint
<-
function
(
TopTowerSW
,
cloud
,
theta_deg
)
{
# Direction vector for the line based on the vertical angle
d
<-
c
(
0
,
-1
,
tan
(
theta_deg
*
pi
/
180
))
# Calculate the endpoint of the line based on the vertical angle
distance_to_edge
<-
TopTowerSW
$
y
-
min
(
cloud
$
y
)
z_change_due_to_angle
<-
distance_to_edge
*
d
[
3
]
end_point
<-
c
(
TopTowerSW
$
x
,
TopTowerSW
$
y
-
distance_to_edge
,
TopTowerSW
$
z
+
z_change_due_to_angle
)
return
(
end_point
)
}
plot3d
(
Data
,
aspect
=
"iso"
,
size
=
0.2
)
# For demonstration: Compute endpoint with a vertical angle of 30 degrees
pan3d
(
3
)
theta_deg
<-
-60
spheres3d
(
Devise
,
radius
=
0.5
,
col
=
"red"
)
end_point
<-
compute_endpoint
(
TopTowerSW
,
cloud
,
theta_deg
)
segments3d
(
rbind
(
TopTowerSW
,
end_point
),
col
=
'red'
)
## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
## Create an artificial tower. Just the corners
## Create an artificial tower. Just the corners
CornA
<-
data.table
(
X
=
3
,
Y
=
3
,
Z
=
40
)
CornA
<-
data.table
(
X
=
3
,
Y
=
3
,
Z
=
40
)
## CornA <- Data[identify3d(Data, n = 1), ]
## CornA <- Data[identify3d(Data, n = 1), ]
spheres3d
(
CornA
,
radius
=
1
,
col
=
"yellow"
)
##
spheres3d(CornA, radius = 1, col = "yellow")
## Left click on a ground point in the point cloud
## Left click on a ground point in the point cloud
## Right click to scape
## Right click to scape
...
...
This diff is collapsed.
Click to expand it.
CodeR/04_Compute_angles_cone_clear.r
0 → 100644
+
147
−
0
View file @
22a52677
rm
(
list
=
ls
())
gc
()
## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
## packages functions
rm
(
list
=
ls
())
library
(
"rlas"
)
library
(
"rgl"
)
library
(
"stringr"
)
library
(
"VoxR"
)
library
(
"data.table"
)
## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
## DataR and CodeR
## DataR <- paste0("/home/nibio/REPOS/holger_tower_pov/DataR")
## DataR <- paste0("/home/pepito/Documents/holger_tower_pov/DataR")
## DataR <- paste0("/home/nica/Documents/Holger_tower_pov/DataR")
## CodeR <- paste0("/home/nica/Documents/Holger_tower_pov/CodeR")
## DataR <- paste0("/home/nibio/REPOS/holger_tower_pov/DataR")
## CodeR <- paste0("/home/nibio/REPOS/holger_tower_pov/CodeR")
CodeR
<-
"/home/pepito/Documents/holger_tower_pov/CodeR/"
DataR
<-
"/home/pepito/Desktop/Temp_Holger_Data/"
## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
## Functions
source
(
paste0
(
CodeR
,
"/Functions/Functions.r"
))
## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
## plot data - full path
DataPaths
<-
list.files
(
DataR
,
## pattern = "plot_",
## pattern = "Plot_",
pattern
=
"merg"
,
## pattern = ".rds",
## pattern = "Tower",
full.names
=
T
)
## Data <- read.las(DataPaths[1])
Data
<-
fst
::
read_fst
(
DataPaths
[
1
],
as.data.table
=
T
)
Data
<-
Data
[,
1
:
3
]
## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
Data
<-
Data
[
sample
(
1
:
nrow
(
Data
),
nrow
(
Data
)
*
0.10
),
]
gc
()
Data
[,
X
:=
Scale
(
Data
$
X
)]
Data
[,
Y
:=
Scale
(
Data
$
Y
)]
## plot3d(Data, aspect = "iso", size = 0.2)
## plot3d(Data, aspect = "iso", size = 3)
## pan3d(3)
## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
## Devise location is not in the area mentioned in the email.
## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
## A point at the topo ofthe tower, SW edge
TopTowerSW
<-
data.table
(
X
=
173.4451
,
Y
=
104.955
,
Z
=
322.6828
)
TopTowerSE
<-
data.table
(
X
=
175.3196
,
Y
=
105.2212
,
Z
=
322.7785
)
## spheres3d(TopTowerSE, radius = 0.3, col = "yellow")
## spheres3d(TopTowerSW, radius = 0.3, col = "red")
## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
# Function to compute the distance between points and the line
distance2line
<-
function
(
cloud
,
TopTowerP
,
theta_deg
)
{
theta_rad
<-
theta_deg
*
pi
/
180
d
<-
c
(
0
,
-1
,
tan
(
theta_rad
))
cloud_south_indices
<-
which
(
cloud
$
Y
<
TopTowerP
$
Y
)
cloud_south
<-
cloud
[
cloud_south_indices
,
]
cloud_matrix
<-
as.matrix
(
cloud_south
)
TopTowerP_matrix
<-
matrix
(
rep
(
as.matrix
(
TopTowerP
),
nrow
(
cloud_south
)),
ncol
=
3
,
byrow
=
TRUE
)
r
<-
cloud_matrix
-
TopTowerP_matrix
cross_product
<-
cbind
(
r
[,
2
]
*
d
[
3
]
-
r
[,
3
]
*
d
[
2
],
r
[,
3
]
*
d
[
1
]
-
r
[,
1
]
*
d
[
3
],
r
[,
1
]
*
d
[
2
]
-
r
[,
2
]
*
d
[
1
])
distances
<-
sqrt
(
rowSums
(
cross_product
^
2
))
/
sqrt
(
sum
(
d
^
2
))
Result
<-
data.table
(
distances
=
distances
,
ids
=
cloud_south_indices
)
gc
()
return
(
Result
)
}
# Example usage: Compute the distance with a vertical angle of 45 degrees
distances
<-
distance2line
(
Data
,
TopTowerSW
,
-40
)
## plot3d(Data, aspect = "iso", size = 0.2)
## spheres3d(TopTowerSW, radius = 1, col = "red")
## plot3d(Data[distances[distances < 40, ]$ids, ], size = 10, col = "red", add = T)
## pan3d(3)
## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
Data
[
distances
$
ids
,
CO
:=
distances
$
distances
]
## plot3d(Data, aspect = "iso", size = 0.5, col = "grey21")
## plot3d(Data[CO < 5 , ], aspect = "iso", size = 4, col = "red", add = T)
## spheres3d(CornA, radius = 1, col = "yellow")
## spheres3d(CornB, radius = 1, col = "red")
## pan3d(3)
Data
[
!
is.na
(
CO
),
H
:=
dist3d2
(
X
,
Y
,
Z
,
TopTowerSW
$
X
,
TopTowerSW
$
Y
,
TopTowerSW
$
Z
),
by
=
seq_len
(
nrow
(
Data
))]
Data
[,
Angle
:=
CO
/
H
]
Data
[,
Angle
:=
(
asin
(
Angle
))
*
180
/
pi
]
OpeningAngl
<-
20
plot3d
(
Data
[
Angle
<
OpeningAngl
,
],
aspect
=
"iso"
,
size
=
1.5
,
col
=
"red"
)
spheres3d
(
TopTowerSW
,
radius
=
1
,
col
=
"yellow"
)
plot3d
(
Data
[
Angle
>
OpeningAngl
,
],
aspect
=
"iso"
,
size
=
0.5
,
col
=
"grey21"
)
spheres3d
(
TopTowerSW
,
radius
=
1
,
col
=
"red"
)
plot3d
(
Data
[
Angle
<
OpeningAngl
,
],
aspect
=
"iso"
,
size
=
0.5
,
col
=
"red"
,
add
=
T
)
pan3d
(
3
)
## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
This diff is collapsed.
Click to expand it.
CodeR/04_Compute_angles_cone_clear.r~
0 → 100644
+
224
−
0
View file @
22a52677
rm(list=ls())
gc()
## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
## packages functions
rm(list=ls())
library("rlas")
library("rgl")
library("stringr")
library("VoxR")
library("data.table")
## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
## DataR and CodeR
## DataR <- paste0("/home/nibio/REPOS/holger_tower_pov/DataR")
## DataR <- paste0("/home/pepito/Documents/holger_tower_pov/DataR")
## DataR <- paste0("/home/nica/Documents/Holger_tower_pov/DataR")
## CodeR <- paste0("/home/nica/Documents/Holger_tower_pov/CodeR")
## DataR <- paste0("/home/nibio/REPOS/holger_tower_pov/DataR")
## CodeR <- paste0("/home/nibio/REPOS/holger_tower_pov/CodeR")
CodeR <- "/home/pepito/Documents/holger_tower_pov/CodeR/"
DataR <- "/home/pepito/Desktop/Temp_Holger_Data/"
## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
## Functions
source(paste0(CodeR, "/Functions/Functions.r"))
objects()
dist3d
## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
## plot data - full path
DataPaths <- list.files(DataR,
## pattern = "plot_",
## pattern = "Plot_",
pattern = "merg",
## pattern = ".rds",
## pattern = "Tower",
full.names = T)
## Data <- read.las(DataPaths[1])
Data <- fst::read_fst(DataPaths[1], as.data.table = T)
Data <- Data[, 1:3]
## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
## 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), ]
gc()
## Data <- rbind(Data, Devise)
Data[, X := Scale(Data$X)]
Data[, Y := Scale(Data$Y)]
## Devise <- tail(Data, 1)
plot3d(Data, aspect = "iso", size = 0.2)
pan3d(3)
## spheres3d(Devise, radius = 0.5, col = "red")
## Devise location is not in the area mentioned in the email.
## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
## Select the top of the tower
## ids <- plot3d(Data, aspect = "iso", size = 4)
## pan3d(3)
## TopTowerIndex <- selectpoints3d(ids["data"], value = FALSE)
## TopTowerIndex[, 2]
## plot3d(Data[-TopTowerIndex[, 2], ], aspect = "iso", size = 3)
## pan3d(3)
## plot3d(Data[TopTowerIndex[, 2], ], aspect = "iso",
## size = 3, col = "red", add = T)
## plot3d(Data[TopTowerIndex[, 2], ], aspect = "iso",
## size = 5, col = "red")
## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
## TopTowerP <- Data[TopTowerIndex[, 2], ]
## saveRDS(TopTowerP, paste0(DataR, "TopTowerP.rds"))
## TopTowerP <- readRDS(paste0(DataR, "TopTowerP.rds"))
## plot3d(TopTowerP, aspect = "iso",
## size = 10, col = "red")
## ToWerCenter <- TopTowerP[10, ]
## spheres3d(ToWerCenter, radius = 0.1, col = "black")
## TopTowerSW <- TopTowerP[306, ]
## spheres3d(TopTowerSW, radius = 0.1, col = "green")
## TopTowerSE <- TopTowerP[168, ]
## spheres3d(TopTowerSE, radius = 0.1, col = "blue")
cloud <- data.frame(x = runif(10000, -20, 20), y = runif(10000, -20, 20), z = runif(10000, 0, 40))
# Define the TopTowerSW point
TopTowerSW <- data.frame(x = 0.5, y = 0.5, x= 40)
# Function to compute the distance between points and the line
distance_to_line <- function(cloud, TopTowerSW, theta_deg) {
# Convert angle to radians
theta_rad <- theta_deg * pi / 180
# Direction vector for the line
d <- c(0, -1, tan(theta_rad))
# Filter the cloud to retain only points to the south of TopTowerSW
cloud_south <- cloud[cloud$y < TopTowerSW$y, ]
# Convert both data.frames to matrices for easier computation
cloud_matrix <- as.matrix(cloud_south)
TopTowerSW_matrix <- matrix(rep(as.matrix(TopTowerSW), nrow(cloud_south)), ncol=3, byrow=TRUE)
# Compute the position vectors relative to TopTowerSW for all points
r <- cloud_matrix - TopTowerSW_matrix
# Compute the cross product between r and d for all points
cross_product <- cbind(r[,2]*d[3] - r[,3]*d[2], r[,3]*d[1] - r[,1]*d[3], r[,1]*d[2] - r[,2]*d[1])
# Compute the distance using the formula mentioned
distances <- sqrt(rowSums(cross_product^2)) / sqrt(sum(d^2))
return(distances)
}
# Example usage: Compute the distance with a vertical angle of 45 degrees
distances <- distance_to_line(cloud, TopTowerSW, -60)
plot3d(cloud, size = 5)
plot3d(TopTowerSW, size = 10, add = T, col = "red")
Puntitos <- as.numeric(names(distances[distances < 4]))
plot3d(cloud[Puntitos, ], size = 10, col = "blue", add = T)
compute_endpoint <- function(TopTowerSW, cloud, theta_deg) {
# Direction vector for the line based on the vertical angle
d <- c(0, -1, tan(theta_deg * pi / 180))
# Calculate the endpoint of the line based on the vertical angle
distance_to_edge <- TopTowerSW$y - min(cloud$y)
z_change_due_to_angle <- distance_to_edge * d[3]
end_point <- c(TopTowerSW$x, TopTowerSW$y - distance_to_edge, TopTowerSW$z + z_change_due_to_angle)
return(end_point)
}
# For demonstration: Compute endpoint with a vertical angle of 30 degrees
theta_deg <- -60
end_point <- compute_endpoint(TopTowerSW, cloud, theta_deg)
segments3d(rbind(TopTowerSW, end_point), col = 'red')
## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
## Create an artificial tower. Just the corners
CornA <- data.table(X = 3, Y = 3, Z = 40)
## CornA <- Data[identify3d(Data, n = 1), ]
## spheres3d(CornA, radius = 1, col = "yellow")
## Left click on a ground point in the point cloud
## Right click to scape
## CornB <- Data[identify3d(Data, n = 1), ]
spheres3d(CornB, radius = 1, col = "red")
## rgl.pop()
Data[, Dista := dist3d(X, Y, Z,
CornA$X, CornA$Y, CornA$Z,
CornB$X, CornB$Y, CornB$Z),
by = seq_len(nrow(Data))]
setnames(Data, "Dista", "CO")
## plot3d(Data, aspect = "iso", size = 0.5, col = "grey21")
## plot3d(Data[CO < 5 , ], aspect = "iso", size = 4, col = "red", add = T)
## spheres3d(CornA, radius = 1, col = "yellow")
## spheres3d(CornB, radius = 1, col = "red")
## pan3d(3)
Data[, H := dist3d2(X, Y, Z,
CornA$X, CornA$Y, CornA$Z),
by = seq_len(nrow(Data))]
Data[, Angle := CO/H]
Data[, Angle := (asin(Angle))*180/pi]
OpeningAngl <- 20
plot3d(Data[Angle < OpeningAngl, ], aspect = "iso", size = 1.5, col = "red")
spheres3d(CornA, radius = 1, col = "yellow")
spheres3d(CornB, radius = 1, col = "red")
plot3d(Data[Angle > OpeningAngl, ], aspect = "iso", size = 0.5, col = "grey21")
spheres3d(CornA, radius = 1, col = "yellow")
spheres3d(CornB, radius = 1, col = "red")
plot3d(Data[Angle < OpeningAngl, ], aspect = "iso", size = 0.5, col = "red", add = T)
pan3d(3)
## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
This diff is collapsed.
Click to expand it.
Preview
0%
Loading
Try again
or
attach a new file
.
Cancel
You are about to add
0
people
to the discussion. Proceed with caution.
Finish editing this message first!
Save comment
Cancel
Please
register
or
sign in
to comment