Plot Raster Data
Overview
Teaching: 25 min
Exercises: 15 minQuestions
How can I create categorized or customized maps of raster data?
How can I customize the color scheme of a raster image?
How can I layer raster data in a single image?
Objectives
Build customized plots for a single band raster using the
ggplot2
package.Layer a raster dataset on top of a hillshade to create an elegant basemap.
Things You’ll Need To Complete This Episode
See the lesson homepage for detailed information about the software, data, and other prerequisites you will need to work through the examples in this episode.
Plot Raster Data in R
This episode covers how to plot a raster in R using the ggplot2
package with customized coloring schemes.
It also covers how to layer a raster on top of a hillshade to produce
an eloquent map. We will continue working with the Digital Surface Model (DSM) raster
for the NEON Harvard Forest Field Site.
At the end of this lesson, we will have made this map:
Plotting Data Using Breaks
In the previous episode, we viewed our data using a continuous color ramp. For
clarity and visibility of the plot, we may prefer to view the data “symbolized” or colored according to ranges of values. This is comparable to a “classified”
map. To do this, we need to tell ggplot
how many groups to break our data into, and
where those breaks should be. To make these decisions, it is useful to first explore the distribution of the data using a bar plot. To begin with, we will use dplyr
’s mutate()
function combined with cut()
to split the data into 3 bins.
HARV_DSM_df <- HARV_DSM_df %>%
mutate(fct_elevation = cut(Altitude, breaks = 3))
ggplot() +
geom_bar(data = HARV_DSM_df, aes(fct_elevation))
If we want to know the cutoff values for the groups, we can ask for the unique values
of fct_elevation
:
unique(HARV_DSM_df$fct_elevation)
[1] (379,416] (342,379] (305,342]
Levels: (305,342] (342,379] (379,416]
And we can get the count of values in each group using dplyr
’s
group_by()
and count()
functions:
HARV_DSM_df %>%
group_by(fct_elevation) %>%
count()
# A tibble: 3 × 2
# Groups: fct_elevation [3]
fct_elevation n
<fct> <int>
1 (305,342] 418891
2 (342,379] 1530073
3 (379,416] 370835
We might prefer to customize the cutoff values for these groups.
Lets round the cutoff values so that we have groups for the ranges of
301–350 m, 351–400 m, and 401–450 m.
To implement this we will give mutate()
a numeric vector of break points instead
of the number of breaks we want.
custom_bins <- c(300, 350, 400, 450)
HARV_DSM_df <- HARV_DSM_df %>%
mutate(fct_elevation_2 = cut(Altitude, breaks = custom_bins))
unique(HARV_DSM_df$fct_elevation_2)
[1] (400,450] (350,400] (300,350]
Levels: (300,350] (350,400] (400,450]
Data Tips
Note that when we assign break values a set of 4 values will result in 3 bins of data.
The bin intervals are shown using
(
to mean exclusive and]
to mean inclusive. For example:(305, 342]
means “from 306 through 342”.
And now we can plot our bar plot again, using the new groups:
ggplot() +
geom_bar(data = HARV_DSM_df, aes(fct_elevation_2))
And we can get the count of values in each group in the same way we did before:
HARV_DSM_df %>%
group_by(fct_elevation_2) %>%
count()
# A tibble: 3 × 2
# Groups: fct_elevation_2 [3]
fct_elevation_2 n
<fct> <int>
1 (300,350] 741815
2 (350,400] 1567316
3 (400,450] 10668
We can use those groups to plot our raster data, with each group being a different color:
ggplot() +
geom_raster(data = HARV_DSM_df , aes(x = x, y = y, fill = fct_elevation_2)) +
coord_quickmap()
The plot above uses the default colors inside ggplot
for raster objects.
We can specify our own colors to make the plot look a little nicer.
R has a built in set of colors for plotting terrain, which are built in
to the terrain.colors()
function.
Since we have three bins, we want to create a 3-color palette:
terrain.colors(3)
[1] "#00A600" "#ECB176" "#F2F2F2"
The terrain.colors()
function returns hex colors -
each of these character strings represents a color.
To use these in our map, we pass them across using the
scale_fill_manual()
function.
ggplot() +
geom_raster(data = HARV_DSM_df , aes(x = x, y = y,
fill = fct_elevation_2)) +
scale_fill_manual(values = terrain.colors(3)) +
coord_quickmap()
More Plot Formatting
If we need to create multiple plots using the same color palette, we can create
an R object (my_col
) for the set of colors that we want to use. We can then
quickly change the palette across all plots by modifying the my_col
object, rather than each individual plot.
We can label the x- and y-axes of our plot too using xlab
and ylab
.
We can also give the legend a more meaningful title by passing a value
to the name
argument of the scale_fill_manual()
function.
my_col <- terrain.colors(3)
ggplot() +
geom_raster(data = HARV_DSM_df , aes(x = x, y = y,
fill = fct_elevation_2)) +
scale_fill_manual(values = my_col, name = "Elevation") +
coord_quickmap()
Or we can also turn off the labels of both axes by passing element_blank()
to
the relevant part of the theme()
function.
ggplot() +
geom_raster(data = HARV_DSM_df , aes(x = x, y = y,
fill = fct_elevation_2)) +
scale_fill_manual(values = my_col, name = "Elevation") +
theme(axis.title = element_blank()) +
coord_quickmap()
Challenge: Plot Using Custom Breaks
Create a plot of the Harvard Forest Digital Surface Model (DSM) that has:
- Six classified ranges of values (break points) that are evenly divided among the range of pixel values.
- Axis labels.
- A plot title.
Answers
HARV_DSM_df <- HARV_DSM_df %>% mutate(fct_elevation_6 = cut(Altitude, breaks = 6)) my_col <- terrain.colors(6) ggplot() + geom_raster(data = HARV_DSM_df , aes(x = x, y = y, fill = fct_elevation_6)) + scale_fill_manual(values = my_col, name = "Elevation") + ggtitle("Classified Elevation Map - NEON Harvard Forest Field Site") + xlab("UTM Easting Coordinate (m)") + ylab("UTM Northing Coordinate (m)") + coord_quickmap()
Layering Rasters
We can layer a raster on top of a hillshade raster for the same area, and use a transparency factor to create a 3-dimensional shaded effect. A hillshade is a raster that maps the shadows and texture that you would see from above when viewing terrain. In our example, the texture will actually be treetops.
We will add a custom color, making the plot grey.
First we need to read in our DSM hillshade data and view the structure:
HARV_hill <-
raster("data/NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_DSMhill.tif")
HARV_hill
class : RasterLayer
dimensions : 1367, 1697, 2319799 (nrow, ncol, ncell)
resolution : 1, 1 (x, y)
extent : 731453, 733150, 4712471, 4713838 (xmin, xmax, ymin, ymax)
crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
source : HARV_DSMhill.tif
names : HARV_DSMhill
values : -0.7136298, 0.9999997 (min, max)
Next we convert it to a dataframe, so that we can plot it using ggplot2
:
HARV_hill_df <- as.data.frame(HARV_hill, xy = TRUE)
str(HARV_hill_df)
'data.frame': 2319799 obs. of 3 variables:
$ x : num 731454 731454 731456 731456 731458 ...
$ y : num 4713838 4713838 4713838 4713838 4713838 ...
$ HARV_DSMhill: num NA NA NA NA NA NA NA NA NA NA ...
Now we can plot the hillshade data:
ggplot() +
geom_raster(data = HARV_hill_df,
aes(x = x, y = y, alpha = HARV_DSMhill)) +
scale_alpha(range = c(0.15, 0.65), guide = "none") +
coord_quickmap()
Data Tips
Turn off, or hide, the legend on a plot by adding
guide = "none"
to ascale_something()
function or by settingtheme(legend.position = "none")
.The alpha value determines how transparent the colors will be (0 being transparent, 1 being opaque).
We can layer another raster on top of our hillshade by adding another call to
the geom_raster()
function. Let’s overlay HARV_DSM
on top of the hill_HARV
.
ggplot() +
geom_raster(data = HARV_DSM_df ,
aes(x = x, y = y,
fill = Altitude)) +
geom_raster(data = HARV_hill_df,
aes(x = x, y = y,
alpha = HARV_DSMhill)) +
scale_fill_viridis_c() +
scale_alpha(range = c(0.15, 0.65), guide = "none") +
ggtitle("Elevation with hillshade") +
coord_quickmap()
Challenge: Create DTM & DSM for San Joaquin Ecological Reserve
Use the files in the
data/NEON-DS-Airborne-Remote-Sensing/SJER/
directory to create a Digital Terrain Model map and Digital Surface Model map of the San Joaquin Experimental Range field site.Make sure to:
- include hillshade in the maps,
- label axes on the DSM map and exclude them from the DTM map,
- include a title for each map,
- experiment with various alpha values and color palettes to represent the data.
Answers
# CREATE DSM MAPS # import DSM data SJER_DSM <- raster("data/NEON-DS-Airborne-Remote-Sensing/SJER/DSM/SJER_dsmCrop.tif") # convert to a df for plotting SJER_DSM_df <- as.data.frame(SJER_DSM, xy = TRUE) # import DSM hillshade SJER_DSM_hill <- raster("data/NEON-DS-Airborne-Remote-Sensing/SJER/DSM/SJER_dsmHill.tif") # convert to a df for plotting SJER_DSM_hill_df <- as.data.frame(SJER_DSM_hill, xy = TRUE) # Build Plot ggplot() + geom_raster(data = SJER_DSM_df , aes(x = x, y = y, fill = SJER_dsmCrop, alpha = 0.8) ) + geom_raster(data = SJER_DSM_hill_df, aes(x = x, y = y, alpha = SJER_dsmHill) ) + scale_fill_viridis_c() + guides(fill = guide_colorbar()) + scale_alpha(range = c(0.4, 0.7), guide = "none") + # remove grey background and grid lines theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + xlab("UTM Easting Coordinate (m)") + ylab("UTM Northing Coordinate (m)") + ggtitle("DSM with Hillshade") + coord_quickmap()
# CREATE DTM MAP # import DTM SJER_DTM <- raster("data/NEON-DS-Airborne-Remote-Sensing/SJER/DTM/SJER_dtmCrop.tif") SJER_DTM_df <- as.data.frame(SJER_DTM, xy = TRUE) # DTM Hillshade DTM_hill_SJER <- raster("data/NEON-DS-Airborne-Remote-Sensing/SJER/DTM/SJER_dtmHill.tif") DTM_hill_SJER_df <- as.data.frame(DTM_hill_SJER, xy = TRUE) ggplot() + geom_raster(data = SJER_DTM_df , aes(x = x, y = y, fill = SJER_dtmCrop, alpha = 2.0) ) + geom_raster(data = DTM_hill_SJER_df, aes(x = x, y = y, alpha = SJER_dtmHill) ) + scale_fill_viridis_c() + guides(fill = guide_colorbar()) + scale_alpha(range = c(0.4, 0.7), guide = "none") + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + theme(axis.title.x = element_blank(), axis.title.y = element_blank()) + ggtitle("DTM with Hillshade") + coord_quickmap()
Key Points
Continuous data ranges can be grouped into categories using
mutate()
andcut()
.Use built-in
terrain.colors()
or set your preferred color scheme manually.Layer rasters on top of one another by using the
alpha
aesthetic.