Datacarpentries 01

Introduction to Geospatial Raster and Vector Data with R

Markdown Author: Jessie Bell

I. Raster Data

#first, preview the metadata. It is ideal to do this before importing your data.

beforereadingdata <- describe("data/NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_dsmCrop.tif")


#you can store the information above from the tif file by using the next code

HARV_dsmCrop_info <- capture.output(describe("data/NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_dsmCrop.tif"))

Metadata

A: Nomenclature

Use the following nomenclature for naming your strings, etc…

filetype _ location

• File Type = Digital Surface Model (DSM)
• Location = Harvard Forest

 → DSM_HARV

B. Summaries

Use summary(raster) to get data range statistics. This function will only summarize a random 100,000 cells from the raster though. To get all cells in summary stat, use the following functions together:

summary(values(DSM_HARV)

# read in the raster using rast() function

#naming below is datatype_HARV. HARV is because the data is from Harvard Forest. The DSM because the data filetype is DSM

DSM_HARV <- rast("data/NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_dsmCrop.tif")

DSM_HARV

## class       : SpatRaster 
## dimensions  : 1367, 1697, 1  (nrow, ncol, nlyr)
## resolution  : 1, 1  (x, y)
## extent      : 731453, 733150, 4712471, 4713838  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 18N (EPSG:32618) 
## source      : HARV_dsmCrop.tif 
## name        : HARV_dsmCrop 
## min value   :       305.07 
## max value   :       416.07

Visualize

A: Terra Package

You can use the terra package to convert tif to dataframe. Once you do this conversion, your map will be plotable.

DSM_HARV_df <- as.data.frame(DSM_HARV, xy = T)

Now you can view the structure of the data and see that it is formatted as a dataframe.

str(DSM_HARV_df)

## 'data.frame':    2319799 obs. of  3 variables:
##  $ x           : num  731454 731455 731456 731457 731458 ...
##  $ y           : num  4713838 4713838 4713838 4713838 4713838 ...
##  $ HARV_dsmCrop: num  409 408 407 407 409 ...

B: ggplot()

Now we can use ggplot to view the data.

ggplot() +
    geom_raster(data = DSM_HARV_df , aes(x = x, y = y, fill = HARV_dsmCrop)) +
    scale_fill_viridis_c() +
    coord_quickmap()

C: Terra’s plot function

plot(DSM_HARV)

CRS

A: Check out project string

We can view the CRS string associated with our R object using the crs() function.

crs(DSM_HARV, proj = T)

## [1] "+proj=utm +zone=18 +datum=WGS84 +units=m +no_defs"

#this project string has a zone - so we automatically know it uses a UTM 

#our units are in meters

Min & Max

It is useful to know the minimum or maximum values of a raster dataset. In this case, given we are working with elevation data, these values represent the min/max elevation range at our site.

Raster statistics are often calculated and embedded in a GeoTIFF for us. We can view these values:

minmax(DSM_HARV)

##     HARV_dsmCrop
## min       305.07
## max       416.07

If the min and max have not been calulcated yet, we can do that by using the following code:

DSM_HARV <- setMinMax(DSM_HARV)

We can see that the elevation at our site ranges from 305.0700073 m to 416.0699768 m.

Bands

A: How many bands?

To check how many bands your raster dataset contains, use the nly() function.

nlyr(DSM_HARV)

## [1] 1

B: Multi-band

Raster data can also be multi-band, meaning that one raster file contains data for more than one variable or time period for each cell. By default the raster() function only imports the first band in a raster regardless of whether it has one or more bands. Jump to a later episode in this series for information on working with multi-band rasters: “Work with Multi-band Rasters in R”.

Missing Data

To highlight NA values in ggplot, alter the scale_fill_*() layer to contain a colour instruction for NA values, like scale_fill_viridis_c(na.value = ‘deeppink’).

Find what the no data value is within a tif by using the following code:

moredescribe <- describe("data/NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_dsmCrop.tif")

showfilesource <- sources(DSM_HARV)

Output of describe prints things like:

[1] “Driver: GTiff/GeoTIFF”
[2] “Files: data/NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_dsmCrop.tif”

[3] “Size is 1697, 1367”
[4] “Coordinate System is:”
[5] “PROJCRS["WGS 84 / UTM zone 18N",”

[6] ”

Bad Data

Sometimes a raster’s metadata will tell us the range of expected values for a raster. Values outside of this range are suspect and we need to consider that when we analyze the data. Sometimes, we need to use some common sense and scientific insight as we examine the data - just as we would for field data to identify questionable values.

Plotting data with appropriate highlighting can help reveal patterns in bad values and may suggest a solution.

Histograms

We can explore the distribution of values contained within our raster using the geom_histogram() function which produces a histogram. Histograms are often useful in identifying outliers and bad data values in our raster data.

ggplot() +
    geom_histogram(data = DSM_HARV_df, aes(HARV_dsmCrop), bins = 40)

Note that the shape of this histogram looks similar to the previous one that was created using the default of 30 bins. The distribution of elevation values for our Digital Surface Model (DSM) looks reasonable. It is likely there are no bad data values in this particular raster.

practice <- describe("data/NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_DSMhill.tif")

II. Plotting Rasters

Setting Breaks

A: Using Data 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.

DSM_HARV_df_mutated <- DSM_HARV_df %>%
  mutate(fct_elevation = cut(HARV_dsmCrop, breaks = 6))


ggplot(data = DSM_HARV_df_mutated, aes(fct_elevation, fill=fct_elevation)) +
    geom_bar(show.legend = F)

plot(DSM_HARV_df_mutated$HARV_dsmCrop)

If we want to know the cutoff values for the groups, we can ask for the unique values of fct_elevation:

unique(DSM_HARV_df_mutated$fct_elevation)

## [1] (398,416] (379,398] (361,379] (342,361] (324,342] (305,324]
## Levels: (305,324] (324,342] (342,361] (361,379] (379,398] (398,416]

And we can get the count of values in each group using dplyr’s group_by() and count() functions:

DSM_HARV_df_mutated %>%
  group_by(fct_elevation) %>%
  count()

## # A tibble: 6 × 2
## # Groups:   fct_elevation [6]
##   fct_elevation      n
##   <fct>          <int>
## 1 (305,324]      30211
## 2 (324,342]     388680
## 3 (342,361]     777816
## 4 (361,379]     752257
## 5 (379,398]     349651
## 6 (398,416]      21184

Make it colorful!

Use groups with different colors to make cute maps!

ggplot()+
  geom_raster(data = DSM_HARV_df_mutated, aes(x, y, fill = fct_elevation))+
  coord_quickmap()

DSM_HARV_df_mutated <- DSM_HARV_df %>%
                mutate(fct_elevation = cut(HARV_dsmCrop, breaks = 10))


ggplot(data = DSM_HARV_df_mutated, aes(fct_elevation, fill=fct_elevation)) +
    geom_bar(show.legend = F)

ggplot()+
  geom_raster(data = DSM_HARV_df_mutated, aes(x, y, fill = fct_elevation))+
  coord_quickmap()+
  theme(legend.position = "bottom")

The plots above use ggplots default colors – we can use our own color palletes

my_col <- terrain.colors(10)

ggplot() +
 geom_raster(data = DSM_HARV_df_mutated , aes(x = x, y = y,
                                      fill = fct_elevation)) + 
    scale_fill_manual(values = terrain.colors(10)) + 
    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 = DSM_HARV_df_mutated , aes(x = x, y = y,
                                      fill = fct_elevation)) + 
    scale_fill_manual(values = my_col, name = "Elevation") +
    theme(axis.title = element_blank()) + 
    coord_quickmap()

Axes Labels

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 = DSM_HARV_df_mutated , aes(x = x, y = y,
                                      fill = fct_elevation)) + 
    scale_fill_manual(values = my_col, name = "Elevation") +
    theme(axis.title = element_blank()) + 
    coord_quickmap()

Challenge

Challenge! Using Custom Breaks

1. Six classified ranges of values (break points) that are evenly divided among the range of pixel values.

2. Axis labels.

3. A plot title.

DSM_HARV_df_mutated <- DSM_HARV_df %>%
                mutate(fct_elevation = cut(HARV_dsmCrop, breaks = 6))

my_col <- terrain.colors(6)

ggplot() + 
  geom_raster(data = DSM_HARV_df_mutated, aes(x = x, y = y, fill = fct_elevation)) + 
    scale_fill_manual(values = my_col, name = "Elevation") +
  xlab("UTM Easting Coordinate (m)")+
  ylab("UTM Northing Coordinate (m)")+
  ggtitle("Classified Elevation Map - NEON Harvard Forest Field Site")+
    theme() + 
    coord_quickmap()

III. Layering Rasters

Adding Hillshade

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. We will add a custom color, making the plot grey. |

DSM_hill_HARV <- rast('data/NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_DSMhill.tif')


DSM_hill_HARV

## class       : SpatRaster 
## dimensions  : 1367, 1697, 1  (nrow, ncol, nlyr)
## resolution  : 1, 1  (x, y)
## extent      : 731453, 733150, 4712471, 4713838  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 18N (EPSG:32618) 
## source      : HARV_DSMhill.tif 
## name        : HARV_DSMhill 
## min value   :   -0.7136298 
## max value   :    0.9999997

Convert to DF

DSM_hill_HARV_df <- as.data.frame(DSM_hill_HARV, xy = TRUE) 

str(DSM_hill_HARV_df)

## 'data.frame':    2313675 obs. of  3 variables:
##  $ x           : num  731455 731456 731457 731458 731459 ...
##  $ y           : num  4713837 4713837 4713837 4713837 4713837 ...
##  $ HARV_DSMhill: num  -0.15567 0.00743 0.86989 0.9791 0.96283 ...

Plot Hillshade

ggplot() +
  geom_raster(data = DSM_hill_HARV_df, 
              aes(x=x, y=y, alpha=HARV_DSMhill))+
      scale_alpha(range=c(0.15, 0.65), 
                  guide = "none")+
  coord_quickmap()

Layering

ggplot()+
  geom_raster(data=DSM_HARV_df, 
              aes(x=x, y=y, 
                  fill= HARV_dsmCrop))+
  geom_raster(data = DSM_hill_HARV_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

Use the files in the data/NEON-DS-Airborne-Remote-Sensing/SJER/ directory to create BOTH 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

DSM

#DSM Data

DSM_SJER <- rast("data/NEON-DS-Airborne-Remote-Sensing/SJER/DSM/SJER_dsmCrop.tif")

R_df <- as.data.frame(DSM_SJER, xy=T)

# DSM Hillshade

DSM_hill_SJER <- 
  rast("data/NEON-DS-Airborne-Remote-Sensing/SJER/DSM/SJER_dsmHill.tif")

DSM_hill_SJER_df <- as.data.frame(DSM_hill_SJER, xy=T)

#plot

ggplot()+
  geom_raster(data=R_df, aes(x, y, fill=SJER_dsmCrop, alpha=0.8))+
  geom_raster(data=DSM_hill_SJER_df, aes(x, y, alpha = SJER_dsmHill))+
  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())+
  xlab("UTM Easting Coordinate (m)")+
  ylab("UTM Northing Coordinate (m)")+
  ggtitle("DSM with Hillshade")

  coord_quickmap()

## <ggproto object: Class CoordQuickmap, CoordCartesian, Coord, gg>
##     aspect: function
##     backtransform_range: function
##     clip: on
##     default: FALSE
##     distance: function
##     expand: TRUE
##     is_free: function
##     is_linear: function
##     labels: function
##     limits: list
##     modify_scales: function
##     range: function
##     render_axis_h: function
##     render_axis_v: function
##     render_bg: function
##     render_fg: function
##     setup_data: function
##     setup_layout: function
##     setup_panel_guides: function
##     setup_panel_params: function
##     setup_params: function
##     train_panel_guides: function
##     transform: function
##     super:  <ggproto object: Class CoordQuickmap, CoordCartesian, Coord, gg>

DTM

DTM_SJER <- rast("data/NEON-DS-Airborne-Remote-Sensing/SJER/DTM/SJER_dtmCrop.tif")

DTM_hill_SJER <- rast("data/NEON-DS-Airborne-Remote-Sensing/SJER/DTM/SJER_dtmHill.tif")

#dataframe em

DTM_hill_SJER_df <- as.data.frame(DTM_hill_SJER, xy=T)

DTM_SJER_df <- as.data.frame(DTM_SJER, xy=T)

#ggplot it

ggplot()+
  geom_raster(data=DTM_SJER_df, aes(x, y, fill=SJER_dtmCrop, alpha=2.0))+
  geom_raster(data=DTM_hill_SJER_df, aes(x, y, alpha=SJER_dtmHill))+
  xlab("")+
  ylab("")+
  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())+
  ggtitle("DTM with Hillshade")+
  coord_quickmap()

Key Points

1. Continuous data ranges can be grouped into categories using mutate() and cut().

2. Use built-in terrain.colors() or set your preferred color scheme manually.

3. Layer rasters on top of one another by using the alpha aesthetic.