tuto_satellite_relief.R

library(tidyverse)
library(sf)
library(raster)
library(rayshader)

#on crée une fonction avec X=longitude, Y=latitude, (les deux en coordonnées géographiques), le titre de la vignette,  une distance autour du point central, par défaut à 2km, et une échelle de mise en relief (zs).

SortLaCarte<-function(X,Y,titre,distanceautourpoint=2000,repartirde,zs=6){
  Lieu_4326<-tibble(nom=titre,y=Y,x= X)%>%
    st_as_sf(coords=c("x","y"),crs=4326)%>%
    st_transform(crs=2154)%>%
    st_buffer(dist = distanceautourpoint)%>%
    st_transform(crs=4326)
  Extent_spdf<-as_Spatial(Lieu_4326)
  
  # bounding box
  bbox_xy_wgs1984 <- st_bbox(Lieu_4326)
  
  # coordonnées des coins
  xmin_bbox <- bbox_xy_wgs1984$xmin %>% as.vector()
  xmax_bbox <- bbox_xy_wgs1984$xmax %>% as.vector()
  ymin_bbox <- bbox_xy_wgs1984$ymin %>% as.vector()
  ymax_bbox <- bbox_xy_wgs1984$ymax %>% as.vector()
  
  #on utilise l'api de l'IGN sur l'élévation
  url_img_elev <- paste0("https://wxs.ign.fr/inspire/inspire/r/wms?LAYERS=EL.GridCoverage&EXCEPTIONS=text/xml&FORMAT=image/jpeg&SERVICE=WMS&VERSION=1.3.0&REQUEST=GetMap&STYLES=&CRS=EPSG:4326&BBOX=",ymin_bbox,",",xmin_bbox,",",ymax_bbox,",",xmax_bbox,"&WIDTH=1024&HEIGHT=1024")
  
  #on télécharge l'image
  download.file(url = url_img_elev,
                destfile = paste0("data/jpg/elevation_pt_",Y,"_",X, ".jpg"),
                mode = 'wb')
  #on la lit
  raster_tel<-raster::raster(paste0("data/jpg/elevation_pt_",Y,"_",X, ".jpg"),layer=0)
  #on lui définit une projection...
  crs(raster_tel)<-"+proj=longlat +datum=WGS84 +no_defs"
  #... des points extremes
  extent(raster_tel)<-c(xmin_bbox,xmax_bbox,ymin_bbox,ymax_bbox)
  proj4_2154<-"+proj=lcc +lat_0=46.5 +lon_0=3 +lat_1=49 +lat_2=44 +x_0=700000 +y_0=6600000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
  #On reprojette en lambert 93
  elevation_2154 = raster::projectRaster(raster_tel, crs = proj4_2154, method = "bilinear")
  # on crée une matrice d'élévation (voir le site de rayshader pour plus de détails)
  elmat = matrix(raster::extract(elevation_2154,raster::extent(elevation_2154),buffer=1000),
                 nrow=ncol(elevation_2154),ncol=nrow(elevation_2154))
  
  #on fait pareil pour l'api des orthophotos
  # téléchargement des images
  url_img_ortho <- paste0("https://wxs.ign.fr/ortho/geoportail/r/wms?LAYERS=ORTHOIMAGERY.ORTHOPHOTOS.BDORTHO&EXCEPTIONS=text/xml&FORMAT=image/jpeg&SERVICE=WMS&VERSION=1.3.0&REQUEST=GetMap&STYLES=&CRS=EPSG:4326&BBOX=",ymin_bbox,",",xmin_bbox,",",ymax_bbox,",",xmax_bbox,"&WIDTH=1256&HEIGHT=1256")
  
  # si l'image existe on la lit, sinon, on la télécharge
  if(!file.exists(paste0("./data/jpg/img_ortho_actu_poi_",Y,"_",X, ".jpg"))){
    download.file(url = url_img_ortho,
                  destfile = paste0("./data/jpg/img_ortho_actu_poi_",Y,"_",X, ".jpg"),
                  mode = 'wb')
  }
  
  ## on met la photo au dessus de l'élévation
  # en lisant les bandes rouges, vertes et bleues successivement
  img_ortho_1 = raster::raster(paste0("./data/jpg/img_ortho_actu_poi_",Y,"_",X, ".jpg"), 
                               band=1)
  img_ortho_2 = raster::raster(paste0("./data/jpg/img_ortho_actu_poi_",Y,"_",X, ".jpg"), 
                               band=2)
  img_ortho_3 = raster::raster(paste0("./data/jpg/img_ortho_actu_poi_",Y,"_",X, ".jpg"), 
                               band=3)
  img_rgb = raster::stack(img_ortho_1, img_ortho_2, img_ortho_3)
  crs(img_rgb)<-"+proj=longlat +datum=WGS84 +no_defs"
  img_rgb_ok<-setExtent(img_rgb,extent(Lieu_4326))
  img_sat = raster::brick(img_rgb_ok)
  
  img_rgb_array = array(0,dim=c(nrow(img_sat),ncol(img_sat),3))
  
  img_rgb_array[,,1] = img_sat$layer.1@data@values/255 #Red layer
  img_rgb_array[,,2] = img_sat$layer.2@data@values/255 #Blue layer
  img_rgb_array[,,3] = img_sat$layer.3@data@values/255 #Green layer
  
  img_rgb_array[,,1][img_rgb_array[,,1]>1]<-1
  img_rgb_array[,,2][img_rgb_array[,,2]>1]<-1
  img_rgb_array[,,3][img_rgb_array[,,3]>1]<-1
  
  img_rgb_array = aperm(img_rgb_array, c(2,1,3))
  
  #On crée un tour autour du centre
  Plans<-tibble(org=1:90,
                zoom=c(seq(from=60, to=30,length.out=45),seq(from=30, to=60,length.out=45)),
                phi=c(seq(from=90,to=30,length.out=30),seq(from=30,to=15,length.out=30),seq(from=15,to=90,length.out=30)),
                theta=c(seq(from=0,to=180,length.out=45),seq(from=180,to=360,length.out=45)))
  rgl::clear3d()
  
  #si on veut tout sortir, on sort toutes les 90 images
  # for (i in 1:nrow(Plans)){
  #si on veut juste voir quelques éléments
  for (i in c(1,20,45,60,80)){
    
    plot_3d(img_rgb_array,elmat, zscale=zs,fov=0,theta=Plans$theta[i],zoom=Plans$zoom[i]/100,phi=Plans$phi[i], windowsize = c(900,900),
            water=FALSE, waterdepth = 0, wateralpha = 0.4,watercolor = "#5C948B",
            solidcolor = "#141E28",waterlinecolor = "white",waterlinealpha = 1,background = "#F2E1D0", shadowcolor = "#523E2B")
    
    
    render_snapshot(title_text = titre,
                    title_bar_color = "#F2E1D0", title_color = "white", title_bar_alpha = 1,
                    filename = paste0("data/sorties/test",titre,"_",str_pad(i,width=3,pad="0")),width = 900,height=900)  
  }
  
}
#on fait ensuite tourner cette fonction sur un point
SortLaCarte(X = -3.9456440913109434,Y = 48.34998510742631,titre = "Mont Saint-Michel, Finistère", distanceautourpoint = 3000,zs=4)