3D Plot Hemis National Park

Getting all the necessary Libraries

library(here)

## here() starts at D:/R projects/Rayshader/3D-Maps-of-Satellite-Images

library(tidyverse)

## -- Attaching packages ----------------------- tidyverse 1.3.0 --

## v ggplot2 3.3.0     v purrr   0.3.4
## v tibble  3.0.1     v dplyr   0.8.5
## v tidyr   1.0.3     v stringr 1.4.0
## v readr   1.3.1     v forcats 0.5.0

## Warning: package 'ggplot2' was built under R version 3.6.3

## Warning: package 'tibble' was built under R version 3.6.3

## Warning: package 'purrr' was built under R version 3.6.3

## Warning: package 'dplyr' was built under R version 3.6.3

## Warning: package 'forcats' was built under R version 3.6.3

## -- Conflicts -------------------------- tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

library(raster)

## Warning: package 'raster' was built under R version 3.6.3

## Loading required package: sp

## Warning: package 'sp' was built under R version 3.6.3

## 
## Attaching package: 'raster'

## The following object is masked from 'package:dplyr':
## 
##     select

## The following object is masked from 'package:tidyr':
## 
##     extract

library(rayshader)

## Warning: package 'rayshader' was built under R version 3.6.3

library(scales)

## Warning: package 'scales' was built under R version 3.6.3

## 
## Attaching package: 'scales'

## The following object is masked from 'package:purrr':
## 
##     discard

## The following object is masked from 'package:readr':
## 
##     col_factor

library(vembedr)

## Warning: package 'vembedr' was built under R version 3.6.3

Elevation data

elevation1 = raster::raster(here("elevation data","N33E077.hgt"))
elevation2 = raster::raster(here("elevation data","N34E077.hgt"))

hemis_elevation = raster::merge(elevation1,elevation2)

height_shade(raster_to_matrix(hemis_elevation)) %>%
  plot_map()

Loading georeferenced satellite imagery

I have no clue what this means.

hemis_r = raster::raster(here("landsat8","LC08_L1TP_147037_20190926_20191017_01_T1_B4.TIF"))
hemis_g = raster::raster(here("landsat8","LC08_L1TP_147037_20190926_20191017_01_T1_B3.TIF"))
hemis_b = raster::raster(here("landsat8","LC08_L1TP_147037_20190926_20191017_01_T1_B2.TIF"))

hemis_rbg = raster::stack(hemis_r, hemis_g, hemis_b)
raster::plotRGB(hemis_rbg, scale=255^2)

Darkness correction

The image in the preceeding section look as dark to as the one in the blog. Will try out the correction and see how it goes.

hemis_rbg_corrected = sqrt(raster::stack(hemis_r, hemis_g, hemis_b))
raster::plotRGB(hemis_rbg_corrected)

This seems too light to me. Let’s see how the contrast adjustment goes(later).

Combining landsat and elevation data

hemis_elevation_utm = raster::projectRaster(hemis_elevation, crs = crs(hemis_r), method = "bilinear")
crs(hemis_elevation_utm)

## CRS arguments:
##  +proj=utm +zone=43 +datum=WGS84 +units=m +no_defs +ellps=WGS84
## +towgs84=0,0,0

Croping the area of the park

bottom_left = c(y=77.135192, x=33.437587)
top_right   = c(y=77.455761, x=33.690335)

extent_latlong = sp::SpatialPoints(rbind(bottom_left, top_right), proj4string=sp::CRS("+proj=longlat +ellps=WGS84 +datum=WGS84"))
extent_utm = sp::spTransform(extent_latlong, raster::crs(hemis_elevation_utm))

e = raster::extent(extent_utm)
e

## class      : Extent 
## xmin       : 698488.5 
## xmax       : 727630.8 
## ymin       : 3701838 
## ymax       : 3730529

There are things that I dont get and there are things that I dont get at all

hemis_rgb_cropped = raster::crop(hemis_rbg_corrected, e)
elevation_cropped = raster::crop(hemis_elevation_utm, e)

names(hemis_rgb_cropped) = c("r","g","b")

hemis_r_cropped = rayshader::raster_to_matrix(hemis_rgb_cropped$r)
hemis_g_cropped = rayshader::raster_to_matrix(hemis_rgb_cropped$g)
hemis_b_cropped = rayshader::raster_to_matrix(hemis_rgb_cropped$b)

hemisel_matrix = rayshader::raster_to_matrix(elevation_cropped)

hemis_rgb_array = array(0,dim=c(nrow(hemis_r_cropped),ncol(hemis_r_cropped),3))

hemis_rgb_array[,,1] = hemis_r_cropped/255 #Red layer
hemis_rgb_array[,,2] = hemis_g_cropped/255 #Blue layer
hemis_rgb_array[,,3] = hemis_b_cropped/255 #Green layer

hemis_rgb_array = aperm(hemis_rgb_array, c(2,1,3))

plot_map(hemis_rgb_array)

Adjusting contrast

hemis_rgb_contrast = scales::rescale(hemis_rgb_array,to=c(0,1))

plot_map(hemis_rgb_contrast)

3D PLot

Commenting the below commands as these are working fine while executed chunk by chunk, but creates a blank plot while kniting. The generated plot is displayed below.

# plot_3d(hemis_rgb_contrast, hemisel_matrix, windowsize = c(1100,900), zscale = 15, shadowdepth = -50,
#         zoom=0.5, phi=45,theta=-45,fov=70, background = "#F2E1D0", shadowcolor = "#523E2B")
# render_snapshot( title_text = "Home of the Sonw Leopards | Imagery: Landsat 8 | DEM: 30m SRTM",
#                 title_bar_color = "#b2905d", title_color = "white", title_bar_alpha = 0.5)

Searching Snow Leopards

Making the video

# angles= seq(0,360,length.out = 1441)[-1]
# for(i in 1:1440) {
#  render_camera(theta=-45+angles[i])
#  render_snapshot(filename = sprintf("hemispark%i.png", i),
#                  title_text = "HOme of the Snow Leopards| Imagery: Landsat 8 | DEM: 30m SRTM",
#                  title_bar_color = "#1f5214", title_color = "white", title_bar_alpha = 1)
# }
# rgl::rgl.close()
# 
# 
# list <- list.files(pattern = 'hemispark')
# images <- image_read(list)
# animationconfirmed <- image_animate(images,fps = 100,optimize=T)
# image_write(animationconfirmed, 'Hemis National Park.gif')

The loop was executed to generate 1440 png images, however I used windws application to make a video out of it. The video is below.