1

I have a RasterLayer object z:

class       : RasterLayer 

dimensions  : 100, 100, 10000  (nrow, ncol, ncell)

resolution  : 0.06534037, 0.03165736  (x, y)

extent      : -96.63908, -90.10505, 40.36051, 43.52624  (xmin, xmax, ymin, ymax)

coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 

data source : in memory

names       : layer 

values      : 0.001829022, 2.227507  (min, max)

The latitude and longitude coordinates shown with extent(z) represents the state of Iowa in the U.S. The function plot(z) produces a plot with the x-y axes consistent with the coordinates represented with extent(z), as well as the coordinate reference with projection(z) shown above.

However, the function plot3D(z) produces a plot with the x-y axes inconsistent with the coordinates represented with extent(z), and the coordinate reference with projection(z) shown above. The function plot3D(z) is apparently using a different coordinate reference system than what is shown above in my RasterLayer object z.

Question: I would like the axes on the 3D plot (using the function plot3D) to continue using the latitude and longitude shown with extent(z), which is also consistent with the reference coordinate system projection(z) shown above.

A reproducible example for illustration:

library(maps)

library(maptools)

library(akima)

library(raster)

library(rasterVis)

library(rgl)

x <- runif(10, -95, -92)

y <- runif(10, 41, 43)

intensity <- runif(10, 0, 2)

iowa <- map("state", "iowa", fill = TRUE, col = "transparent", plot = TRUE)

iowa.spat <- map2SpatialPolygons(iowa, IDs = c("iowa"))

interp <- interp.new(x, y, intensity,

                     xo = seq(min(iowa$x, na.rm = TRUE),

                              max(iowa$x, na.rm = TRUE), length = 100),

                     yo = seq(min(iowa$y, na.rm = TRUE),

                              max(iowa$y, na.rm = TRUE), length = 100),

                     extrap = TRUE)

interp.raster <- raster(interp)

z <- mask(interp.raster, iowa.spat)

z

plot(z)

rgl::open3d()

plot3D(z)

decorate3d()
  • Screenshot? Reproducible example? It is not just an aspect ratio thing is it? – mdsumner Apr 29 '16 at 4:31
  • In response to your timely comment, I modified my question to include a reproducible example. Look forward to hearing your thoughts. – user72225 Apr 29 '16 at 15:51
2

It's because plot3D is just scaling the z values to deal with the radically different scales in plotting "heights" against longitude/latitude values.

It's better to scale the aspect ratio since you don't then have to fake anything.

Try this, see how everything is perfect with adjust = FALSE - it's just that degrees and metres don't belong in the same space.

plot3D(z, adjust = FALSE)
decorate3d()

Set the aspect ratio (rather than munge the data):

aspect3d(1, 1, 1e-1)

fwiw, rgl has a lot more power than some of these high level plot wrappers expose - if you are really keen, try installing quadmesh.

## note, this did use "gris", which was never generally available
library(quadmesh) ## install.packages("quadmesh")
## create quad mesh
qm <- quadmesh(z, z)
scl <- function(x) (x - min(x, na.rm = TRUE))/diff(range(x, na.rm = TRUE))
shade3d(qm, col = terrain.colors(64)[scl(qm$vb[3,]) * 64 + 1][qm$ib])
aspect3d(1, 1, 1e-1)

Better still, because we have a quad mesh rather than an inflexible affine raster, we can reproject the x/y vertices and repeat without any aspect ratio ad hockery:

library(rgdal)
## work transpose
qm$vb[1:2, ] <- t(project(t(qm$vb[1:2, ]), "+proj=laea +ellps=WGS84 +lon_0=-93 +lat_0=42"))
open3d()
shade3d(qm, col = terrain.colors(64)[scl(qm$vb[3,]) * 64 + 1][qm$ib])

This is not a terribly exciting data set, but might be worth trying with your real data.

  • Thank you for pointing that out, and in particular for the helpful context. That's so important for my edification. I was so sure that the issue pertained to extent() and projection() that I didn't allow myself to step back to see the forest from the trees. In the end, I decided to move forward with your gris package, implementing the code that you had graciously provided on my data set. I am very happy with the results. Admittedly, I do like the name of the R package too. A thousand thanks. – user72225 Apr 29 '16 at 21:42
  • Just a note, I updated the answer to use CRAN's 'quadmesh' rather than the old 'gris' that was only on Github. – mdsumner Jul 18 '17 at 9:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.