Strange 3D output from rasterVis::plot3D when plotting a digital elevation model

Question

I'm trying to plot a few 3D digital elevation models using rasterVis::plot3D.

The DEMs are originally geo tifs, they are all DEMs of the same location just different areas in that location and they all have the same coordinate reference system ("+proj=utm +zone=21 +south +datum=WGS84 +units=m +no_defs"), the only thing that differs is the spatial extent and resolution. One is an area of approx. 200 x 700 meters at 2 meter resolution and the rest are 1 x 1m at 1cm resolution. The x and y values are coordinates and z is elevation. They all plot just fine and look normal when plotting them in 2D.

However, when I plot the larger area 2m resolution DEM using rasterVis::plot3D I get a beautiful 3D plot:

just as expected, but when I plot any of the smaller, finer resolution DEMs I get a very strange output

.

I can't figure out what is causing this strange output.

Onedrive link to download the geo tif file

This is what Q1cm.tif looks like when plotted in 2D (the square is a plastic quadrat - we did a vegetation survey within the quadrat).

library(raster)
library(rasterVis)
library(rgl)

dem2m<-raster("lidar.tif")
dem1cm<-raster("Q1cm.tif")

dem2m
class      : RasterLayer 
dimensions : 197, 292, 57524  (nrow, ncol, ncell)
resolution : 2, 2  (x, y)
extent     : 360047, 360631, 3080519, 3080913  (xmin, xmax, ymin, ymax)
crs        : +proj=utm +zone=21 +south +datum=WGS84 +units=m +no_defs 
source     : lidar.tif 
names      : lidar 
values     : 18.7875, 59.18345  (min, max)

dem1cm
class      : RasterLayer 
dimensions : 314, 336, 105504  (nrow, ncol, ncell)
resolution : 0.01, 0.01  (x, y)
extent     : 360091.7, 360095, 3080597, 3080601  (xmin, xmax, ymin, ymax)
crs        : +proj=utm +zone=21 +south +datum=WGS84 +units=m +no_defs 
source     : memory
names      : Q1cm 
values     : 54.84027, 55.15403  (min, max)

rgl::open3d()
plot3D(dem1cm)

Any ideas on what's going on or how to get around the issue?

Thanks in advance!

Answer 1

The problem is rounding errors. Most rgl graphics uses single precision floating point, and that's not enough for your data. For example, the y values range from 10716468 to 10716479. In single precision there appear to be just 12 distinguishable values in that range (i.e. the whole numbers). That's why your plot looks so weird: whole ranges of y values are overplotted at the same location because of the rounding to single precision.

Your actual data is stored in double precision, which is fine. I think you could get a plot if you subtracted the mean of the x and y values (or the min, or max: basically anything in that range) from the ones extracted from your dataset. The subtraction needs to happen before you send the data to rgl, and then when rgl rounds it to single precision, you won't lose so much.

Unfortunately I don't know enough about raster and rasterVis to know how to do that. The raster::scale() function appears to do that for the z values, but they're not the problem here.

EDITED to add: @BenBolker provided the solution. You can modify the "extent" of the raster using the extent<- function. For example,

E <- as.vector(extent(dem1cm))
extent(dem1cm) <- E - rep(E[c(1, 3)], each = 2)
plot3D(dem1cm)
decorate3d()    # adds the scale

which produces a plot that looks like this (after some rotation):