added a function to animate on rayshader plot3d

zappingseb · Nov 4, 2019 · 7c97e8a · 7c97e8a
1 parent aec6047
commit 7c97e8a
Show file tree

Hide file tree

Showing 15 changed files with 334 additions and 13 deletions.
diff --git a/R/get_elevdata_from_bbox.R b/R/get_elevdata_from_bbox.R
@@ -28,9 +28,7 @@ get_elevdata_from_bbox <- function(bbox, type = c("SRTM", "EUDEM"), file = NULL)
 
   elevation_matrix <- .get_elmatdata_matrix_from_raw(data_raw)
 
-  if (type=="EUDEM") {
-    elevation_matrix <- t(elevation_matrix)
-  }
+  elevation_matrix <- t(elevation_matrix)
 
   if (type == "EUDEM") {
     elevation_labels <- .get_elmat_labels_bbox_size(bbox, dim(elevation_matrix))

diff --git a/R/get_elevdata_long.R b/R/get_elevdata_long.R
@@ -8,9 +8,13 @@ get_elevdata_long <- function(elevdata) {
 
   elmat_filtered_long <- reshape2::melt(elevdata, id.vars=c("deg_elmat_lat"))
 
+  elmat_filtered_long$variable <- as.numeric(as.character(elmat_filtered_long$variable))
+  elmat_filtered_long$deg_elmat_lat <- as.numeric(as.character(elmat_filtered_long$deg_elmat_lat))
   elmat_filtered_long <- elmat_filtered_long[
-    which(elmat_filtered_long$deg_elmat_lat != min(elmat_filtered_long$deg_elmat_lat)),
+    which(elmat_filtered_long$deg_elmat_lat != min(elmat_filtered_long$variable)),
+    ]
+  elmat_filtered_long <- elmat_filtered_long[
+    which(elmat_filtered_long$deg_elmat_lat != max(elmat_filtered_long$variable)),
     ]
-  elmat_filtered_long$variable <- as.numeric(as.character(elmat_filtered_long$variable))
   return(elmat_filtered_long)
 }
diff --git a/R/get_elevdata_long_water.R b/R/get_elevdata_long_water.R
@@ -0,0 +1,31 @@
+
+get_elevdata_long_water <- function(elevdata) {
+  dtc <- detect_water(elevdata %>% unlabel_elevdata)
+  elmat_df <- as.data.frame(dtc)
+  elmat_df <- cbind(elmat_df, elevdata[1:(nrow(elevdata) -1), ncol(elevdata)]) %>% as.data.frame
+  colnames(elmat_df) <- c(colnames(elevdata)[(ncol(elevdata) - 1):1], "deg_elmat_lat")
+  elmat_df <- rbind(elmat_df, elevdata[nrow(elevdata), c((ncol(elevdata) -1):1, ncol(elevdata))]) %>% as.data.frame
+
+
+  water_data <- get_elevdata_long(elmat_df)
+  real_data <- get_elevdata_long(elevdata)
+  real_data$value[which(water_data$value == 1)] <- min(real_data$value)
+
+  browser()
+  my_plot <- ggplot() +
+    geom_tile(
+      data = real_data,
+      aes_string("as.numeric(as.character(variable))","deg_elmat_lat",  fill = "value"),
+      alpha = 0.75) +
+    scale_y_continuous("Latitude", expand = c(0,0)) +
+    scale_fill_gradientn("Elevation", colours = terrain.colors(10)) +
+    # annotate(geom = 'tile', x =   water_data$variable, y =    water_data$deg_elmat_lat, 
+             # fill = (c("#FFFFFF00", "#63C600FF"))(water_data$value)) +
+    coord_fixed()
+  plot_gg(my_plot, shadow_intensity = 0.7, width = 5, height = 5, multicore = TRUE, scale = 350,
+          zoom = 0.5,
+          theta = 30,
+          phi = 60, windowsize = c(800, 800), 
+          raytrace = TRUE, saved_shadow_matrix = shadow_mat)
+  return(get_elevdata_long(elmat_df))
+}
diff --git a/R/read_arcgis_elev.R b/R/read_arcgis_elev.R
@@ -48,7 +48,7 @@ eudem_image_create <- function(bbox_arcgis, plot2d = TRUE, plot3d = FALSE) {
 
   legend_elevation <- seq(legend_max, legend_min, by = - (legend_max - legend_min) / (nrow(legend_matrix)-1))
 
-  img_elev <- get_arcgis_map_image(type = "elev", bbox = bbox_arcgis, width = 522, height = 800)
+    img_elev <- get_arcgis_map_image(type = "elev", bbox = bbox_arcgi, width = 522, height = 800)
   image_elev <- png::readPNG(img_elev)
 
   rgb_image <- image_elev[1, 1, c(1:3)]
@@ -77,13 +77,14 @@ eudem_image_create <- function(bbox_arcgis, plot2d = TRUE, plot3d = FALSE) {
 
   if (plot3d) {
 
+    # elevation_matrix <- elevation_matrix[, ncol(elevation_matrix):1]
+
     # Plotting the matrix as 3d
     elevation_matrix %>%
       sphere_shade(texture = "desert") %>%
       add_water(detect_water(elevation_matrix), color = "desert") %>%
       add_shadow(ray_shade(elevation_matrix, zscale = 3, maxsearch = 300), 0.5) %>%
-      plot_3d(elevation_matrix, zscale = 10, fov = 0, theta = 135, zoom = 0.75, phi = 45, windowsize = c(1000, 800))
-
+      plot_3d(elevation_matrix, zscale = 15, fov = 0, theta = 135, zoom = 0.75, phi = 45, windowsize = c(1000, 800))
   }
 
 }

diff --git a/R/video_animation.R b/R/video_animation.R
@@ -34,9 +34,21 @@ video_animation <- function(gpx_table = NULL, elevdata_long = NULL, number_of_sc
   video_indeces <- get_video_indeces(time_data = gpx_table$time_right, number_of_screens = number_of_screens)
 
   if (number_of_screens > 50) {
-    warning("theta and zoom parameters will be ignored")
-    theta_angles <- rev(30 - 50 * 1/(1 + exp(seq(-5, 6, length.out = length(video_indeces)))))
-    zoom_scale <- 0.5 + 0.5 * 1/(1 + exp(seq(-5, 5, length.out = length(video_indeces))))
+    if (length(theta) != 2) {
+
+      warning("theta parameters will be ignored")
+      theta_angles <- rev(30 - 50 * 1/(1 + exp(seq(-5, 6, length.out = length(video_indeces)))))
+    }else {
+      theta_angles <- rev(theta[1] - theta[2] * 1/(1 + exp(seq(-3, 3, length.out = length(video_indeces)))))
+    }
+
+    if (length(zoom) != 2) {
+      warning("zoom parameters will be ignored")
+      zoom_scale <- 0.5 + 0.5 * 1/(1 + exp(seq(-5, 5, length.out = length(video_indeces))))
+    } else {
+      zoom_scale <- zoom[1] + zoom[2] * 1/(1 + exp(seq(-5, 5, length.out = length(video_indeces))))
+
+    }
   } else{
     theta_angles <- rep(theta, length(video_indeces))
     zoom_scale <- rep(zoom, length(video_indeces))
@@ -94,7 +106,8 @@ video_animation <- function(gpx_table = NULL, elevdata_long = NULL, number_of_sc
     all_paths <- tempfile(fileext = ".txt")
 
     writeLines(con = all_paths,
-               paste0("file '",tempdir(), "\\video", 1:length(video_indeces), ".png'")
+               paste0("file '",tempdir(), "\\video", c(1:length(video_indeces),rep(length(video_indeces), 24)), ".png'")
+
     )
 
     outputfile <- tempfile(fileext = ".mp4")

diff --git a/R/video_animation_rayshade.R b/R/video_animation_rayshade.R
@@ -0,0 +1,96 @@
+video_animation_rayshade <- function(gpx_table, elevdata, number_of_screens = 500, theta = c(360, 340),
+                                     output_file_loc = tempfile(fileext = ".mp4"), places = NULL) {
+
+  video_indeces <- get_video_indeces(time_data = gpx_table$time_right, number_of_screens = number_of_screens)
+  elev_data_old <- elevdata
+
+  elevdata <- elevdata[c((nrow(elevdata) - 1):1, nrow(elevdata)), c((ncol(elevdata) - 1):1, ncol(elevdata))]
+  elevdata[nrow(elevdata), ] <- elevdata[nrow(elevdata), c((ncol(elevdata) - 1):1,ncol(elevdata)) ]
+  elevdata[, ncol(elevdata)] <- elevdata[c((nrow(elevdata)-1):1,nrow(elevdata)) , ncol(elevdata)]
+
+  colnames(elevdata) <- colnames(elevdata)[c((ncol(elevdata) - 1):1, ncol(elevdata))]
+
+  lon_elevdata <- as.numeric(colnames(elevdata)[(ncol(elevdata) - 1):1])
+  lat_elevdata <- as.numeric(elevdata$deg_elmat_lat[1:(nrow(elevdata) -1)])
+
+  gpx_tab_filtered <- gpx_table[video_indeces, ]
+
+  gpx_tab_filtered$lon_idx <- vapply(gpx_tab_filtered$lon, function(x) which.min(abs(x - lon_elevdata)), numeric(1))
+  gpx_tab_filtered$lat_idx <- vapply(gpx_tab_filtered$lat, function(x) which.min(abs(x - lat_elevdata)), numeric(1))
+
+  elevation_matrix <- elevdata %>% unlabel_elevdata() %>% t
+
+  gpx_tab_filtered$rel_speed_col <- scales::colour_ramp(viridisLite::viridis(10, option = "A", begin = 1, end = 0))(-gpx_tab_filtered$rel_speed / -max(gpx_tab_filtered$rel_speed))
+
+  gpx_tab_filtered$label <- rep(NA, nrow(gpx_tab_filtered))
+
+  place_indeces <- c()
+
+  if (!is.null(places) && number_of_screens > 50) {
+
+    for (row_index in 1:nrow(places)) {
+
+      place_indeces <- c(place_indeces, which.min(
+
+         sqrt((gpx_tab_filtered$lat - places$lat[row_index]) ^ 2 + (gpx_tab_filtered$lon - places$lon[row_index]) ^ 2))
+      )
+
+      gpx_tab_filtered$label[place_indeces[length(place_indeces)]] <- as.character(places$label[row_index])
+    }
+
+  }
+
+  image_size <- define_image_size(bbox_arcgis, major_dim = 600)
+  overlay_file <- tempfile(fileext = ".png")
+  get_arcgis_map_image(bbox_arcgis, map_type = "World_Topo_Map", file = overlay_file,
+                       width = image_size$width, height = image_size$height, 
+                       sr_bbox = 4326)
+  overlay_file_rot <- tempfile(fileext = ".png")
+  magick::image_write(magick::image_flop(magick::image_flip(image_read(path = overlay_file))), overlay_file_rot)
+  overlay_img <- png::readPNG(overlay_file_rot)
+
+  # Plotting the matrix as 3d
+  elevation_matrix %>%
+    sphere_shade(texture = "desert") %>%
+    add_water(detect_water(elevation_matrix), color = "desert") %>%
+    add_shadow(ray_shade(elevation_matrix, zscale = 3, maxsearch = 300), 0.5) %>%
+    add_overlay(overlay_img, alphalayer = 0.5) %>%
+    plot_3d(elevation_matrix, zscale = 15, fov = 1, theta = 135, zoom = 0.75, phi = 45, windowsize = c(1000, 800))
+
+  theta_angles <- rev(theta[1] - theta[2] * 1/(1 + exp(seq(-3, 3, length.out = length(video_indeces)))))
+
+  for (i in 1:nrow(gpx_tab_filtered)) {
+
+    render_label(elevation_matrix, x = gpx_tab_filtered[i, "lon_idx"], y = gpx_tab_filtered[i, "lat_idx"], z = 100, 
+                 zscale = 15, text = NULL, textsize = 15, linewidth = 6, freetype = FALSE, color = gpx_tab_filtered[i, "rel_speed_col"]) 
+    render_camera(theta = theta_angles[i])
+
+    if (!is.na(gpx_tab_filtered[i, "label"])) {
+      render_label(elevation_matrix, x = gpx_tab_filtered[i, "lon_idx"], y = gpx_tab_filtered[i, "lat_idx"], z = 600, 
+                   zscale = 15, text = gpx_tab_filtered[i, "label"],
+                   textsize = 15, linewidth = 5, freetype = FALSE, color = "black") 
+
+    }
+
+    render_snapshot(filename = file.path(tempdir(), paste0("video_rayshade_two", i, ".png")), clear = FALSE)
+  }
+  # ------ make it a movie -------
+  all_paths <- tempfile(fileext = ".txt")
+
+  writeLines(con = all_paths,
+             paste0("file '",tempdir(), "\\video_rayshade_two", c(1:length(video_indeces), rep(length(video_indeces), 48)), ".png'")
+
+  )
+
+  outputfile <- tempfile(fileext = ".mp4")
+
+  system(intern = TRUE,
+         paste0("ffmpeg ",
+                ifelse(overwrite, "-y", "-n"),
+                " -f concat -r 24 -safe 0 -i \"",
+                all_paths,
+                "\" -vf \"fps=24,format=yuv420p\" ", outputfile))
+  cat(outputfile)
+  file.copy(from = outputfile, to = output_file_loc, overwrite = overwrite)
+  return(output_file_loc)
+}
diff --git a/readme.txt b/readme.txt
@@ -0,0 +1,51 @@
+PROCESSED SRTM DATA VERSION 4.1
+
+The data  distributed here  are in  ARC GRID,  ARC ASCII  and Geotiff format, in
+decimal degrees and datum WGS84.  They are derived from the USGS/NASA SRTM data.
+CIAT  have  processed  this  data  to  provide  seamless  continuous  topography
+surfaces.  Areas with  regions of no  data in the  original SRTM data  have been
+filled using interpolation methods described by Reuter et al. (2007).
+
+Version 4.1 has the following enhancements over V4.0:
+- Improved ocean mask used, which includes some small islands  previously  being 
+  lost in the cut data.
+- Single no-data line of pixels along meridians fixed.
+- All GeoTiffs with 6000 x 6000 pixels.
+- For ASCII format files the projection definition is included in .prj files.
+- For GeoTiff format files the projection definition is in the .tfw  (ESRI TIFF 
+  World) and a .hdr file that reports PROJ.4 equivelent projection definitions.
+
+DISTRIBUTION
+
+Users are  prohibited from  any commercial,  non-free resale,  or redistribution
+without explicit written permission from CIAT. Users should acknowledge CIAT  as
+the source used  in the creation  of any reports,  publications, new data  sets,
+derived products, or services resulting from the use of this data set. CIAT also
+request  reprints of  any publications  and notification  of any  redistributing
+efforts.   For commercial  access to  the data,  send requests  to Andy   Jarvis
+(a.jarvis@cgiar.org).
+
+NO WARRANTY OR LIABILITY
+
+CIAT provides  these data  without any  warranty of  any kind whatsoever, either
+express or implied,  including warranties of  merchantability and fitness  for a
+particular purpose. CIAT shall not  be liable for incidental, consequential,  or
+special damages arising out of the use of any data.
+
+ACKNOWLEDGMENT AND CITATION
+
+We kindly ask  any users to  cite this data  in any published  material produced
+using this data,  and if possible  link web pages  to the CIAT-CSI  SRTM website
+(http://srtm.csi.cgiar.org).
+
+Citations should be made as follows:
+
+Jarvis A., H.I. Reuter, A.  Nelson, E. Guevara, 2008, Hole-filled  seamless SRTM
+data V4, International  Centre for Tropical  Agriculture (CIAT), available  from
+http://srtm.csi.cgiar.org.
+
+REFERENCES
+
+Reuter  H.I,  A.  Nelson,  A.  Jarvis,  2007,  An  evaluation  of  void  filling
+interpolation  methods  for  SRTM  data,  International  Journal  of  Geographic
+Information Science, 21:9, 983-1008.
diff --git a/srtm_39_03.hdr b/srtm_39_03.hdr
@@ -0,0 +1,32 @@
+Geotiff_Information:
+   Version: 1
+   Key_Revision: 1.0
+   Tagged_Information:
+      ModelTiepointTag (2,3):
+         0                0                0                
+         10               50               0                
+      ModelPixelScaleTag (1,3):
+         0.000833333333   0.000833333333   0                
+      End_Of_Tags.
+   Keyed_Information:
+      GTModelTypeGeoKey (Short,1): ModelTypeGeographic
+      GTRasterTypeGeoKey (Short,1): RasterPixelIsArea
+      GeographicTypeGeoKey (Short,1): GCS_WGS_84
+      GeogCitationGeoKey (Ascii,7): "WGS 84"
+      GeogAngularUnitsGeoKey (Short,1): Angular_Degree
+      End_Of_Keys.
+   End_Of_Geotiff.
+
+GCS: 4326/WGS 84
+Datum: 6326/World Geodetic System 1984
+Ellipsoid: 7030/WGS 84 (6378137.00,6356752.31)
+Prime Meridian: 8901/Greenwich (0.000000/  0d 0' 0.00"E)
+
+PROJ.4 Definition: +proj=latlong +ellps=WGS84 +to_meter=1.0000000000
+
+Corner Coordinates:
+Upper Left    (10.0000000,50.0000000)
+Lower Left    (10.0000000,45.0000000)
+Upper Right   (15.0000000,50.0000000)
+Lower Right   (15.0000000,45.0000000)
+Center        (12.5000000,47.5000000)
diff --git a/srtm_39_03.tfw b/srtm_39_03.tfw
@@ -0,0 +1,6 @@
+            0.0008333333
+            0.0000000000
+            0.0000000000
+           -0.0008333333
+           10.0004166667
+           49.9995833333
diff --git a/srtm_39_03.tif b/srtm_39_03.tif
diff --git a/vignettes/create_video.Rmd b/vignettes/create_video.Rmd
@@ -113,7 +113,7 @@ into a `mp4` by `ffmpeg`. Therefore `ffmpeg` needs to be installed and inside th
 
 
 ```{r, eval = FALSE, echo = TRUE}
-output_gif <- video_animation(gpx_table = gpx_table, elevdata_long = elmat_long, make_gif = TRUE, number_of_screens = 5,
+output_gif <- video_animation(gpx_table = gpx_table, elevdata_long = elmat_long, make_gif = TRUE, number_of_screens = 4,
                 output_file_loc = tempfile(fileext = ".gif"))
 
 Sys.setenv("PATH" = paste0(Sys.getenv("PATH"), ";", "<Where is ffmpeg>"))

diff --git a/vignettes/readme.txt b/vignettes/readme.txt
@@ -0,0 +1,51 @@
+PROCESSED SRTM DATA VERSION 4.1
+
+The data  distributed here  are in  ARC GRID,  ARC ASCII  and Geotiff format, in
+decimal degrees and datum WGS84.  They are derived from the USGS/NASA SRTM data.
+CIAT  have  processed  this  data  to  provide  seamless  continuous  topography
+surfaces.  Areas with  regions of no  data in the  original SRTM data  have been
+filled using interpolation methods described by Reuter et al. (2007).
+
+Version 4.1 has the following enhancements over V4.0:
+- Improved ocean mask used, which includes some small islands  previously  being 
+  lost in the cut data.
+- Single no-data line of pixels along meridians fixed.
+- All GeoTiffs with 6000 x 6000 pixels.
+- For ASCII format files the projection definition is included in .prj files.
+- For GeoTiff format files the projection definition is in the .tfw  (ESRI TIFF 
+  World) and a .hdr file that reports PROJ.4 equivelent projection definitions.
+
+DISTRIBUTION
+
+Users are  prohibited from  any commercial,  non-free resale,  or redistribution
+without explicit written permission from CIAT. Users should acknowledge CIAT  as
+the source used  in the creation  of any reports,  publications, new data  sets,
+derived products, or services resulting from the use of this data set. CIAT also
+request  reprints of  any publications  and notification  of any  redistributing
+efforts.   For commercial  access to  the data,  send requests  to Andy   Jarvis
+(a.jarvis@cgiar.org).
+
+NO WARRANTY OR LIABILITY
+
+CIAT provides  these data  without any  warranty of  any kind whatsoever, either
+express or implied,  including warranties of  merchantability and fitness  for a
+particular purpose. CIAT shall not  be liable for incidental, consequential,  or
+special damages arising out of the use of any data.
+
+ACKNOWLEDGMENT AND CITATION
+
+We kindly ask  any users to  cite this data  in any published  material produced
+using this data,  and if possible  link web pages  to the CIAT-CSI  SRTM website
+(http://srtm.csi.cgiar.org).
+
+Citations should be made as follows:
+
+Jarvis A., H.I. Reuter, A.  Nelson, E. Guevara, 2008, Hole-filled  seamless SRTM
+data V4, International  Centre for Tropical  Agriculture (CIAT), available  from
+http://srtm.csi.cgiar.org.
+
+REFERENCES
+
+Reuter  H.I,  A.  Nelson,  A.  Jarvis,  2007,  An  evaluation  of  void  filling
+interpolation  methods  for  SRTM  data,  International  Journal  of  Geographic
+Information Science, 21:9, 983-1008.