Devoted to open data and open source in science and education.

View All Tutorials

This tutorial is a part of a series!

Click below to view all lessons in the series!


LiDAR (9)
R programming (70)
Remote Sensing (12)
Data Visualization (4)
Hyperspectral Remote Sensing (7)
Hierarchical Data Formats (HDF5) (24)
Spatial Data & GIS (18)
Time Series (15)
Phenology (7)
Raster Data (8)
Vector Data (6)
Metadata (1)
Git & GitHub (6)
(1) (1) (1)

Tutorial by R Package

dplyr (8)
ggplot2 (17)
h5py (1)
lubridate (time series) (6)
maps (1)
maptools (3)
plyr (2)
raster (32)
rasterVis (raster time series) (3)
rgdal (GIS) (23)
rgeos (5)
rhdf5 (21)
sp (7)
scales (4)
gridExtra (4)
ggtheme (0)
grid (2)
reshape2 (3)
plotly (6)

View ALL Tutorial Series

Twitter Youtube Github


R Bloggers


In this tutorial, we will learn how to process data that cover different spatial extents in R. Note that this assumes the data are on the same raster grid and same resolution.

First, let’s load the required libraries.

# load libraries

# be sure to set your working directory
# setwd("~/Documents/data/NEONDI-2016") # Mac
# setwd("~/data/NEONDI-2016")  # Windows

## import functions
# install devtools (only if you have not previously intalled it)
# call devtools library

# install from github
# call library

# source("/Users/lwasser/Documents/GitHub/neon-aop-package/neonAOP/R/aop-data.R")

Import LiDAR data

To begin, let’s open the NEON LiDAR Digital Surface and Digital Terrain Models (DSM and DTM) which are in GeoTIFF format.

# read aspect data from previous lesson
TEAK_nsAspect <- raster("outputs/TEAK/TEAK_nsAspect.tif")

Mask Data

Once we have created a threhold classified raster, we can use it for different things. One application is to use it as an analysis mask for another dataset.

Let’s try to find all pixels that have an NDVI value >.6 and are north facing.

# Define the file name to be opened
f <- "NEONdata/D17-California/TEAK/2013/spectrometer/reflectance/Subset3NIS1_20130614_100459_atmcor.h5"

# define the CRS in EPGS format for the file
epsg <- 32611
# create a list of bands
bands <- c(60,83)

ndvi.stack <- create_stack(f, 

# calculate NDVI
ndvi <- (ndvi.stack[[2]]-ndvi.stack[[1]]) / (ndvi.stack[[2]]+ndvi.stack[[1]])
names(ndvi) <- "TEAK_hsiNDVI"

# let's test this out

# let's create a mask
ndvi[ndvi<.6] <- NA

# force the two to have the same CRS
crs(ndvi) <- crs(TEAK_nsAspect)

Create Raster Stack

Next, let’s create a raster stack of ndvi and aspect.

new.stack <- stack(TEAK_nsAspect, ndvi)

## Error in compareRaster(x): different extent

Notice we get an error. Why?

Let’s compare the extents of the two objects.

# view extents of both objects

## class       : Extent 
## xmin        : 325963 
## xmax        : 326507 
## ymin        : 4102904 
## ymax        : 4103482


## class       : Extent 
## xmin        : 325963 
## xmax        : 326506 
## ymin        : 4102905 
## ymax        : 4103482

# are the extents the same?
extent(ndvi) == extent(TEAK_nsAspect)

## [1] FALSE

The extents are slightly different. They are one pixel apart in ymin and xmax. Thus, when we try to create a stack, we get an error. All layers in a stack need to be of the same extent.

We can create an if statement that checks the extent and crops both rasters to the overlap. Let’s try it.

# check the extents of the two layers -- if they are different
# crop both datasets
if (extent(ndvi) == extent(TEAK_nsAspect)){
  print("Extents are the same, no need to crop")
  } else {
  # calculate overlap between the two datasets
  overlap <- intersect(extent(ndvi), extent(TEAK_nsAspect))
  # now let's crop both datasets to the overlap region
  ndvi <- crop(ndvi, overlap)
  asp.ns <- crop(TEAK_nsAspect, overlap)
  print("Extents are different, cropping data")

## [1] "Extents are different, cropping data"

# let's try to create a stack again.
new.stack <- stack(TEAK_nsAspect, ndvi)

# mask out only pixels that are north facing and NDVI >.6
nsFacing.ndvi <- mask(new.stack[[1]], new.stack[[2]])
nsFacing.ndvi[nsFacing.ndvi==0] <- NA

Create Final Plot

# plot extent
plot.extent <- extent(nsFacing.ndvi)

# plot 
     main="North & South Facing pixels, NDVI > .6",

# allow legend to plot outside of bounds

legend((par()$usr[2] + 20), plot.extent@ymax-200, # set x,y legend location
       legend = c("North", "South"),
       fill = c("blue", "green"), 
       bty="n") # turn off border

Export Classified Raster

# export geotiff 
            overwrite = TRUE,
            NAflag = -9999)

Get Lesson Code

(some browsers may require you to right click.)