[R] 3d plot of earth with cut

[Duncan Murdoch]

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 using 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 using 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.
>