Tag: sentinel 2

I am inspired by Tyler Morgan-Wall tutorial on Rayshader with satellite imagery overlay published few days ago. I have been trying this for some times, but somehow I always fail. So when Tyler shares his code, I immediately follow and try it on my own.

The code is so elegant, one can simply follow and have their own beautiful rendered 3d images with satellite imagery overlies. I really suggest that you visit Tyler’s blog. There are plenty of amazing tutorials there.

In this post, I would like to add something that I believe is useful when anyone wants to create a 3d images with satellite imagery overlies with rayshader.

I add two things that a little bit different than Tyler’s step to create the images, which are: 1) download sentinel-2 images from Sentinel-Hub apps , and 2) download SRTM images from R . These two steps are crucial and can save some time and avoiding us to download a huge data of satellite images.

First we will download the images of Sentinel-2. The regions that I want to share is Karangetang Volcano in North Sulawesi, Indonesia. This volcano is interesting because it has been erupting for the past few years and we can see the trace of the eruption finely from Sentinel-2 images.

Sentinel-2 is a satellite program by European Space Agency. It has been acquiring images from the orbit since 2015 with spatial resolution of 10m. We can access Sentinel-2 images through this apps https://apps.sentinel-hub.com/eo-browser/ . You can register to it. It is a free and a very user-friendly apps.

Above is the interface of EO Browser after you sign-in. On the left side is the menu for selecting the satellite, cloud coverage, and time range. You can select Sentinel-2 and pick L2A as this images is already processed. Set the cloud coverage to the lowest so you can get a free cloud image.

On the right side, there is a marker symbol. This is a pin to select point of interest. It also allows us to make an area of interest. So I make a rectangle. Afterward you can click “Search” and EO Browser will find you list of available data set.

After selecting an images, EO Browser will show you available image processing method, which can be used for many purpose. You can also develop your own custom script if you will. You can check for available script here: https://custom-scripts.sentinel-hub.com/ .

For our purpose, lets just do True Color image. You can download the image on the right side, just two symbol below the marker symbol, called “Download image”. Unlike Earth Explorer that force us to download all the data in bulk, EO Browser offers us the chance to download only the data we need, and we can crop it to the area we made earlier. You can select the Image format to Tiff 8 bit, image resolution high, and coordinate system to WGS84. Select layer True Color.

Afterward you can immediately download the file, and the size will not be huge (still in the order of tens of MB). Now we can call it in R.

 library(raster)
 library(sp)
 library(rayshader)
 library(elevatr)
 library(leaflet)  

 satellite_images <- raster::brick("Sentinel-2 L2A from 2019-09-22.tiff")
 plotRGB(satellite_images)

So our data consist of one tiff file. We call it using brick. To be honest I am not really sure the reason why raster::brick work, and raster::raster doesn’t. So just use raster::brick for this purpose.

Second, we will create our Area of Interest. We will create a simple AoI, a rectangle with 4 edges. For that I make a function, so you just need to add xleft, ytop, xright, and ybot, or two coordinates of top left and bottom right edges.

AOI <- function(p, q, r, s) {
   xleft <- p
   ytop <- q
   xright <- r
   ybot <- s
   x_coord <- c(xleft,xleft,xright,xright)
   y_coord <- c(ytop,ybot, ybot, ytop)
   xym <- cbind(x_coord, y_coord)
   p <- Polygon(xym)
   ps = Polygons(list(p),1)
   sps = SpatialPolygons(list(ps))
   proj4string(sps) = CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
   return(sps)
 }

You can call your AoI and then afterward check it on leaflet to see where it lies. If it is ok then you can continue to the next step. Also make sure that your Sentinel-2 images are completely within your study area.

mask <- AOI(xleft, ytop, xright, ybot)
leaflet(mask) %>% 
   addTiles() %>% 
   addPolygons() 

Afterward we will find our digital elevation model from SRTM using elevatr developed by Jeffrey Hollister from US EPA (https://github.com/jhollist/elevatr ). This packages allows us to download elev_point and elev_raster provided by USGS and USNOAA. The DEM needs to be projected so it the projection system will be the same as our satellite image. We can plot it on the map using leaflet package (https://rstudio.github.io/leaflet/ ).

 Elevation_File <- get_elev_raster(AOI1, z=11)
 projectRaster(Elevation_File, 
               crs= "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")

 leaflet(mask) %>% 
   addTiles() %>% 
   addPolygons() %>% 
   addRasterImage(Elevation_File)

I cropped my file using crop function. I clipped it using AoI which I named “mask”. Afterward I plotted it.

elevation_mask <- crop(Elevation_File, mask)
 satellite_imagery_mask <- crop(satellite_imagery, mask)
 plot1 <- leaflet() %>% 
   addTiles() %>% 
   addRasterImage(elevation_mask, opacity =0.5)
 plot1

From here and after I just follow Tyler’s tutorial. So I make my data to be matrix. Satellite data has to be differentiate as r,g,b because that is how color has to be treated (color is combination of red, green, and blue bands).

Because it is too complicated for me to explain, so I will just quote Tyler Morgan-Wall himself for this section:

“Now we’ll crop our datasets to the same region, and create an 3-layer RGB array of the image intensities. This is what rayshader needs as input to drape over the elevation values. We also need to transpose the array, since rasters and arrays are oriented differently in R, because of course they are🙄. We do that with the aperm() function, which performs a multi-dimensional transpose. We’ll also convert our elevation data to a base R matrix, which is what rayshader expects for elevation data.”

Tyler Morgan-Wall

names(satellite_imagery_mask) = c("r","g","b")
 satellite_imagery_mask_r = rayshader::raster_to_matrix(satellite_imagery_mask$r)
 satellite_imagery_mask_g = rayshader::raster_to_matrix(satellite_imagery_mask$g)
 satellite_imagery_mask_b = rayshader::raster_to_matrix(satellite_imagery_mask$b)
 satellite_imagery_mask
 elevation_matrix = rayshader::raster_to_matrix(elevation_mask)
 elevation_mask
 sentinel_mask_array = array(0,dim=c(nrow(satellite_imagery_mask_r),ncol(satellite_imagery_mask_r),3))
 sentinel_mask_array[,,1] = satellite_imagery_mask_r/255 #Red layer
 sentinel_mask_array[,,2] = satellite_imagery_mask_g/255 #Blue layer
 sentinel_mask_array[,,3] = satellite_imagery_mask_b/255 #Green layer
 sentinel_mask_array = aperm(sentinel_mask_array, c(2,1,3))
 plot_3d(sentinel_mask_array, elevation_matrix, windowsize = c(1100,900), zscale = 30, shadowdepth = -50,
         zoom=0.5, phi=45,theta=-45,fov=70, background = "#F2E1D0", shadowcolor = "#523E2B")
 render_snapshot(title_text = "Karangetang Volcano, Indonesia | Imagery: Sentinel-2 | DEM: 30m SRTM",
                   title_bar_color = "#1f5214", title_color = "white", title_bar_alpha = 1)

Voilla, a beautiful 3D image are made, easily. Thanks to amazing Rayshader package! In his tutorial, Tyler Morgan-Wall can make a gif, however it is not possible for me due to my laptop limitation. Rendering process is painful and unbearable for her. But you can try.

You can check the code of this work here.

Have fun with Rayshader!