6

I am currently trying to create an interactive 3D surface plot using plotly in R. My dataset includes a matrix of Latitude and Longitude and I have a DEM for the region I would like to make this surface plot for. However, I am unable to do the same.

I have a dataframe Points.sp that contains a list of Lat/Long along with elevation data for areas where I have sampled.

> head(Points.sp)
                         Points     Long      Lat Elevation
1                                    NA       NA        NA
2                      Nair212       NA       NA        NA
4                      CESF543 75.72130 12.88470   876.000
5                     CESF1587 75.93559 11.93745  1494.523
6          tysoni_WILD15AMP650 75.82200 12.38000  1051.000
7           tysoni_SDBDU201274 75.65570 12.22010  1176.000

Here's what I have tried:

Points_new <- data.frame(Points.sp$Long, Points.sp$Lat, Points.sp$Elevation)
Points_new <- as.matrix(Points_new)
p <- plot_ly(z = ~Points_new) %>% add_surface()
p

enter image description here

In the above plot, I manually loaded the elevation values for each of these points, but the problem being that I need these points on a DEM for the region, and everytime I try to load the DEM, I get an error saying - 'z' must be a numeric matrix.

The DEM here is a raster for the region and I obtained it from the SRTM dataset and cropped it to contain my region only.

 elev <- raster("E:\\MasterThesis\\Absence Approach\\Data\\Elevation_Data\\alt\\alt")
 cr1 <- crop(elev, extent(WG.new)) #WG.new is a mask for a smaller region
 newelev<- mask(cr1, WG.new)

 > newelev
  class       : RasterLayer 
  dimensions  : 1577, 624, 984048  (nrow, ncol, ncell)
  resolution  : 0.008333334, 0.008333334  (x, y)
  extent      : 72.87501, 78.07501, 8.11667, 21.25834  (xmin, xmax, ymin, ymax)
  coord. ref. : +proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs 
  data source : in memory
  names       : alt 
  values      : -14, 2625  (min, max)
  attributes  :
         ID COUNT
  from: -431     2
  to  : 8233     1

If I try loading just the DEM alone by converting the raster to a matrix, I am able to view the DEM, but now I need to add the Lat and Long for locations from the Points.sp dataframe.

newelev <- as.matrix(newelev)
p <- plot_ly(z = ~newelev) %>% add_surface()
p

enter image description here

Any suggestions? I would ideally want points as small circles on the DEM for the region along with the DEM smoothed and stretched out.

*****UPDATE*****

Based on @wittich's elegant solution, I tried what he posted, but I am getting an image like the screenshot posted below.

enter image description here

  • Please either supply your data or at least some summaries. We have no idea what your z is, or your Points.sp data frame. This is clearly a problem arising from the format of your data and without that we can't help. – Spacedman Jun 29 '17 at 9:35
  • @Spacedman I have included the data. Please let me know if you need any other information. – Vijay Ramesh Jun 29 '17 at 12:43
  • I don't think I can justify the time it would take me to work with a 777million cell raster - can you create a smaller raster, or give some code we can cut and paste to create a raster to illustrate the problem. – Spacedman Jun 29 '17 at 12:49
  • @Spacedman Made changes to include a smaller region of interest, along with sending you a link to the Data, if needed. – Vijay Ramesh Jun 29 '17 at 13:07
  • 1
    @VijayRamesh I updated my answer, is that what you were looking for? – wittich Jul 12 '17 at 20:41
1

The problem is here, that you have to translate your own coordinates to relative coordinates.

The way you use ploty you totally ignore the absolute coordinates. For plotting the DEM SRTM data you don't have the problem as it use the X, Y of the TIFF file.

Your newelev.tif example has the size of 624 x 1577 Pixel so your relative coordinates has to bee between

P1 <- c(0, 0)
P2 <- c(624, 1577)

In absolute values that would be:

P1.wgs84 <- c(8.1166702191853, 72.8750131881337)
P2.wgs84 <- c(21.2583375712263, 78.0750134593281)

I hope I could show you the problem. If I find some time I will look for a programmable solution for you.


Update

As promised here a way how you could do it... there is probably a better way, but anyway it works.

library("plotly")
library("rgdal")
library("raster")

# Point <- Points.sp[,c("Long", "Lat", "Elevation")] # to get ride of the first column
Point <- rbind(c(75.72130, 12.88470, 876.000),
                c(75.93559, 11.93745, 1494.523),
                c(75.82200, 12.38000, 1051.000),
                c(75.65570, 12.22010, 1176.000))
colnames(Point) <- c("X", "Y", "Z")
rownames(Point) <- c(1:4)
head(Point)

Point.extent <- extent(Point[,1:2])
Point.extent <- Point.extent + 0.5 # grow extend 0.5 degree
Point.extent

setwd("/home/wittich/Downloads/GISStackExchange/")
elev <- raster("newelev.tif")
elev.new1 <- crop(elev, Point.extent) # use point extent
elev.new1

d <- dim(elev.new1) # get same dimension
point.raster <- raster(nrow=d[1], ncol=d[2], extent(elev.new1))
cells <- cellFromXY(point.raster, Point[,1:2])
point.raster[cells] <- Point[,3] # set height value to cells
crs(point.raster) <- crs(elev.new1) # copy crs just in case
point.raster

elev.new2 <- as.matrix(elev.new1)

# here comes the tricky part where it gets the row and col number by using the cell number and set it to a data frame
point.raster2 <- as.matrix(point.raster)
point.cell.value <- as.data.frame(point.raster, na.rm=TRUE )
point.df <- data.frame(rowColFromCell(point.raster, as.numeric(rownames(point.cell.value))), height=point.cell.value[,1])
head(point.df)

p <- plot_ly(z = elev.new2, type = "surface") %>% 
         add_trace(data = point.df, x = ~col, y = ~row, z = ~height, mode = "markers", type = "scatter3d", 
             marker = list(size = 5, color = "red", symbol = 104))
p

output screenshot

Hope that was what you were looking for.


Update II

I had also trouble to display the plot correctly, I thought that is a problem of my ubuntu R/rstudio setup. What helped was to open the plot in a Firefox browser window:

open in new window

  • Using the absolute coordinates would be also an option! – wittich Jul 12 '17 at 20:43
  • Thanks a ton @wittich. I shall play around with the code and try to get at a slightly different angle, smaller size of points as well as smoothen out the DEM. This is a great start! – Vijay Ramesh Jul 13 '17 at 3:58
  • Still having some trouble @wittich. Can I e-mail you my code so far? Perhaps you might know where I am going wrong. I am getting an error message, when I try setting height to cells : point.raster[cells] <- Point[,3] # set height value to cells – Vijay Ramesh Jul 13 '17 at 5:11
  • Check if Point is a matrix and has 3 columns, my using is.matrix(Point). I used rbind to create the matrix. I created my code in a way that you can run it completely... so just compare every step, and output of each variable. In your case it should work to create the matrix Point like this: Point <- Points.sp[,c("Long", "Lat", "Elevation")]. Try check for example if head(Point[,3]) (which is the 3rd column) is similar to your data. And let me know if you need further help... – wittich Jul 13 '17 at 8:20
  • So, it is a matrix.. I think I sent you an invite to a chat where I have mentioned the other error. It looks like it works, but I am unable to see anything. There is a single dot in the background that is supposedly the DEM as well as points. Let me know if you got my chat invite. – Vijay Ramesh Jul 13 '17 at 10:58
-2

'z' must be a numeric matrix. Perhaps you need to set the matrix to number.

Convert your data frame dat2 to matrix with function as.matrix()? I got the example from a stack over flow question

  • 1
    That aspect is taken care of when the elevation raster is loaded, as the DEM has a z value. Your answer does not help me with an output that has a DEM with points shown on the DEM. – Vijay Ramesh Jul 6 '17 at 4:10
  • the point was to set the value to numeric matrix. – lewis Sep 18 '17 at 19:55

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.