[R] 3d plot of earth with cut

[Duncan Murdoch]

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