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]]
Good to hear you've made such progress. Just a couple of comments: - You should use points3d() rather than rgl.points(). The latter is a low level function that may have unpleasant side effects, especially mixing it with other *3d() functions like persp3d(). - There are several ways to draw a flat surface to illustrate your data. Which one to use really depends on the form of data. With randomly distributed points, yours is as good as any. If the values are a known function of location, there are probably better ones. Duncan Murdoch On 23/10/2020 12:18 p.m., Balint Radics wrote:> 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]] > > ______________________________________________ > R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. >
One addition to this thread: I've just committed contourLines3d() to rgl (in version 0.102.28). It will let you display contour lines of a function on any surface in a scene. Duncan Murdoch On 23/10/2020 2:30 p.m., Duncan Murdoch wrote:> Good to hear you've made such progress. Just a couple of comments: > > - You should use points3d() rather than rgl.points(). The latter is a > low level function that may have unpleasant side effects, especially > mixing it with other *3d() functions like persp3d(). > - There are several ways to draw a flat surface to illustrate your data. > Which one to use really depends on the form of data. With randomly > distributed points, yours is as good as any. If the values are a known > function of location, there are probably better ones. > > Duncan Murdoch > > > On 23/10/2020 12:18 p.m., Balint Radics wrote: >> 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]] >> >> ______________________________________________ >> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see >> https://stat.ethz.ch/mailman/listinfo/r-help >> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html >> and provide commented, minimal, self-contained, reproducible code. >> >