Dear All,
Thanks a lot for the useful help again. I manage to get it done up to a
point where I think I
just need to apply some smoothing/interpolation to get denser points, to
make it nice.
Basically, I started from Duncen's script to visualize and make the
clipping along a plane
at a slice.
Then I map my data points' values to a color palette and just plot them as
points on this
plane. Since I have already the (x,y,z) coordinates for my points in the
slice's plane
I just plot them directly. I copied the code below..
To make it nicer would be able to make a real "smooth" map on the 2D
surface, rather
than plotting all points (e.g. using polygons?).
Best,
Balint
############################################################
# Construct a plane at a given longitude
r <- 6378.1 # radius of Earth in km
fixlong <- 10.0*pi/180.0 # The longitude slice
# Construct a plane in 3D:
# Let vec(P0) = (P0x, P0y, P0z) be a point given in the plane of the
longitude
# let vec(n) = (nx, ny, nz) an orthogonal vector to this plane
# then vec(P) = (Px, Py, Pz) will be in the plane if (vec(P) - vec(P0)) *
vec(n) = 0
# We pick 2 arbitrary vectors in the plane out of 3 points
p0x <- r*cos(2)*cos(fixlong)
p0y <- r*cos(2)*sin(fixlong)
p0z <- r*sin(2)
p1x <- r*cos(2.5)*cos(fixlong)
p1y <- r*cos(2.5)*sin(fixlong)
p1z <- r*sin(2.5)
p2x <- r*cos(3)*cos(fixlong)
p2y <- r*cos(3)*sin(fixlong)
p2z <- r*sin(3)
# Make the vectors pointing to P and P0
v1x <- p1x - p0x # P
v1y <- p1y - p0y
v1z <- p1z - p0z
v2x <- p2x - p0x # P0
v2y <- p2y - p0y
v2z <- p2z - p0z
# The cross product will give a vector orthogonal to the plane, (nx, ny, nz)
nx <- v1y*v2z - v1z*v2y;
ny <- v1z*v2x - v1x*v2z;
nz <- v1x*v2y - v1y*v2x;
# normalize
nMag <- sqrt(nx*nx + ny*ny + nz*nz);
nx <- nx / nMag;
ny <- ny / nMag;
nz <- nz / nMag;
# Plane equation (vec(P) - vec(P0)) * vec(n) = 0, with P=(x, y, z), P0=(x0,
y0, z0),
# giving a*(x-x0)+b*(y-y0)+c*(z-z0) = 0, where x,x0 are two points in the
plane
# a, b, c are the normal vector coordinates
a <- -nx
b <- -ny
c <- -nz
d <- -(a*v2x + b*v2y + c*v2z )
open3d()
# Plot the globe - from Duncan
# points of a sphere
lat <- matrix(seq(90, -90, len = 50)*pi/180, 50, 50, byrow = TRUE)
long <- matrix(seq(-180, 180, len = 50)*pi/180, 50, 50)
x <- r*cos(lat)*cos(long)
y <- r*cos(lat)*sin(long)
z <- r*sin(lat)
# Plot with texture
ids <- persp3d(x, y, z, col = "white",
texture = system.file("textures/world.png", package
"rgl"),
specular = "black", axes = FALSE, box = FALSE, xlab =
"",
ylab = "", zlab = "", normal_x = x, normal_y
= y, normal_z
= z)
# Plot the plane across the longitude slice
#planes3d(a, b, c, d, alpha = 0.6) # optionally visualize the plane
# Apply clipping to only one side of the plane using the normal vector
clipplanes3d(a, b, c, d)
# Map something onto this plane - how? Let's try with rgl.points and
mapping the colors
# The data is: data_activity and variables are $X, $Y, $Z, $Ar
library(leaflet)
# map the colors to the data values
pal <- colorNumeric(
palette = "Blues",
domain = data_activity$Ar) #
# plot the points and the mapped colors
rgl.points( data_activity$X, data_activity$Y, data_activity$Z, color
pal(data_activity$Ar), size=3)
############################################################
On Fri, Oct 23, 2020 at 1:50 AM aBBy Spurdle, ?XY <spurdle.a at gmail.com>
wrote:
> > It should be a 2D slice/plane embedded into a 3D space.
>
> I was able to come up with the plot, attached.
> My intention was to plot national boundaries on the surface of a sphere.
> And put the slice inside.
> However, I haven't (as yet) worked out how to get the coordinates for
> the boundaries.
>
> Let me know, if of any value.
> And I'll post the code.
> (But needs to be polished first)
>
[[alternative HTML version deleted]]