Adding point locations to a 3D DEM plot in R

Question

I have some point locations which include UTMs and Elevation as a data frame I also have a DEM layer.

I have figured out how to plot the DEM in 3D using plot3D in rgl.

I can also plot the points in 3D using points3d.

I have been able to put them in the same plot using points3d with add=TRUE however the points and DEM are radically far away from each other.

In the code below I also tried to change this to a spatial data frame but rgl doesn't seem to like that.

Is it possible to plot them together with the points laying over the DEM?

I have searched and searched for a solution to this.

Here is the R code I have used so far:

> library(raster)
> library(rgdal)
> library(maptools)
> library(rgeos)
> library(lattice)
> library(latticeExtra)
> library(sp)
> library(rasterVis)
> library(rgl)
> 
> # taking data read from a .csv of UTM and elevation values
> 
> Points.Sp <- data.frame(Points=Rawdata$PointName, UTM.N=Rawdata$UTM.N, UTM.W=Rawdata$UTM.W, Elevation=Rawdata$Elevation)
> Points.Sp <- unique(Points.Sp) #weeding out duplicates
> Points.Sp <- Points.Sp[,c(3,2,4)] #getting rid of point names # I realize this looks messy but it gets what I want
> head(Points.Sp)
    UTM.W   UTM.N Elevation
1  275815 3879223      1340
8  274813 3879727      1325
29 275312 3879727      1258
45 275812 3879724      1169
66 276313 3879727      1067
75 276813 3879727      1208
> 
> dem.in <- raster("D:/Thesis/SouthernApps/Coweeta/Coweeta/DEM_30m_wgs84.img") # reading in DEM
> plot(dem.in) # check in 2D # takes a long time very large, need to crop
> 
> dem.crop <- crop(dem.in, c(272000, 282000, 3878000, 3884000))
> plot(dem.crop) # check in 2D, looks good.
> 
> plot3D(dem.crop) # plot in 3D looks like exactly what I want
> 
> points3d(Points.Sp, pch=19, cex=2, col="black", add=TRUE) # adds the points to plot but in wrong place
> 
> #attempting to set a CRS in case this is the problem.
> coordinates(Points.Sp)=c(1,2)
> proj4string(Points.Sp)=CRS("++proj=utm +zone=17") # set CRS
> str(Points.Sp)
Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
  ..@ data       :'data.frame': 71 obs. of  1 variable:
  .. ..$ Elevation: int [1:71] 1340 1325 1258 1169 1067 1208 1256 1089 1031 959 ...
  ..@ coords.nrs : num [1:2] 1 2
  ..@ coords     : num [1:71, 1:2] 275815 274813 275312 275812 276313 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:71] "1" "8" "29" "45" ...
  .. .. ..$ : chr [1:2] "UTM.W" "UTM.N"
  ..@ bbox       : num [1:2, 1:2] 274309 3878440 279876 3883732
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:2] "UTM.W" "UTM.N"
  .. .. ..$ : chr [1:2] "min" "max"
  ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
  .. .. ..@ projargs: chr "+proj=utm +zone=17 +ellps=WGS84"
> 
> # trying this a different way after setting CRS
> x <- Points.Sp@coords[1:71,1]
> y <- Points.Sp@coords[1:71,2]
> z <- Points.Sp@data$Elevation
> m <- data.frame(x=x,y=y,z=z)
> 
> plot3D(dem.crop) #again, plot in 3D looks like exactly what I want
> points3d(m, pch=19, cex=2, col="black", add=TRUE) # still adds the points to plot but in wrong place

This code reproduces the problem.

## define a Raster object
data(volcano)
r <- raster(volcano)
extent(r) <- c(0, 610, 0, 870)

## extract sample points
xy <- sampleRandom(r1, 100, xy = TRUE)     
r1<-data.frame(x=seq(0, 500, length=(71)), y=seq(0, 500, length=(71)), z=seq(0,500, length=(71)))

## display them
plot3D(r, adjust = FALSE)

points3d(r1, add=TRUE)

This looks like a bug or design flaw in `plot3D`, which isn't an `rgl` function. I'm not sure which package you got it from. The `rgl` function is `plot3d`, which likely won't work on your `dem.crop` object. — user2554330, Aug 13 '16 at 10:58
Please, post a [reproducible question](http://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example). — Oscar Perpiñán, Aug 13 '16 at 11:30
I would like the points to lie flat on the 3d raster...I just found a solution to this which I will post below as soon as I find a way to reproduce it with the volcano example — CAWA, Aug 14 '16 at 16:25
Sorary, I don't have a proper understanding of your problem. — cuttlefish44, Aug 14 '16 at 23:20

score 2 · Answer 1 · answered Aug 13 '16 at 11:27

2

As documented in the help page, both the x-axis and y-axis are adjusted with the z values. You can disable this default setting with adjust = FALSE:

library(rgl)
library(rasterVis)

## define a Raster object
data(volcano)
r <- raster(volcano)
extent(r) <- c(0, 610, 0, 870)

## extract sample points
xy <- sampleRandom(r, 100, xy = TRUE)     

## display them
plot3D(r, adjust = FALSE)
points3d(xy)

answered Aug 13 '16 at 11:27

Oscar Perpiñán

4,491
17
28

Suggestion for `rasterVis`: define S3 methods for `rgl::plot3d`, rather than a new incompatible function `plot3D`. Use the `rgl::aspect3d` function to adjust the display rather than modifying the data. – user2554330 Aug 13 '16 at 11:36
I just tried the posted solution of adjust=FALSE and, while this sort of helped to align the points with the DEM it did not set them over the landscape, rather the points were still very far away and even when adjusted in the viewing window, did not line up all that well. – CAWA Aug 13 '16 at 23:31
Excuse me, I should say that your example ran fine, however I could not get my data to do the same thing – CAWA Aug 14 '16 at 00:08
Once again, post a reproducible question. – Oscar Perpiñán Aug 14 '16 at 09:19
@user2554330 thanks for the suggestions. I will consider them with the original developer of this function. – Oscar Perpiñán Aug 14 '16 at 09:22

CAWA · Accepted Answer · 2016-08-14T19:35:29.280

## define a Raster object
data(volcano)
r <- raster(volcano)
extent(r) <- c(0, 610, 0, 870)

## extract sample points
xy <- sampleRandom(r1, 100, xy = TRUE)     

#must extract the data from the raster and recombine with the xy data.
#I don't know why this is different than simply using the raw values but it
#provides the desired effect. 
  
r1<-data.frame(x=seq(0, 500, length=(71)), y=seq(0, 500, length=(71)))
z<-extract(r, r1)
r1$z<-z
## display them
plot3D(r, adjust = FALSE)

points3d(r1, add=TRUE)
  
#points now lie flat on 3d image.

Points flush to 3d Image

Image for original problem

Adding point locations to a 3D DEM plot in R

2 Answers2