Ever since I released rayshader to the public, there's been one question that comes up time and time again: "How do I use rayshader to overlay satellite imagery onto a 3D surface?" And I can see why: 3D maps with real satellite images are super cool, and there are few detailed tutorials out there on how it's done in R. Additionally, the process behind making this type of map can be intimidating for non-GIS experts, as it exposes you to all complexities of the GIS field: combining different datasets from separate sources—often each with their own distinct coordinate systems—into a single map. Not only that, but you have to also know where and how to get this data—not obvious to someone who doesn't do this regularly. So I'm going to walk you through how to obtain the data required to make these types of maps, as well as the R code used to generate them. In my opinion, the hardest part is all the manual work in downloading the data! Once that's done, the code required is short and straightforward, and rayshader makes visualizing elevation data in 3D a breeze. Let’s get started!
First, we'll download elevation data. A good place to start for reasonable resolution elevation data is the Shuttle Radar Topography Mission (SRTM) dataset. It’s a global 30 meter resolution dataset (meaning, 30m between each point), and although there are more precise and higher resolution datasets out there, few have the global coverage of the SRTM–and it’s dead easy to get the data. Specifically, there’s simple tile selectors out there to download the data by just clicking on a map, which is a far easier process than many other online data repositories.
We are going to be visualizing Zion National Park, UT, which has beautiful, dramatic topography. I took a trip with my Dad1 out to Zion after I finished my PhD, and the gorgeous landscape made up of deep valleys and sheer cliffs makes for an excellent 3D visualization. Let’s start by grabbing elevation data from Derek Watkin’s SRTM tile downloader. We’ll zoom into tiles N37W113.hgt
and N37W114.hgt
, since the park overlaps both of them. To download: click the tiles, click download, and unzip the files (they should have the file type “hgt”).

Now, we need to get satellite imagery data. This is a bit more involved. There are plenty of sources for this type of data, both paid and freely available. We are going to be taking data from the Landsat 8 project, which does require registering a free USGS account, but then we get access to up-to-date high-quality satellite imagery, for free. So let me walk you through the the step-by-step process to obtain this data.
We are going to go to the USGS Earth Explorer. Once you register an account, you can download datasets. It takes a few minutes to receive your confirmation email, but it comes. Once you get your confirmation email, you can dive right in.

Here’s the interface, with Utah in the viewer. We’re going to zoom in closer to the bottom left corner, which is the location of Zion National Park. Let’s do that:

Now that we’re zoomed into the region we want, we’re going to hit the “Use Map” button to mark the current map view as the bounding box for the area we want to search for. We’ll then click “Data Sets” in the bottom left corner. This will query the database and return datasets that include the region we have selected. We specifically want Landsat 8 data, so we’ll select the following:

Now, we’ll click “Results”, giving us a list of data products we can select from. If you click the little image icon (second to the left), it will show you a preview of the data in the window to the right. Bluish areas indicate cloud cover, so we’ll flip through the search results until we get a date with a clear shot of the terrain. The first one I find is “LC08_L1TP_038034_20191101_20191114_01_T1”, taken on November 1st, 2019 (or, about 500 years ago in corona-time). Let’s download the full resolution Level 1 GeoTIFF data product (it's big, almost a gigabyte) and unzip it.

And now–we write our script to turn this into a 3D model! Getting the data was the hardest part. That is, unless you can’t install the following packages–if that's the case, godspeed, it’s a rite of passage to have to fight with your local GDAL install.
First, let’s load the packages we need. We’ll need rayshader
(of course) for 3D plotting, raster
for loading and manipulating the data, scales
to rescale the color channels to adjust image contrast, and sp
to transform some point coordinates between coordinate systems. FYI, I’m using the namespace operator ::
throughout the code just to be clear what packages I’m calling from, but you don’t need to do that after you’ve loaded the package with library()
.
# install.packages(c("rayshader", "raster", "sp"))
library(rayshader)
library(sp)
library(raster)
library(scales)
Let’s load our data. We’ll start with loading the elevation dataset, since it’s an easier process. We’ll use the raster::raster()
function to load our SRTM hgt
files. This gives us two raster objects. We then combine the two with the raster::merge()
function. Since they’re coming from the same data source, we don’t have to worry about transforming coordinate systems to match–they should just merge together seamlessly. We’ll plot the elevation using rayshader::height_shade()
so you can see what it looks like:
elevation1 = raster::raster("LC08_L1TP_038034_20191117_20191202_01_T1/N37W113.hgt")
elevation2 = raster::raster("~/Desktop/LC08_L1TP_038034_20191117_20191202_01_T1/N37W114.hgt")
zion_elevation = raster::merge(elevation1,elevation2)
height_shade(raster_to_matrix(zion_elevation)) %>%
plot_map()

height_shade()
function. Zion National Park is in the center.Great! Now, let’s load our georeferenced satellite imagery. The red, blue, and green bands on Landsat 8 data are bands B4, B3, and B2, respectively. You can delete the rest of the bands if you want to free 600 megabytes of space on your hard drive. Let’s then plot it with raster::plotRGB()
.
zion_r = raster::raster("LC08_L1TP_038034_20191117_20191202_01_T1_B4.TIF")
zion_g = raster::raster("LC08_L1TP_038034_20191117_20191202_01_T1_B3.TIF")
zion_b = raster::raster("LC08_L1TP_038034_20191117_20191202_01_T1_B2.TIF")
zion_rbg = raster::stack(zion_r, zion_g, zion_b)
raster::plotRGB(zion_rbg, scale=255^2)

Why is it so dark? We need to apply a gamma correction to the the imagery, which corrects raw linear intensity data for our non-linear perception of darkness. We do that simply by taking the square root of the data (we can also remove the scale argument).
zion_rbg_corrected = sqrt(raster::stack(zion_r, zion_g, zion_b))
raster::plotRGB(zion_rbg_corrected)

Great! The contrast is muted, but we’ll fix that later. Now, let’s combine the two datasets. Here we encounter our first GIS difficulty: the coordinate reference systems of our elevation data and imagery data don’t match! Oh no.
raster::crs(zion_r)
## CRS arguments:
## +proj=utm +zone=12 +datum=WGS84 +units=m +no_defs +ellps=WGS84
## +towgs84=0,0,0
raster::crs(zion_elevation)
## CRS arguments:
## +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
Our imagery data is given in UTM coordinates, while our elevation is in long/lat. Thankfully, it’s fairly straightforward to transform the elevation data from long/lat to UTM with the raster::projectRaster()
function. We use the “bilinear” method for interpolation since elevation is a continuous variable.
crs(zion_r)
## CRS arguments:
## +proj=utm +zone=12 +datum=WGS84 +units=m +no_defs +ellps=WGS84
## +towgs84=0,0,0
zion_elevation_utm = raster::projectRaster(zion_elevation, crs = crs(zion_r), method = "bilinear")
crs(zion_elevation_utm)
## CRS arguments:
## +proj=utm +zone=12 +datum=WGS84 +units=m +no_defs +ellps=WGS84
## +towgs84=0,0,0
Now, let’s crop the region down to the park itself. To figure out what long/lat values enclosed our area, I went to Google Maps and double clicked on the map to extract the bottom left and top right corners of the bounding box for the final map. We’ll use the sp
package to transform these long/lat coords into UTM coordinates.
bottom_left = c(y=-113.155277, x=37.116253)
top_right = c(y=-112.832502, x=37.414948)
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(zion_elevation_utm))
e = raster::extent(extent_utm)
e
## class : Extent
## xmin : 308511.9
## xmax : 337834.2
## ymin : 4109943
## ymax : 4142481
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.
zion_rgb_cropped = raster::crop(zion_rbg_corrected, e)
elevation_cropped = raster::crop(zion_elevation_utm, e)
names(zion_rgb_cropped) = c("r","g","b")
zion_r_cropped = rayshader::raster_to_matrix(zion_rgb_cropped$r)
zion_g_cropped = rayshader::raster_to_matrix(zion_rgb_cropped$g)
zion_b_cropped = rayshader::raster_to_matrix(zion_rgb_cropped$b)
zionel_matrix = rayshader::raster_to_matrix(elevation_cropped)
zion_rgb_array = array(0,dim=c(nrow(zion_r_cropped),ncol(zion_r_cropped),3))
zion_rgb_array[,,1] = zion_r_cropped/255 #Red layer
zion_rgb_array[,,2] = zion_g_cropped/255 #Blue layer
zion_rgb_array[,,3] = zion_b_cropped/255 #Green layer
zion_rgb_array = aperm(zion_rgb_array, c(2,1,3))
plot_map(zion_rgb_array)

We will also now scale our data to improve the contrast and make the image more vibrant. I’m going to use the scales
package to rescale the image to use the full range of color.
zion_rgb_contrast = scales::rescale(zion_rgb_array,to=c(0,1))
plot_map(zion_rgb_contrast)

Excellent. Now we just input this image into plot_3d()
along with our elevation data. For a realistic landscape, we should set zscale = 30
(since the elevation data was taken at 30 meter increments), but we’re going to set zscale = 15
to give the landscape 2x exaggeration. If you’re following along, you should now have an rgl
window open with a 3D model you can rotate around.
plot_3d(zion_rgb_contrast, zionel_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 = "Zion National Park, Utah | Imagery: Landsat 8 | DEM: 30m SRTM",
title_bar_color = "#1f5214", title_color = "white", title_bar_alpha = 1)

The static snapshot above is nice, but 3D is best seen in motion, so let’s create a movie of us rotating around the scene. We’ll move the camera with render_camera()
, and generate 1440 frames with render_snapshot()
. At 60 frames per second, this will generate a 24 second video. You can do this via the av
package, but I’m going to call ffmpeg
directly via a system()
call, since the output of av
doesn’t seem to play nicely with embeddable web videos (for some reason).
angles= seq(0,360,length.out = 1441)[-1]
for(i in 1:1440) {
render_camera(theta=-45+angles[i])
render_snapshot(filename = sprintf("zionpark%i.png", i),
title_text = "Zion National Park, Utah | Imagery: Landsat 8 | DEM: 30m SRTM",
title_bar_color = "#1f5214", title_color = "white", title_bar_alpha = 1)
}
rgl::rgl.close()
#av::av_encode_video(sprintf("zionpark%d.png",seq(1,1440,by=1)), framerate = 30,
# output = "zionpark.mp4")
rgl::rgl.close()
system("ffmpeg -framerate 60 -i zionpark%d.png -pix_fmt yuv420p zionpark.mp4")
And we’re done! What to do now? Well, this satellite imagery and elevation data is relatively low resolution–you could now search for higher resolution imagery and elevation data to produce a more detailed figure. The general process is the same for any data source: find the data, transform the data into a common coordinate system, crop the data to the desired region, and then combine the data in plot_3d()
. You could also try using the render_highquality()
function in rayshader
to make beautiful, pathtraced maps (like in the featured animation at the top of the page). Or you could learn more about using rayshader
(include transforming LIDAR point data to extremely hi-res elevation data) in my masterclass: the materials are free and open source. And this blog post is by no means the authoritative source on how to do this: there are lots of packages and free books on remote sensing and GIS. Have fun!
And whatever you end up doing, if you create something cool be sure to tweet it with the hashtags #rstats and #rayshader! We’d all love to see what you’ve created.
Tyler is a data scientist, physicist, writer, and programmer.
Hey Tyler,
Awesome content as always, one thing I am curious about is how accurately do both data sources need to line up? Like once we get the SRTM tile does the landsat imagery need to be the exact same coordinates and elevation in order for everything to turn out as awesome as yours? If so is there any easy way to determine the coordinates? It appears the SRTM zipfile contains one set of coordinates is there a way to use those to figure out a bounding box?
Thanks again for everything you do
Thanks! The
projectRaster
function lines up the coordinates between the two, and thecrop
function makes sure they cover the exact same extent. I’d recommend the process I showed where you figure out your bounding box in lat/long, and use thesp::spTransform
function to convert them to the SRTM coordinates.Hello,
Can I translate this work to Arabic language and share it in public because many student do not good in english with the credit for you?
Hi Tyler,
Thank you so much this illustration.
Sorry for my silly Q. I have lost to find “LC08_L1TP_038034_20191117_20191202_01_T1” used in for generating elevation1 and elevation2 at the beginning. I could not find it on downloads from Derek Watkin’s SRTM tile downloader and USGS Earth Explorer. What did I miss?
So many thanks
Hi Samu,
That should be available on the USGS Earth Explorer–if you follow my instructions exactly, it should be on the second or third page of search results for Landsat 8 imagery in that part of Utah.
Hi Tyler
Thank you.
I got all the files from file from USGS and Derek Watkin’s SRTM tile downloader except LC08_L1TP_038034_20191117_20191202_01_T1 file used to extract elevation1 and elevation2 in the demo. I missed to get where to get from this specific file. Sorry!
I elevation1 = raster::raster(“LC08_L1TP_038034_20191117_20191202_01_T1/N37W113.hgt”).
With kind regards
Hello Samu,
Assign the elevation1 and elevation2 variables to the file path where the .hgt is located.
For instance, on my computer:
elevation1 = raster::raster(“E:/Data/Rayshader_Tutorial/N37W113.SRTMGL1.hgt/N37W113.hgt”)
elevation2 = raster::raster(“E:/Data/Rayshader_Tutorial/N37W114.SRTMGL1.hgt/N37W114.hgt”)
I’m not sure how these were defined under the LC08 structure in the tutorial notes as LC08 is the file structure that USGS uses that designates the imagery as Landsat 8.
You should have success with the elevation data variable if you make this change.
Hi Mike,
Thank you. Figured out now assigning the assigning path to N37W113.hgt and N37W114.hgt files. I am using a R project structure and put all the files to be called there so that it can be accessed directly. In this case:
elevation1 = raster::raster(“N37W113.hgt”)
elevation2 = raster::raster(“N37W114.hgt”)
I got the 2D Figure 8 done (plot_map(zion_rgb_array)!
Sadly, it throws error for cropping and fails to make the 3D plot.
> elevation_cropped = raster::crop(zion_elevation_utm, e)
Error in .local(x, y, …) : extents do not overlap
I don’t what did go wrong. Anything I am missing?
Thank you again.
This is simply brilliant
[…] Tyler Morgan-Wall […]
Loving the tutorial, and got a 3d map, however my image is very whitewashed even after increasing contrast. Is there anything else I can do to help this?
I also tried adjusting it in Photoshop, but doing that deletes the coordinate system in the file, and I can not overlay that with the .hgt files.
Hi Tyler –
Great work! I enjoyed using your workflow to create a rotating model of nearby (to me) Mt. Rainier.
Have you had a chance to work with outputs from drone imagery, with resolution in the 1cm/px range? I’ve been trying to apply this workflow to the outputs of Open Drone Map (ODM – geotiff orthophotos and DSMs and DTMs), but haven’t been wholly successful yet. The .hgt files from your example work fine, but I get null-type results using the ODM DTMs .tif files.
Any help would be greatly appreciated! Thanks!
Hi Tyler,
I appear to be the first true novice commenting here. I would really like to use your tool to 3d print a couple of terrains that I have visited.
I followed everything necessary to get my tif file but that is where I had to stop. If possible could you right up a side tutorial or even point to one that gets me to where I can run the following on my maching? I have pulled your project down from github.
# To install the latest version from Github:
# install.packages(“devtools”)
devtools::install_github(“tylermorganwall/rayshader”)
I am thinking that there is a dependency or package manager that I need and don’t know about. I have read other peoples experiences using your tool but they all start their journey from the same spot that I don’t know how to get to.
Hi Nate,
Have you tried just straight up pull it down with:
install.packages(“rayshader”)
Worked for me. Cheers!
Hey Tyler, Just curious if you might know why I am getting a mismatching of the solid-fill when using plot_3d() function? It just looks like the solid fill underneath is pivoted off the origin of the elevation layer, and even using the solid = FALSE argument in plot_3d seems to have little to no effect on the rendered outcome.
my code can be seen here on GitHub:
https://github.com/samshott/Creating-a-3d-map-in-R/blob/master/R_3dMap
Can you link to an image showing the issue?
Hi, awesome tutorial. I am trying to add some point locations, but have failed to do so. My approach is to use extract() from the DEM and add the altitude (z) as a column to my csv wit xy data (lat,long). But I don´t seem to be able to ploit it. Any advice?
Hi, really cool post. I would like to know your input on plotting some points or a shp on top of the 3d map. I have being trying some stuff, but nothing seems to work. Thanks again for the tutorial!
Hi Reinaldo, check out the `generate_polygon_overlay()` and `render_points()` functions
This tutorial is fantastic – worked for my own area of interest on the first try! Thanks Tyler. Keep up the good work.
Hi, I got a problem with this:
Error in .local(x, y, …) : extents do not overlap
In adition, my final map in 3D is B/W, I can’t see any colour. Why? 🙁
Hi Sebastian, that error means it’s likely your elevation data and the imagery don’t share the same CRS. For the lack of color, try plotting the raster on it’s own and make sure it’s multiband
How do you add a scale bar and compas directions?
Cheers
Rayshader has both a `render_scalebar()` and a `render_compass()` function that do precisely those two things 🙂
https://www.rayshader.com/reference/render_scalebar.html
https://www.rayshader.com/reference/render_compass.html
This is a really good step to step guide! Any idea how could I stitch together a few images? my interested area happened to be in the corner of landset images. So no single image will cover the entire area, I have to use multiple images and combined them.
See the “Mapping the Effects of Rising Sea Levels” section in my masterclass (https://github.com/tylermorganwall/MusaMasterclass), there’s an example there of merging several rasters to make one combined raster. It’s for elevation data but imagery should follow the same process.
Thank you Tyler. I’ve been going through your guide and recreating your tutorial this evening. It works! I’m thinking of doing the area around Uluru (Ayers Rock) next.
Rendering is quite slow as it appears to be GPU bound. Would you know how I might be able to render using my GPU (CUDA) instead?
Hi Steve, the package that powers the 3D rendering (rgl) does in fact use the GPU, but building the actual 3D model can only be done in the CPU. Try setting “triangulate = TRUE” and play with the “max_error” argument and see if that speeds up the rendering process.
I meant CPU-bound
Hi, appreciate if you can have a look at why the texture of 3D model is missing after I apply an overlay. Below is my code:
data %>%
sphere_shade(texture = “desert”,zscale = 0.1) %>%
add_shadow(ray_shade(data, zscale = 0.1), 0.5) %>%
add_overlay(height_shade(data2, texture = topo.colors(256))) %>%
plot_3d(data, zscale = 0.1, fov = 0, theta = 135, zoom = 0.75, phi = 45, windowsize = c(1000, 800),background = “#F2E1D0”, shadowcolor = “#523E2B”)
render_highquality(lightdirection = c(-45,45), lightaltitude = 30)
Link to rendered image:
https://user-images.githubusercontent.com/8111532/99611011-bde04b00-2a77-11eb-96df-d9b11723112a.png
Thanks for the great tutorial!!
Hi,
can I used an FCC, or a landcover classified image instead of the RGB bands? These data are already cropped to the required extent.
Regards
Anjan
They just need to be converted to 3-layer RGB arrays to be used in rayshader (see the output of
sphere_shade()
for an example)Hi,
Thanks for the reply. I used the bands and reached this step
plot_3d(zion_rgb_contrast, zionel_matrix, windowsize = c(1100,900), zscale = 15, shadowdepth = -50,
zoom=0.5, phi=45,theta=-45,fov=70, background = “#F2E1D0”, shadowcolor = “#523E2B”)
However, for all subsequent steps, I get this error:
Error in loadNamespace(name) : there is no package called ‘rayimage’.
When i try to install rayimage, it does not install and in the end says that the installation of the package had a non-zero exit status.
Would you have a workaround to this?
Regards,
Anjan
solved.
reinstalled rayshader and it worked!
Thank for well explained sequential tutorial…
Is there a way to use a google or mapbox satellite layer as the base layer?
Regards,
Anjan
https://rdrr.io/github/poissonconsulting/poisspatial/src/R/basemap.R
^^^check out ggmap_to_raster()… returns a raster stack which you can transpose like the raster stack in the example. use crs(ras) <- 4326 just before the return of the function to add a projection to the map for cropping later..
Thanks for the great tutorial Tyler !
I followed the step-by-step process to render the Blue Mountains region in Sydney.
While the elevation looks apt, the satellite image overlay seems very hazy and low resolution.
https://mapping.s3.ap-south-1.amazonaws.com/Blue+mountains+hazy.JPG
I downloaded it from https://earthexplorer.usgs.gov/ and the B2, B4, B3 bands were all ~80MB in size.
Based on result I am getting, wondering if you have any hypothesis why the 3D render is not as clear as yours ?
Thanks !
Saurabh
This is stunning stuff. I am struggling to reproduce it with my own stuff; my resulting images look extremely narrow and tall. I wonder – do the elevation data and satellite imagery need to be at the same resolution? My elevation raster is 90 m; I don’t know what the USGS-EROS satellite data are, but seem to be higher res.
I am having trouble unzipping my files. They were downloaded as a GZ file instead of a ZIP file.
I use 7-zip to extract the tar .gz file.
Thanks Tyler – great tutorial! Followed all the steps successfully, except until the very last one where we create the movie.
I’m copying and pasting your precise code into R Studio and get the following error message:
Error in png::writePNG(hillshade, filename, text = c(source = “rayshader”)) :
unable to create zionpark1.png
Do you have any idea why please?
Thanks.
[kindly ignore above question. the issue was related to incorrectly setting my working directory :-)]
Awesome tutorial, thanks so much for everything you’ve done on this!
For some reason I’m having trouble adding labels to this map using render_label, but likely because there are many parts of the package that I don’t understand yet. I have no problem adding labels to the montereybay example in the documentation and I *am* able to add a compass to this map with no problem. Any insight into what my problem with labels might be? I’ve tried following along on a few different examples online and it seems many of them are using xs and ys from a previous version of rayshader (but I will keep looking).
Thanks again for a great package!
Hi Tyler, When I download the sat imagery file (GeoTiff), it unpacks as a single file. I’m expecting to see many files with the different bands as suggested, but I see only one file. I’ve tested the other files that are also available to download, but I don’t see files B2-B4.tiff band files. Am I doing something wrong?
My bad…I see than I need to run a 2-step extraction of that .gz file. Upon looking at the ‘archive’ , I see the Band files.
Whenever i run the elevation1 code, i get this error below, i cant seem to get passed this point
Error in .local(.Object, …) :
Error in .rasterObjectFromFile(x, band = band, objecttype = “RasterLayer”, :
Cannot create a RasterLayer object from this file. (file does not exist)
I get the error: Error in .local(x, y, …) : extents do not overlap.
Could you give me a help?
Thank you!
Hi Tyler,
Awesome tutorial, you do amazing stuff.
I have a question about your Lat/Long crop and conversion to UTM. I’m getting an odd output.
Is the Proj.4 representation depreciated?
For example:
Deprecated Proj.4 representation:
+proj=utm +zone=32 +datum=WGS84 +units=m +no_defs
Would this have to be updated to
WKT2 2019 representation:
???
Everything else in the tutorial seems to work great.
Hi Tyler, there is a way to put a shapefile above the 3d satelite/dem map? maybe a basin or city borders?
Hi Tyler – thank you for this post. It is very helpful and an excellent resource.
It would be great if labels could be placed on the maps, perhaps outside the map with stems connected to a point or region to be labelled… Also, I would love be able to change the data mapped on the DEM. For instance, instead of mapping RGB data, it would be great to be able to show trends in, say, precipitation from a univariate raster.
A minor thing – in some of your assigned variables, you use rbg and then rgb. Maybe you did it on purpose though… Thanks again for this great package and tutorial!