Khumbu icefall in 4D
This video shows all available Venµs images since launch over the Khumbu icefall near Mt Everest (109 images). The images are RGB composites of Level 1C bands 7,4,3 (5 m resolution) that were draped onto the HMA 8m DEM.
How to do that? First download the Venµs images using theia_download.py, extract the region of interest and make color composites using gdal_translate as explained here. With a recent GDAL installation you can even read ZIP archives on-the-fly without decompressing them. Then, prepare an elevation raster from the HMA 8m DEM so that it has the same resolution and extent as the extracted Venµs images. Finally make the 3D view in R using rayshader.
library(rayshader)
library(abind)
library(raster)
p="./"
fdem=paste0(p,"HMA_DEM8m_MOS_20170716_tile-677_UTM45_filled_crop_Venus.tif")
pimg=paste0(p,"/VIS1C/")
dem=raster(fdem)
elmat=matrix(extract(dem, extent(dem), buffer = 1000),nrow = ncol(dem), ncol = nrow(dem))
files=list.files(path=pimg, pattern="jpg$", full.names=FALSE, recursive=FALSE)
# rotating view
th=200
k=0
for (x in files) {
fimg=paste0(pimg,x)
fout=paste0(p,"/3Drot/",x)
print(x)
txt=substr(x, 36, 43)
k=k+1
if (k<length(files)/2){
th=th+0.5
}else{
th=th-0.5
}
imgr = raster(fimg,band=1)
imgg = raster(fimg,band=2)
imgb = raster(fimg,band=3)
r = t(matrix(extract(imgr, extent(imgr), buffer = 1000),nrow = ncol(imgr), ncol = nrow(imgr)))
g = t(matrix(extract(imgg, extent(imgr), buffer = 1000),nrow = ncol(imgr), ncol = nrow(imgr)))
b = t(matrix(extract(imgb, extent(imgr), buffer = 1000),nrow = ncol(imgr), ncol = nrow(imgr)))
col=abind(r,g,b,along=3)/255
plot_3d(col, elmat, zscale = 5, fov = 0, theta = th, zoom = 0.7, phi = 45, windowsize = c(1000, 800))
render_snapshot(fout,clear = TRUE,title_text = txt, title_color = "black",
title_font = "Arial", gravity = "NorthEast", title_offset = c(0,0))
}
Merci à Gérard Dedieu pour l’idée !
Here is the same animation without rotation: