3d plot of earth with cut

Hello,

Could someone suggest a package/way to make a 3D raster plot of the Earth
(with continent boundaries), and then make a "cut" or "slice" of it such
that one can also visualize some scalar quantity as a function of the
Radius/Depth across that given slice ?

Formally, I would have a given, fixed longitude, and a list of vectors
{latitude, radius, Value}
that would show the distribution of the quantity "Value" at various depths
and latitudes in 3D.

Thanks a lot in advance,
Balint

        [[alternative HTML version deleted]]

Hello,

Could someone suggest a package/way to make a 3D raster plot of the Earth
(with continent boundaries), and then make a "cut" or "slice" of it such
that one can also visualize some scalar quantity as a function of the
Radius/Depth across that given slice ?

Formally, I would have a given, fixed longitude, and a list of vectors
{latitude, radius, Value}
that would show the distribution of the quantity "Value" at various depths
and latitudes in 3D.

Hello,

Could someone suggest a package/way to make a 3D raster plot of the Earth
(with continent boundaries), and then make a "cut" or "slice" of it such
that one can also visualize some scalar quantity as a function of the
Radius/Depth across that given slice ?

Formally, I would have a given, fixed longitude, and a list of vectors
{latitude, radius, Value}
that would show the distribution of the quantity "Value" at various depths
and latitudes in 3D.

Thanks a lot in advance,
Balint

        [[alternative HTML version deleted]]

1. Have you looked here:
https://cran.r-project.org/web/views/Spatial.html
(I assume you have done some web searches on possible terms like "3D Earth
Data R" or whatever)

2. You might try posting on the r-sig-geo list rather than here, where
relative expertise may more likely be available.

Cheers,

Bert Gunter

"The trouble with having an open mind is that people keep coming along and
sticking things into it."
-- Opus (aka Berkeley Breathed in his "Bloom County" comic strip )


On Thu, Oct 22, 2020 at 11:41 AM Balint Radics <balintradics at gmail.com>
wrote:

Hello,

Could someone suggest a package/way to make a 3D raster plot of the Earth
(with continent boundaries), and then make a "cut" or "slice" of it such
that one can also visualize some scalar quantity as a function of the
Radius/Depth across that given slice ?

Formally, I would have a given, fixed longitude, and a list of vectors
{latitude, radius, Value}
that would show the distribution of the quantity "Value" at various depths
and latitudes in 3D.

Thanks a lot in advance,
Balint

        [[alternative HTML version deleted]]

On 21/10/2020 8:45 a.m., Balint Radics wrote:

Hello,

Could someone suggest a package/way to make a 3D raster plot of the Earth
(with continent boundaries), and then make a "cut" or "slice" of it such
that one can also visualize some scalar quantity as a function of the
Radius/Depth across that given slice ?

Formally, I would have a given, fixed longitude, and a list of vectors
{latitude, radius, Value}
that would show the distribution of the quantity "Value" at various

depths

and latitudes in 3D.

The rgl package has a full sphere of the Earth with (obsolete) political
boundaries in example(persp3d).  To cut it in half along the plane
through a given longitude (and that longitude + 180 deg), you could use
clipPlanes3d, or clipObj3d.  For example,

library(rgl)
lat <- matrix(seq(90, -90, len = 50)*pi/180, 50, 50, byrow = TRUE)
long <- matrix(seq(-180, 180, len = 50)*pi/180, 50, 50)

r <- 6378.1 # radius of Earth in km
x <- r*cos(lat)*cos(long)
y <- r*cos(lat)*sin(long)
z <- r*sin(lat)

open3d()
ids <- persp3d(x, y, z, col = "white",
         texture = system.file("textures/worldsmall.png", package = "rgl"),
         specular = "black", axes = FALSE, box = FALSE, xlab = "", ylab
= "", zlab = "",
         normal_x = x, normal_y = y, normal_z = z)

clipFn <- function(coords) {
   pmax(coords[,1], coords[,2]) # Just an example...
}
clipObj3d(ids["surface"], clipFn)


Filling in the exposed surface could be done with polygon3d(), with some
work to construct the polygon. Displaying the function across that
surface could be done in a couple of ways, either by using a texture map
(like for the map), or subdividing the polygon and setting colour by the
coordinates of each vertex.

Note that the clipObj3d function isn't on CRAN yet; there you'd have to
use clipPlanes3d.  You can get the newer version from R-forge.

Duncan Murdoch

If you have "value" as a function of latitude and radius, isn't that a
2D (not 3D) scalar field?
Which can be plotted using a regular heatmap.

If you want a curved edge where depth=0 (radius=?), that's not too
difficult to achieve.
Not quite sure what continent boundaries mean in this context, but
that could possibly be added to.

Or do you want a 2D slice superimposed within a 3D space?

And if so, does it need to be dynamic?
(i.e. Rotate the whole thing, with a trackball/mouse).

Note that the further you go down the list above, the more work is
required.


On Fri, Oct 23, 2020 at 7:41 AM Balint Radics <balintradics at gmail.com>
wrote:

Hello,

Could someone suggest a package/way to make a 3D raster plot of the Earth
(with continent boundaries), and then make a "cut" or "slice" of it such
that one can also visualize some scalar quantity as a function of the
Radius/Depth across that given slice ?

Formally, I would have a given, fixed longitude, and a list of vectors
{latitude, radius, Value}
that would show the distribution of the quantity "Value" at various

depths

and latitudes in 3D.

Thanks a lot in advance,
Balint

        [[alternative HTML version deleted]]

It should be a 2D slice/plane embedded into a 3D space.

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)

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