Making 3D Maps of Sentinel-2 Imagery with Rayshader

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 . You can register to it. It is a free and a very user-friendly apps.

User interface of the EO Browser

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: .

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.


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

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.

Plot of sentinel-2 image

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")

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() %>% 

Afterward we will find our digital elevation model from SRTM using elevatr developed by Jeffrey Hollister from US EPA ( ). 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 ( ).

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

 leaflet(mask) %>% 
   addTiles() %>% 
   addPolygons() %>% 
Downloaded elevation raster and mask area

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)
Cropped elevation raster

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)
 elevation_matrix = rayshader::raster_to_matrix(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)
Render_snapshot result of Karangetang Volcano, North Sulawesi, Indonesia
render_highquality result gives more dramatic effect

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!

Memantau Kebocoran Minyak di Laut Jawa

Minyak bocor di Laut Jawa! Kasus ini diungkap oleh Direktur Hulu PT Pertamina Persero Dharmawan Samsu bahwa insiden ini telah terjadi sejak 2 minggu lalu, 12 Juli 2019 dini hari. Ketika itu muncul gelembung gas, yang menyebabkan operasi pengeboran berhenti dan tanggal 14 Juli Pertamina memutuskan evakuasi pegawai.

Pertamina pertama membuka rilis insiden ini pada tanggal 18 Juli dan mengatakan telah melakukan upaya maksimum untuk menghentikan kebocoran gas dan minyak.

Menurut rilis CNCB Indonesia, per tanggal 22 Juli, bocoran minyak telah meyebar hingga ke pantai di pesisir Karawang. Pertamina juga telah mengerahkan 27 kapal yang dilengkapi dengan oil boom dan dispersant ke sekitar anjungan dan titik-titik yang terindikasi terdapat tumpahan minyak di laut.

Yang menarik peristiwa kebocoran minyak ini bisa kita amati jelas dari citra satelit dengan teknik penginderaan jauh.

Dengan menggunakan citra satelit Sentinel 2 dari European Satelite Agency (ESA), kita bisa jelas perkirakan kebocorannya. Saya sebenarnya tidak tahu pasti di mana titik kebocoran dan titik anjungan. Jadi titik ini hanya dugaan, berdasarkan fakta-fakta visual yang jelas.

Dengan piranti lunak Sentinel Hub, Saya mengarahkan pencarian saya ke Laut Jawa di pesisir Karawang. Kemudian mengatur tanggal antara awal Juli hingga sekarang. Dari hasil ini, saya menduga-duga di mana titik anjungan berada. Hasilnya adalah titik dugaan kebocoran ini berada sekitar 10 kilometer arah utara dari Muara Cilamaya. Kemudian saya membuat gambar selang waktunya.

Perhatikan bagian tengah kanan gambar. Sebaran minyak mulai terlihat pada tanggal 14 Juli 2019, kemudian kita lihat banyak titik-titik yang kemungkinan besar adalah kapal-kapal pembersih minyak. Bagian pojok kiri bawah peta adalah pantai di Karawang.
Perbesaran gambar sebelumnya. Kita bisa amati lebih jelas kebocoran minyak dan upaya pembersihan minyak pada tanggal 24 Juli 2019.

Dengan menggunakan citra Sentinel 1 yang berbasis radar, kita bisa perkirakan luas penyebaran tumpahan minyak. Ini memungkinkan karena minyak memiliki respon berbeda terhadap gelombang radar dibandingkan air.

Citra Sentinel 1 Ground Range Detected mode IW tanggal 18 Juli 2019. Mohon maaf tidak ada skala.
Citra Sentinel 1 Ground Range Detected mode IW tanggal 30 Juli 2019. Skala tidak sama dengan gambar tanggal 18 Juli.
Timelapse sebelum kebocoran dan sesudah kebocoran. Perhatikan warna gelap yang merupakan refleksi dari minyak.

Di gambar kita bisa lihat minyak kemungkinan besar yang berwarna gelap, sementara air laut bertekstur kasar. Titik awal warna gelap berasal adalah sama dengan titik di gambar selang waktu sebelumnya. Titik-titik terang kemungkinan besar adalah kapal-kapal yang berseliweran membersihkan tumpahan minyak.

Dengan memanfaatkan citra satelit, kita bisa menghitung seberapa luas pelamparan minyak. Kita juga bisa menggunakan data ini untuk merencanakan strategi penanganan dampak sehingga upaya kita bisa maksimal.

Semoga tulisan ini mampu memberikan informasi yang mencerahkan. Tulisan ini tidak saya maksudkan untuk menyudutkan siapa-siapa. Hanya semata-mata keinginan untuk berbagi cara melihat bumi dari sudut pandang yang lain. Untuk itu, informasi ini mohon untuk tidak dikutip dan mari kita menunggu rilis resmi dari LAPAN sebagai pihak yang berkepentingan di bidang ini.

Saya mendoakan yang terbaik untuk segenap jajaran Pertamina dan kru yang bertugas. Semoga kejadian ini yang terakhir dan tidak terulang lagi.

Semoga penanganan kebocoran dan pembersihan laut terlaksana sebaik-baiknya.

Gambar terbaru dari Sentinel 2 tanggal 29 Juli 2019!