Version 0.4.0: 3D mapping update

-Added plot_3d, a function to take hillshades/textures and visualize them with an elevation matrix in 3D. Includes water visualization and the ability to view the map as a solid object.

-Added zscale argument to detect_water.

-Changed default orientation to match the orientation of the matrix.

-Changed detect_water and calculate_normals to not show a progress bar by default.

-Changed ambient_shade to default to 1-45 degrees, instead of 0-45.

-Changed write_png to stop when user doesn’t provide a filename.

-Updated README

DESCRIPTION

@@ -5,7 +5,7 @@ Version: 0.4.0
 Date: 2018-06-30
 Author: Tyler Morgan-Wall
 Maintainer: Tyler Morgan-Wall <[email protected]>
-Description: Uses ray tracing to produce global illumination maps of elevation matrices.
+Description: Uses a combination of raytracing, spherical UV mapping, Lambertian reflectance, and ambient occlusion to produce hillshades of elevation matrices. Includes water detection and layering, programmable color palette generation, several built-in textures, and 2D and 3D plotting options.
 License: GPL-3
 Imports: doParallel,
          foreach,
@@ -17,7 +17,8 @@ Imports: doParallel,
          png,
          akima,
          magrittr,
-         rgl
+         rgl,
+         imager
 LinkingTo: Rcpp, progress
 RoxygenNote: 6.0.1
 URL: https://github.com/tylermorganwall/rayshader

NAMESPACE

@@ -14,6 +14,7 @@ export(sphere_shade)
 export(write_png)
 import(doParallel)
 import(foreach)
+import(imager)
 import(parallel)
 import(progress)
 import(rgl)

R/ambient_shade.R

@@ -2,7 +2,7 @@
 #'
 #'@description Calculates Ambient Occlusion Shadow Map
 #'
-#'@param heightmap  two-dimensional matrix, where each entry in the matrix is the elevation at that point. All points are assumed to be evenly spaced.
+#'@param heightmap  A two-dimensional matrix, where each entry in the matrix is the elevation at that point. All points are assumed to be evenly spaced.
 #'@param anglebreaks The angle(s), in degrees, as measured from the horizon from which the light originates.
 #'@param sunbreaks Default 12. Number of rays to be sent out in a circle, evenly spaced, around the point being tested.
 #'@param maxsearch Default 20. The maximum distance that the system should propogate rays to check for .
@@ -21,7 +21,7 @@
 #'amb = ambient_shade(heightmap = volcano, 
 #'    sunbreaks = 15, 
 #'    maxsearch = 100)
-ambient_shade = function(heightmap, anglebreaks = seq(0,45,15), sunbreaks = 12, 
+ambient_shade = function(heightmap, anglebreaks = seq(1,46,15), sunbreaks = 12, 
                         maxsearch=20, multicore=FALSE, zscale=1, remove_edges=TRUE, cache_mask = NULL, 
                         shadow_cache=NULL, progbar=TRUE, ...) {
   if(sunbreaks < 3) {

R/calculate_normal.R

@@ -5,14 +5,16 @@
 #'@param heightmap A two-dimensional matrix, where each entry in the matrix is the elevation at that point. All points are assumed to be evenly spaced.
 #'@param zscale Default 1.
 #'@param remove_edges Default TRUE. Slices off artifacts on the edge of the shadow matrix.
-#'@param progbar Default `TRUE`. If `FALSE`, turns off progress bar.
+#'@param progbar Default `FALSE`. If `TRUE`, turns on progress bar.
 #'@return Matrix of light intensities at each point.
 #'@export
 #'@examples
 #'#Here we produce a light intensity map of the `volcano` elevation map.
-calculate_normal = function(heightmap, zscale=1, remove_edges = FALSE, progbar=TRUE) {
+calculate_normal = function(heightmap, zscale=1, remove_edges = FALSE, progbar=FALSE) {
+  flipud = function(x) {
+    x[,ncol(x):1]
+  }
   heightmap = heightmap / zscale
-  heightmap = heightmap[1:nrow(heightmap),ncol(heightmap):1]
   matrices = calculate_normal_cpp(heightmap=heightmap, progbar = progbar)
   returnnormal = list()
   if(remove_edges) {

R/detect_water.R

@@ -3,11 +3,14 @@
 #'@description Detects bodies of water (of a user-defined minimum size) within an elevation matrix.
 #'
 #'@param heightmap A two-dimensional matrix, where each entry in the matrix is the elevation at that point. All grid points are assumed to be evenly spaced.
+#'@param zscale Default `1`. The ratio between the x and y spacing (which are assumed to be equal) and the z axis. For example, if the elevation levels are in units
+#'of 1 meter and the grid values are separated by 10 meters, `zscale` would be 10.
 #'@param cutoff Default `0.999`. The lower limit of the z-component of the unit normal vector to be classified as water.
 #'@param min_area Default length(heightmap)/400. Minimum area (in units of the height matrix x and y spacing) to be considered a body of water.
 #'@param normalvectors Default `NULL`. Pre-computed array of normal vectors from the `calculate_normal` function. Supplying this will speed up water detection.
 #'@param remove_edges Default `TRUE`. Slices off artifacts on the edge of the shadow matrix.
 #'@param keep_groups Default `FALSE`. If `TRUE`, the matrix returned will retain the numbered grouping information.
+#'@param progbar Default `FALSE`. If `TRUE`, turns on progress bar.
 #'@return Matrix indicating whether water was detected at that point. 1 indicates water, 0 indicates no water.
 #'@export
 #'@examples
@@ -21,8 +24,8 @@
 #'  sphere_shade(texture="imhof3") %>%
 #'  add_water(detect_water(island_volcano, min_area = 400),color="imhof3") %>%
 #'  plot_map()
-detect_water = function(heightmap, cutoff = 0.999, min_area=length(heightmap)/400 ,normalvectors=NULL,
-                        remove_edges=TRUE, keep_groups=FALSE) {
+detect_water = function(heightmap, zscale = 1, cutoff = 0.999, min_area=length(heightmap)/400 ,normalvectors=NULL,
+                        remove_edges=TRUE, keep_groups=FALSE, progbar = TRUE) {
   if(!is.null(normalvectors)) {
     zmatrix = abs(normalvectors$z)
     zmatrix = abs(zmatrix)
@@ -33,7 +36,7 @@ detect_water = function(heightmap, cutoff = 0.999, min_area=length(heightmap)/40
     zmatrix[nrow(zmatrix),] = 0
     zmatrix[,ncol(zmatrix)] = 0
   } else {
-    zmatrix = calculate_normal(heightmap)$z
+    zmatrix = calculate_normal(heightmap,zscale=zscale,progbar=progbar)$z
     zmatrix = abs(zmatrix)
     zmatrix[zmatrix < cutoff] = 0
     zmatrix[zmatrix >= cutoff] = 1

R/lamb_shade.R

@@ -22,7 +22,10 @@ lamb_shade = function(heightmap, rayangle=45, sunangle=315, zscale = 1, zero_neg
   sunang_rad = (-sunangle)*pi/180;
   rayang_rad = rayangle*pi/180;
   rayvector = c(sin(sunang_rad)*cos(rayang_rad),cos(sunang_rad)*cos(rayang_rad),-sin(rayang_rad))
-  heightmap = t(heightmap) / zscale;
+  flipud = function(x) {
+    x[,ncol(x):1]
+  }
+  heightmap = t(flipud(heightmap)) / zscale;
   shadowmatrix = lambshade_cpp(heightmap = heightmap, rayvector = rayvector)
   shadowmatrix = scales::rescale_max(shadowmatrix,c(0,1))
   if(zero_negative) {

R/make_base.R

@@ -0,0 +1,47 @@
+#'@title make_base
+#'
+#'@description Makes the base below the 3D elevation map.
+#'
+#'@param heightmap A two-dimensional matrix, where each entry in the matrix is the elevation at that point. All points are assumed to be evenly spaced.
+#'@param basedepth Default `0`.
+#'@param basecolor Default `grey20`.
+#'@param zscale Default `1`. The ratio between the x and y spacing (which are assumed to be equal) and the z axis. For example, if the elevation levels are in units
+#'of 1 meter and the grid values are separated by 10 meters, `zscale` would be 10.
+#'@keywords internal
+make_base = function(heightmap,basedepth=0,basecolor="grey20",zscale=1) {
+  heightmap = heightmap[,ncol(heightmap):1]/zscale
+  heightmap1 = heightmap[1,]
+  heightmap2 = heightmap[,1]
+  heightmap3 = heightmap[nrow(heightmap),]
+  heightmap4 = heightmap[,ncol(heightmap)]
+  heightlist1 = list()
+  heightlist2 = list()
+  heightlist3 = list()
+  heightlist4 = list()
+  heightlist5 = list()
+  for(i in 1:(length(heightmap1)-1)) {
+    heightlist1[[i]] = matrix(c(1,1,1, heightmap1[i],basedepth,basedepth, i,i,i+1),3,3)
+    heightlist1[[i+length(heightmap1)]] = matrix(c(1,1,1, heightmap1[i],basedepth,heightmap1[i+1], i,i+1,i+1),3,3)
+  }
+  heightmat1 = do.call(rbind,heightlist1)
+  for(i in 1:(length(heightmap2)-1)) {
+    heightlist2[[i]] = matrix(c(i,i+1,i,  heightmap2[i],basedepth,basedepth, 1,1,1),3,3)
+    heightlist2[[i+length(heightmap2)]] = matrix(c(i,i+1,i+1,  heightmap2[i],heightmap2[i+1],basedepth, 1,1,1),3,3)
+  }
+  heightmat2 = do.call(rbind,heightlist2)
+  for(i in 1:(length(heightmap3)-1)) {
+    heightlist3[[i]] = matrix(c(nrow(heightmap),nrow(heightmap),nrow(heightmap),  heightmap3[i],basedepth,basedepth, i,i+1,i),3,3)
+    heightlist3[[i+length(heightmap3)]] = matrix(c(nrow(heightmap),nrow(heightmap),nrow(heightmap),  heightmap3[i],heightmap3[i+1],basedepth, i,i+1,i+1),3,3)
+  }
+  heightmat3 = do.call(rbind,heightlist3)
+  for(i in 1:(length(heightmap4)-1)) {
+    heightlist4[[i]] = matrix(c(i,i,i+1, heightmap4[i],basedepth,basedepth, ncol(heightmap),ncol(heightmap),ncol(heightmap)),3,3)
+    heightlist4[[i+length(heightmap4)]] = matrix(c(i,i+1,i+1, heightmap4[i],basedepth,heightmap4[i+1], ncol(heightmap),ncol(heightmap),ncol(heightmap)),3,3)
+  }
+  heightmat4 = do.call(rbind,heightlist4)
+  heightlist5[[1]] =  matrix(c(0,nrow(heightmap),nrow(heightmap),basedepth,basedepth,basedepth,0,0,ncol(heightmap)),3,3)
+  heightlist5[[2]] =   matrix(c(0,0,nrow(heightmap),basedepth,basedepth,basedepth,ncol(heightmap),0,ncol(heightmap)),3,3)
+  heightmat5 = do.call(rbind,heightlist5)
+  fullsides = rbind(heightmat1,heightmat2,heightmat3,heightmat4,heightmat5)
+  rgl::triangles3d(fullsides,lit=FALSE,color=basecolor,front="fill",back="culled")
+}

R/make_lines.R

@@ -0,0 +1,27 @@
+#'@title make_lines
+#'
+#'@description Makes the lines in the corner of the base.
+#'
+#'@param heightmap A two-dimensional matrix, where each entry in the matrix is the elevation at that point. All points are assumed to be evenly spaced.
+#'@param basedepth Default `0`.
+#'@param linecolor Default `grey40`. 
+#'@param zscale Default `1`. The ratio between the x and y spacing (which are assumed to be equal) and the z axis. For example, if the elevation levels are in units
+#'of 1 meter and the grid values are separated by 10 meters, `zscale` would be 10.
+#'@param alpha Default `1`. Transparency.
+#'@param linewidth Default `2`. Linewidth
+#'@keywords internal
+make_lines = function(heightmap,basedepth=0,linecolor="grey40",zscale=1,alpha=1,linewidth = 2) {
+  heightmap = heightmap[,ncol(heightmap):1]/zscale
+  heightval1 = heightmap[2,2]
+  heightval2 = heightmap[nrow(heightmap)-1,2]
+  heightval3 = heightmap[2,ncol(heightmap)-1]
+  heightval4 = heightmap[nrow(heightmap)-1,ncol(heightmap)-1]
+  heightlist = list()
+  heightlist[[1]] = matrix(c(1,1,basedepth,heightval1,1,1),2,3)
+  heightlist[[2]] = matrix(c(nrow(heightmap),nrow(heightmap),basedepth,heightval2,1,1),2,3)
+  heightlist[[3]] = matrix(c(1,1,basedepth,heightval3,ncol(heightmap),ncol(heightmap)),2,3)
+  heightlist[[4]] = matrix(c(nrow(heightmap),nrow(heightmap),basedepth,heightval4,ncol(heightmap),ncol(heightmap)),2,3)
+  for(i in 1:4) {
+    rgl::lines3d(heightlist[[i]],color=linecolor,lwd=linewidth,alpha=alpha)
+  }
+}

R/make_shadow.R

@@ -0,0 +1,20 @@
+#'@title make_shadow
+#'
+#'@description Makes the base below the 3D elevation map.
+#'
+#'@param rows A two-dimensional matrix, where each entry in the matrix is the elevation at that point. All points are assumed to be evenly spaced.
+#'@param cols Default `0`.
+#'@param basedepth Default `grey20`.
+#'@param shadowwidth Default `50`. Shadow width.
+#'@import imager
+#'@keywords internal
+make_shadow = function(rows, cols, basedepth, shadowwidth) {
+  basedepth = matrix(basedepth,nrow = rows+shadowwidth*2, ncol = cols+shadowwidth*2)
+  imagemat = matrix(1,nrow = rows+shadowwidth*2, ncol = cols+shadowwidth*2)
+  imagemat[shadowwidth:(rows-shadowwidth),shadowwidth:(cols-shadowwidth)] = 0
+  shadowarray = as.array(imager::isoblur(imager::as.cimg(imagemat),sigma=shadowwidth/2))[,,1,1]
+  browser()
+  tempmap = tempfile()
+  write_png(shadowarray,tempmap)
+  rgl.surface((-shadowwidth+1):(rows+shadowwidth),(-shadowwidth):(cols+shadowwidth-1),basedepth,texture=paste0(tempmap,".png"),lit=FALSE)
+}

R/make_water.R

@@ -0,0 +1,52 @@
+#'@title make_water
+#'
+#'@description Makes the base below the 3D elevation map.
+#'
+#'@param heightmap A two-dimensional matrix, where each entry in the matrix is the elevation at that point. All points are assumed to be evenly spaced.
+#'@param waterheight Default `0`.
+#'@param watercolor Default `blue`.
+#'@param zscale Default `1`. The ratio between the x and y spacing (which are assumed to be equal) and the z axis. For example, if the elevation levels are in units
+#'of 1 meter and the grid values are separated by 10 meters, `zscale` would be 10.
+#'@param wateralpha Default `0.5`. Water transparency.
+make_water = function(heightmap,waterheight=mean(heightmap),watercolor="lightblue",zscale=1,wateralpha=0.5) {
+  heightmap = heightmap[,ncol(heightmap):1]/zscale
+  heightmap1 = heightmap[1,]
+  heightmap2 = heightmap[,1]
+  heightmap3 = heightmap[nrow(heightmap),]
+  heightmap4 = heightmap[,ncol(heightmap)]
+  heightlist1 = list()
+  heightlist2 = list()
+  heightlist3 = list()
+  heightlist4 = list()
+  heightlist5 = list()
+  for(i in 1:(length(heightmap1)-1)) {
+    if(heightmap1[i] < waterheight) {
+      heightlist1[[i]] = matrix(c(1,1,1, heightmap1[i],waterheight,waterheight, i,i,i+1),3,3)[c(1,3,2),]
+      heightlist1[[i+length(heightmap1)]] = matrix(c(1,1,1, heightmap1[i],waterheight,heightmap1[i+1], i,i+1,i+1),3,3)[c(1,3,2),]
+    }
+  }
+  heightmat1 = do.call(rbind,heightlist1)
+  for(i in 1:(length(heightmap2)-1)) {
+    if(heightmap2[i] < waterheight) {
+      heightlist2[[i]] = matrix(c(i,i+1,i,  heightmap2[i],waterheight,waterheight, 1,1,1),3,3)[c(1,3,2),]
+      heightlist2[[i+length(heightmap2)]] = matrix(c(i,i+1,i+1,  heightmap2[i],heightmap2[i+1],waterheight, 1,1,1),3,3)[c(1,3,2),]
+    }
+  }
+  heightmat2 = do.call(rbind,heightlist2)
+  for(i in 1:(length(heightmap3)-1)) {
+    if(heightmap3[i] < waterheight) {
+      heightlist3[[i]] = matrix(c(nrow(heightmap),nrow(heightmap),nrow(heightmap),  heightmap3[i],waterheight,waterheight, i,i+1,i),3,3)[c(1,3,2),]
+      heightlist3[[i+length(heightmap3)]] = matrix(c(nrow(heightmap),nrow(heightmap),nrow(heightmap),  heightmap3[i],heightmap3[i+1],waterheight, i,i+1,i+1),3,3)[c(1,3,2),]
+    }
+  }
+  heightmat3 = do.call(rbind,heightlist3)
+  for(i in 1:(length(heightmap4)-1)) {
+    if(heightmap4[i] < waterheight) {
+      heightlist4[[i]] = matrix(c(i,i,i+1, heightmap4[i],waterheight,waterheight, ncol(heightmap),ncol(heightmap),ncol(heightmap)),3,3)[c(1,3,2),]
+      heightlist4[[i+length(heightmap4)]] = matrix(c(i,i+1,i+1, heightmap4[i],waterheight,heightmap4[i+1], ncol(heightmap),ncol(heightmap),ncol(heightmap)),3,3)[c(1,3,2),]
+    }
+  }
+  heightmat4 = do.call(rbind,heightlist4)
+  fullsides = rbind(heightmat1,heightmat2,heightmat3,heightmat4)
+  rgl::triangles3d(fullsides,lit=FALSE,color=watercolor,alpha=wateralpha,front="fill",back="culled")
+}

R/make_waterlines.R

@@ -0,0 +1,76 @@
+#'@title make_waterlines
+#'
+#'@description Makes the edge lines of 
+#'
+#'@param heightmap A two-dimensional matrix, where each entry in the matrix is the elevation at that point. All points are assumed to be evenly spaced.
+#'@param waterdepth Default `0`.
+#'@param linecolor Default `grey40`. 
+#'@param zscale Default `1`. The ratio between the x and y spacing (which are assumed to be equal) and the z axis. For example, if the elevation levels are in units
+#'of 1 meter and the grid values are separated by 10 meters, `zscale` would be 10.
+#'@param alpha Default `1`. Transparency of lines.
+#'@param lwd Default `2`. Water line width.
+#'@keywords internal
+make_waterlines = function(heightmap,waterdepth=0,linecolor="grey40",zscale=1,alpha=1,lwd=2) {
+  heightmap = heightmap[,ncol(heightmap):1]/zscale
+  heightval1 = heightmap[1,]
+  heightval2 = heightmap[,1]
+  heightval3 = heightmap[nrow(heightmap),]
+  heightval4 = heightmap[,ncol(heightmap)]
+  heightlist = list()
+  counter = 1
+  drawing = FALSE
+  startcoord = 1
+  for(i in 1:length(heightval1)) {
+    if(heightval1[i] < waterdepth && !drawing) { 
+      startcoord = i 
+      drawing = TRUE
+    } 
+    if((heightval1[i] > waterdepth || i == length(heightval1)) && drawing) {
+      drawing = FALSE
+      heightlist[[counter]] = matrix(c(1,1,waterdepth,waterdepth,startcoord,i),2,3)
+      counter = counter + 1
+    }
+  }
+  drawing = FALSE
+  startcoord = 1
+  for(i in 1:length(heightval2)) {
+    if(heightval2[i] < waterdepth && !drawing) { 
+      startcoord = i 
+      drawing = TRUE
+    } 
+    if((heightval2[i] > waterdepth || i == length(heightval2)) && drawing) {
+      drawing = FALSE
+      heightlist[[counter]] = matrix(c(startcoord,i,waterdepth,waterdepth,1,1),2,3)
+      counter = counter + 1
+    }
+  }
+  drawing = FALSE
+  startcoord = 1
+  for(i in 1:length(heightval3)) {
+    if(heightval3[i] < waterdepth && !drawing) { 
+      startcoord = i 
+      drawing = TRUE
+    } 
+    if((heightval3[i] > waterdepth || i == length(heightval3)) && drawing) {
+      drawing = FALSE
+      heightlist[[counter]] = matrix(c(nrow(heightmap),nrow(heightmap),waterdepth,waterdepth,startcoord,i),2,3)  
+      counter = counter + 1
+    }
+  }
+  drawing = FALSE
+  startcoord = 1
+  for(i in 1:length(heightval4)) {
+    if(heightval4[i] < waterdepth && !drawing) { 
+      startcoord = i 
+      drawing = TRUE
+    } 
+    if((heightval4[i] > waterdepth || i == length(heightval4)) && drawing) {
+      drawing = FALSE
+      heightlist[[counter]] = matrix(c(startcoord,i,waterdepth,waterdepth,ncol(heightmap),ncol(heightmap)),2,3)
+      counter = counter + 1
+    }
+  }
+  for(i in 1:length(heightlist)) {
+    rgl::lines3d(heightlist[[i]],color=linecolor,lwd=lwd,alpha=alpha,depth_mask=FALSE)
+  }
+}