This tutorial covers how to clone a github.com repo to your computer so that you can work locally on files within the repo.

Clone - Copy Repo To Your Computer

In the previous tutorial, we used the github.com interface to fork the central NEON repo. By forking the NEON repo, we created a copy of it in our github.com account.

Now we will learn how to create a local version of our forked repo on our laptop, so that we can efficiently add to and edit repo content.

Copy Repo URL

Start from the github.com interface:

1. Navigate to the repo that you want to clone (copy) to your computer -- this should be YOUR-USER-NAME/DI-NEON-participants.
2. Click on the Clone or Download dropdown button and copy the URL of the repo.

Then on your local computer:

1. Your computer should already be setup with Git and a bash shell interface. If not, please refer to the Institute setup materials before continuing.
2. Open bash on your computer and navigate to the local GitHub directory that you created using the Set-up Materials.

To do this, at the command prompt, type:

$ cd ~/Documents/GitHub

Note: If you have stored your GitHub directory in a location that is different

• i.e. it is not /Documents/GitHub, be sure to adjust the above code to represent the actual path to the GitHub directory on your computer.

Now use git clone to clone, or create a copy of, the entire repo in the GitHub directory on your computer.

# clone the forked repo to our computer
$ git clone https://github.com/neon/DI-NEON-participants.git

The output shows you what is being cloned to your computer.

Cloning into 'DI-NEON-participants.git'...
remote: Counting objects: 3808, done.
remote: Total 3808 (delta 0), reused 0 (delta 0), pack-reused 3808
Receiving objects: 100% (3808/3808), 2.92 MiB | 2.17 MiB/s, done.
Resolving deltas: 100% (2185/2185), done.
Checking connectivity... done.
$

Note: The output numbers that you see on your computer, representing the total file size, etc, may differ from the example provided above.

View the New Repo

Next, let's make sure the repository is created on your computer in the location where you think it is.

At the command line, type ls to list the contents of the current directory.

# view directory contents
$ ls

Next, navigate to your copy of the data institute repo using cd or change directory:

# navigate to the NEON participants repository
$ cd DI-NEON-participants

# view repository contents
$ ls

404.md			_includes		code
ISSUE_TEMPLATE.md	_layouts		images
README.md		_posts			index.md
_config.yml		_site			institute-materials
_data			assets			org

Alternatively, we can view the local repo DI-NEON-participants in a finder (Mac) or Windows Explorer (Windows) window. Simply open your Documents in a window and navigate to the new local repo.

Using either method, we can see that the file structure of our cloned repo exactly mirrors the file structure of our forked GitHub repo.

Summary Workflow -- Create a Local Repo

In the github.com interface:

• Copy URL of the repo you want to work on locally

In shell:

• git clone URLhere

Note: that you can copy the URL of your repository directly from GitHub.

Got questions? No problem. Leave your question in the comment box below. It's likely some of your colleagues have the same question, too! And also likely someone else knows the answer.

In this tutorial, we will fork, or create a copy in your github.com account, an existing GitHub repository. We will also explore the github.com interface.

Create An Account

If you do not already have a GitHub account, go to GitHub and sign up for your free account. Pick a username that you like! This username is what your colleagues will see as you work with them in GitHub and Git.

Take a minute to setup your account. If you want to make your account more recognizable, be sure to add a profile picture to your account!

If you already have a GitHub account, simply sign in.

Navigate GitHub

Repositories, AKA Repos

Let's first discuss the repository or "repo". (The cool kids say repo, so we will jump on the git cool kid bandwagon) and use "repo" from here on in. According to the GitHub glossary:

A repository is the most basic element of GitHub. They're easiest to imagine as a project's folder. A repository contains all of the project files (including documentation), and stores each file's revision history. Repositories can have multiple collaborators and can be either public or private.

In the Data Institute, we will share our work in the DI-NEON-participants repo.

Find an Existing Repo

The first thing that you'll need to do is find the DI-NEON-participants repo. You can find repos in two ways:

1. Type “DI-NEON-participants” in the github.com search bar to find the repository.
2. Use the repository URL if you have it - like so: https://github.com/NEONScience/DI-NEON-participants.

Navigation of a Repo Page

Once you have found the Data Institute participants repo, take 5 minutes to explore it.

Git Repo Names

First, get to know the repository naming convention. Repository names all take the format:

OrganizationName/RepositoryName

So the full name of our repository is:

NEONScience/DI-NEON-participants

Header Tabs

At the top of the page you'll notice a series of tabs. Please focus on the following 3 for now:

• Code: Click here to view structure & contents of the repo.
• Issues: Submit discussion topics, or problems that you are having with the content in the repo, here.
• Pull Requests: Submit changes to the repo for review / acceptance. We will explore pull requests more in the Git 06 tutorial.

Other Text Links

A bit further down the page, you'll notice a few other links:

• commits: a commit is a saved and documented change to the content or structure of the repo. The commit history contains all changes that have been made to that repo. We will discuss commits more in Git 05: Git Add Changes -- Commits .

Fork a Repository

Next, let's discuss the concept of a fork on the github.com site. A fork is a copy of the repo that you create in your account. You can fork any repo at any time by clicking the fork button in the upper right hand corner on github.com.

Check Out Your Data Institute Fork

Now, check out your new fork. Its name should be:

YOUR-USER-NAME/DI-NEON-participants.

It can get confusing sometimes moving between a central repo:

• https://github.com/NEONScience/DI-NEON-participants

and your forked repo:

• https://github.com/YOUR-USER-NAME/DI-NEON-participants

A good way to figure out which repo you are viewing is to look at the name of the repo. Does it contain your username? Or your colleagues'? Or NEON's?

Your Fork vs the Central Repo

Your fork is an exact copy, or completely in sync with, the NEON central repo. You could confirm this by comparing your fork to the NEON central repository using the pull request option. We will learn about pull requests in Git06: Sync GitHub Repos with Pull Requests. For now, take our word for it.

The fork will remain in sync with the NEON central repo until:

1. You begin to make changes to your forked copy of the repo.
2. The central repository is changed or updated by a collaborator.

If you make changes to your forked repo, the changes will not be added to the NEON central repo until you sync your fork with the NEON central repo.

Summary Workflow -- Fork a GitHub Repository

On the github.com website:

• Navigate to desired repo that you want to fork.
• Click Fork button.

Have questions? No problem. Leave your question in the comment box below. It's likely some of your colleagues have the same question, too! And also likely someone else knows the answer.

In this page, you will be introduced to the importance of version control in scientific workflows.

The text and graphics in the first three sections were borrowed, with some modifications, from Software Carpentry's Version Control with Git lessons.

What is Version Control?

A version control system maintains a record of changes to code and other content. It also allows us to revert changes to a previous point in time.

Types of Version control

There are many forms of version control. Some not as good:

• Save a document with a new date (we’ve all done it, but it isn’t efficient)
• Google Docs "history" function (not bad for some documents, but limited in scope).

Some better:

• Mercurial
• Subversion
• Git - which we’ll be learning much more about in this series.

More Resources:

Visit the version control Wikipedia list of version control platforms.

• Read the Git documentation explaining the progression of version control systems.

Why Version Control is Important

Version control facilitates two important aspects of many scientific workflows:

1. The ability to save and review or revert to previous versions.
2. The ability to collaborate on a single project.

This means that you don’t have to worry about a collaborator (or your future self) overwriting something important. It also allows two people working on the same document to efficiently combine ideas and changes.

How Version Control Systems Works

Simple Version Control Model

A version control system keeps track of what has changed in one or more files over time. The way this tracking occurs, is slightly different between various version control tools including git, mercurial and svn. However the principle is the same.

Version control systems begin with a base version of a document. They then save the committed changes that you make. You can think of version control as a tape: if you rewind the tape and start at the base document, then you can play back each change and end up with your latest version.

Once you think of changes as separate from the document itself, you can then think about “playing back” different sets of changes onto the base document. You can then retrieve, or revert to, different versions of the document.

The benefit of version control when you are in a collaborative environment is that two users can make independent changes to the same document.

If there aren’t conflicts between the users changes (a conflict is an area where both users modified the same part of the same document in different ways) you can review two sets of changes on the same base document.

A version control system is a tool that keeps track of these changes for us. Each version of a file can be viewed and reverted to at any time. That way if you add something that you end up not liking or delete something that you need, you can simply go back to a previous version.

Git & GitHub - A Distributed Version Control Model

GitHub uses a distributed version control model. This means that there can be many copies (or forks in GitHub world) of the repository.

Have a look at the graphic below. Notice that in the example, there is a "central" version of our repository. Joe, Sue and Eve are all working together to update the central repository. Because they are using a distributed system, each user (Joe, Sue and Eve) has their own copy of the repository and can contribute to the central copy of the repository at any time.

Create A Working Copy of a Git Repo - Fork

There are many different Git and GitHub workflows. In the NEON Data Institute, we will use a distributed workflow with a Central Repository. This allows us all (all of the Institute participants) to work independently. We can then contribute our changes to update the Central (NEON) Repository. Our collaborative workflow goes like this:

1. NEON "owns" the Central Repository.
2. You will create a copy of this repository (known as a fork) in your own GitHub account.
3. You will then clone (copy) the repository to your local computer. You will do your work locally on your laptop.
4. When you are ready to submit your changes to the NEON repository, you will:
   • Sync your local copy of the repository with NEON's central repository so you have the most up to date version, and then,
   • Push the changes you made to your local copy (or fork) of the repository to NEON's main repository.

Each participant in the institute will be contributing to the NEON central repository using the same workflow! Pretty cool stuff.

Let's get some terms straight before we go any further.

• Central repository - the central repository is what all participants will add to. It is the "final working version" of the project.
• Your forked repository - is a "personal” working copy of the central repository stored in your GitHub account. This is called a fork. When you are happy with your work, you update your repo from the central repo, then you can update your changes to the central NEON repository.
• Your local repository - this is a local version of your fork on your own computer. You will most often do all of your work locally on your computer.

Additional Resources:

Further documentation for and how-to-use direction for Git, is provided by the Git Pro version 2 book by Scott Chacon and Ben Straub , available in print or online. If you enjoy learning from videos, the site hosts several.

This tutorial builds upon the previous tutorial, to work with shapefile attributes in R and explores how to plot multiple shapefiles using base R graphics. It then covers how to create a custom legend with colors and symbols that match your plot.

Load the Data

To work with vector data in R, we can use the rgdal library. The raster package also allows us to explore metadata using similar commands for both raster and vector files.

We will import three shapefiles. The first is our AOI or area of interest boundary polygon that we worked with in Open and Plot Shapefiles in R. The second is a shapefile containing the location of roads and trails within the field site. The third is a file containing the Harvard Forest Fisher tower location. These latter two we worked with in the Explore Shapefile Attributes & Plot Shapefile Objects by Attribute Value in R tutorial.

# load packages
# rgdal: for vector work; sp package should always load with rgdal. 
library(rgdal)  
# raster: for metadata/attributes- vectors or rasters
library(raster)   

# set working directory to data folder
# setwd("pathToDirHere")

# Import a polygon shapefile 
aoiBoundary_HARV <- readOGR("NEON-DS-Site-Layout-Files/HARV",
                            "HarClip_UTMZ18", stringsAsFactors = T)

## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/olearyd/Git/data/NEON-DS-Site-Layout-Files/HARV", layer: "HarClip_UTMZ18"
## with 1 features
## It has 1 fields
## Integer64 fields read as strings:  id

# Import a line shapefile
lines_HARV <- readOGR( "NEON-DS-Site-Layout-Files/HARV", "HARV_roads", stringsAsFactors = T)

## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/olearyd/Git/data/NEON-DS-Site-Layout-Files/HARV", layer: "HARV_roads"
## with 13 features
## It has 15 fields

# Import a point shapefile 
point_HARV <- readOGR("NEON-DS-Site-Layout-Files/HARV",
                      "HARVtower_UTM18N", stringsAsFactors = T)

## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/olearyd/Git/data/NEON-DS-Site-Layout-Files/HARV", layer: "HARVtower_UTM18N"
## with 1 features
## It has 14 fields

Plot Data

In the Explore Shapefile Attributes & Plot Shapefile Objects by Attribute Value in R tutorial we created a plot where we customized the width of each line in a spatial object according to a factor level or category. To do this, we create a vector of colors containing a color value for EACH feature in our spatial object grouped by factor level or category.

# view the factor levels
levels(lines_HARV$TYPE)

## [1] "boardwalk"  "footpath"   "stone wall" "woods road"

# create vector of line width values
lineWidth <- c(2,4,3,8)[lines_HARV$TYPE]
# view vector
lineWidth

##  [1] 8 4 4 3 3 3 3 3 3 2 8 8 8

# create a color palette of 4 colors - one for each factor level
roadPalette <- c("blue","green","grey","purple")
roadPalette

## [1] "blue"   "green"  "grey"   "purple"

# create a vector of colors - one for each feature in our vector object
# according to its attribute value
roadColors <- c("blue","green","grey","purple")[lines_HARV$TYPE]
roadColors

##  [1] "purple" "green"  "green"  "grey"   "grey"   "grey"   "grey"  
##  [8] "grey"   "grey"   "blue"   "purple" "purple" "purple"

# create vector of line width values
lineWidth <- c(2,4,3,8)[lines_HARV$TYPE]
# view vector
lineWidth

##  [1] 8 4 4 3 3 3 3 3 3 2 8 8 8

# in this case, boardwalk (the first level) is the widest.
plot(lines_HARV, 
     col=roadColors,
     main="NEON Harvard Forest Field Site\n Roads & Trails \nLine Width Varies by Type Attribute Value",
     lwd=lineWidth)

Add Plot Legend

In the the previous tutorial, we also learned how to add a basic legend to our plot.

• bottomright: We specify the location of our legend by using a default keyword. We could also use top, topright, etc.
• levels(objectName$attributeName): Label the legend elements using the categories of levels in an attribute (e.g., levels(lines_HARV$TYPE) means use the levels boardwalk, footpath, etc).
• fill=: apply unique colors to the boxes in our legend. palette() is the default set of colors that R applies to all plots.

Let's add a legend to our plot.

plot(lines_HARV, 
     col=roadColors,
     main="NEON Harvard Forest Field Site\n Roads & Trails\n Default Legend")

# we can use the color object that we created above to color the legend objects
roadPalette

## [1] "blue"   "green"  "grey"   "purple"

# add a legend to our map
legend("bottomright", 
       legend=levels(lines_HARV$TYPE), 
       fill=roadPalette, 
       bty="n", # turn off the legend border
       cex=.8) # decrease the font / legend size

However, what if we want to create a more complex plot with many shapefiles and unique symbols that need to be represented clearly in a legend?

Plot Multiple Vector Layers

Now, let's create a plot that combines our tower location (point_HARV), site boundary (aoiBoundary_HARV) and roads (lines_HARV) spatial objects. We will need to build a custom legend as well.

To begin, create a plot with the site boundary as the first layer. Then layer the tower location and road data on top using add=TRUE.

# Plot multiple shapefiles
plot(aoiBoundary_HARV, 
     col = "grey93", 
     border="grey",
     main="NEON Harvard Forest Field Site")

plot(lines_HARV, 
     col=roadColors,
     add = TRUE)

plot(point_HARV, 
     add  = TRUE, 
     pch = 19, 
     col = "purple")

# assign plot to an object for easy modification!
plot_HARV<- recordPlot()

Customize Your Legend

Next, let's build a custom legend using the symbology (the colors and symbols) that we used to create the plot above. To do this, we will need to build three things:

1. A list of all "labels" (the text used to describe each element in the legend to use in the legend.
2. A list of colors used to color each feature in our plot.
3. A list of symbols to use in the plot. NOTE: we have a combination of points, lines and polygons in our plot. So we will need to customize our symbols!

Let's create objects for the labels, colors and symbols so we can easily reuse them. We will start with the labels.

# create a list of all labels
labels <- c("Tower", "AOI", levels(lines_HARV$TYPE))
labels

## [1] "Tower"      "AOI"        "boardwalk"  "footpath"   "stone wall"
## [6] "woods road"

# render plot
plot_HARV

# add a legend to our map
legend("bottomright", 
       legend=labels, 
       bty="n", # turn off the legend border
       cex=.8) # decrease the font / legend size

Now we have a legend with the labels identified. Let's add colors to each legend element next. We can use the vectors of colors that we created earlier to do this.

# we have a list of colors that we used above - we can use it in the legend
roadPalette

## [1] "blue"   "green"  "grey"   "purple"

# create a list of colors to use 
plotColors <- c("purple", "grey", roadPalette)
plotColors

## [1] "purple" "grey"   "blue"   "green"  "grey"   "purple"

# render plot
plot_HARV

# add a legend to our map
legend("bottomright", 
       legend=labels, 
       fill=plotColors,
       bty="n", # turn off the legend border
       cex=.8) # decrease the font / legend size

Great, now we have a legend! However, this legend uses boxes to symbolize each element in the plot. It might be better if the lines were symbolized as a line and the points were symbolized as a point symbol. We can customize this using pch= in our legend: 16 is a point symbol, 15 is a box.

# create a list of pch values
# these are the symbols that will be used for each legend value
# ?pch will provide more information on values
plotSym <- c(16,15,15,15,15,15)
plotSym

## [1] 16 15 15 15 15 15

# Plot multiple shapefiles
plot_HARV

# to create a custom legend, we need to fake it
legend("bottomright", 
       legend=labels,
       pch=plotSym, 
       bty="n", 
       col=plotColors,
       cex=.8)

Now we've added a point symbol to represent our point element in the plot. However it might be more useful to use line symbols in our legend rather than squares to represent the line data. We can create line symbols, using lty = (). We have a total of 6 elements in our legend:

1. A Tower Location
2. An Area of Interest (AOI)
3. and 4 Road types (levels)

The lty list designates, in order, which of those elements should be designated as a line (1) and which should be designated as a symbol (NA). Our object will thus look like lty = c(NA,NA,1,1,1,1). This tells R to only use a line element for the 3-6 elements in our legend.

Once we do this, we still need to modify our pch element. Each line element (3-6) should be represented by a NA value - this tells R to not use a symbol, but to instead use a line.

# create line object
lineLegend = c(NA,NA,1,1,1,1)
lineLegend

## [1] NA NA  1  1  1  1

plotSym <- c(16,15,NA,NA,NA,NA)
plotSym

## [1] 16 15 NA NA NA NA

# plot multiple shapefiles
plot_HARV

# build a custom legend
legend("bottomright", 
       legend=labels, 
       lty = lineLegend,
       pch=plotSym, 
       bty="n", 
       col=plotColors,
       cex=.8)

In this tutorial, we will cover the R knitr package that is used to convert R Markdown into a rendered document (HTML, PDF, etc).

Share & Publish Results Directly from Your Code!

The knitr package allow us to:

• Publish & share preliminary results with collaborators.
• Create professional reports that document our workflow and results directly from our code, reducing the risk of accidental copy and paste or transcription errors.
• Document our workflow to facilitate reproducibility.
• Efficiently change code outputs (figures, files) given changes in the data, methods, etc.

Publish from Rmd files with knitr

To complete this tutorial you need:

1. The R knitr package to complete this tutorial. If you need help installing packages, visit the R packages tutorial.
2. An R Markdown document that contains a YAML header, code chunks and markdown segments. If you don't have an .Rmd file, visit the R Markdown tutorial to create one.

How to Knit

To knit in RStudio, click the knit pull down button. You want to use the
knit HTML for this lesson.

When you click the Knit HTML button, a window will open in your console titled R Markdown. This pane shows the knitting progress. The output (HTML in this case) file will automatically be saved in the current working directory. If there is an error in the code, an error message will appear with a line number in the R Console to help you diagnose the problem.

View the Output

When knitting is complete, the new HTML file produced will automatically open.

Notice that information from the YAML header (title, author, date) is printed at the top of the HTML document. Then the HTML shows the text, code, and results of the code that you included in the RMD document.

You will need to have the rmarkdown and knitr packages installed on your computer prior to completing this tutorial. Refer to the setup materials to get these installed.

Learning Objectives

At the end of this activity, you will:

• Know how to create an R Markdown file in RStudio.
• Be able to write a script with text and R code chunks.
• Create an R Markdown document ready to be ‘knit’ into an HTML document to share your code and results.

Things You’ll Need To Complete This Tutorial

You will need the most current version of R and, preferably, RStudio loaded on your computer to complete this tutorial.

Install R Packages

• knitr: install.packages("knitr")
• rmarkdown: install.packages("rmarkdown")
• raster: install.packages("raster")
• rgdal: install.packages("rgdal")

Download Data

Download NEON Teaching Data Subset: TEAK-Data Institute 2016

The LiDAR and imagery data used to create this raster teaching data subset were collected over the National Ecological Observatory Network's (NEON) Lower Teakettle field site and processed at NEON headquarters. The entire dataset can be accessed by request from the NEON Data Portal.

Download Dataset

You will want to create a data directory for all the Data Institute teaching datasets. We suggest the pathway be ~/Documents/data/NEONDI-2016 or the equivalent for your operating system. Once you've downloaded and unzipped the dataset, move it to this directory.

Additional Resources

• R Markdown Cheatsheet: a very handy reference for using R Markdown
• R Markdown Reference Guide: a more expensive reference for R Markdown
• Introduction to R Markdown by Garrett Grolemund: a tutorial for learning R Markdown

Create an Rmd File

RMarkdown in RStudio Video

Our goal in this series is to document our workflow. We can do this by

1. Creating an R Markdown (RMD) file in R studio and
2. Rendering that RMD file to HTML using knitr.

Watch this 6:38 minute video below to learn more about how you can convert an R Markdown file to HTML (or other formats) using knitr in RStudio. The text size in the video is small so you may want to watch the video in full screen mode.

Now that you have a sense of how R Markdown can be used in RStudio, you are ready to create your own RMD document. Do the following:

1. Create a new R Markdown file and choose HTML as the desired output format.
2. Enter a Title (Explore NEON LiDAR Data) and Author Name (your name). Then click OK.
3. Save the file using the following format: LastName-institute-week3.rmd NOTE: The document title is not the same as the file name.
4. Hit the knit button in RStudio (as is done in the video above). What happens?

If everything went well, you should have an HTML format (web page) output after you hit the knit button. Note that this HTML output is built from a combination of code and documentation that was written using markdown syntax.

Next, we'll break down the structure of an R Markdown file.

Understand Structure of an R Markdown file

Let's next review the structure of an R Markdown (.Rmd) file. There are three main content types:

• Header: the text at the top of the document, written in YAML format.
• Markdown sections: text that describes your workflow written using markdown syntax.
• Code chunks: Chunks of R code that can be run and also can be rendered using knitr to an output document.

Next let's explore each section type.

Header -- YAML block

An R Markdown file always starts with header written using YAML syntax. There are four default elements in the RStudio generated YAML header:

• title: the title of your document. Note, this is not the same as the file name.
• author: who wrote the document.
• date: by default this is the date that the file is created.
• output: what format will the output be in. We will use HTML.

A YAML header may be structured differently depending upon how your are using it. Learn more on the R Markdown documentation page.

R Markdown Text/Markdown Blocks

An RMD document contains a mixture of code chunks and markdown blocks where you can describe aspects of your processing workflow. The markdown blocks use the same markdown syntax that we learned last week in week 2 materials. In these blocks you might describe the data that you are using, how it's being processed and and what the outputs are. You may even add some information that interprets the outputs.

When you render your document to HTML, this markdown will appear as text on the output HTML document.

Learn More about RStudio Markdown Basics

Explore Your R Markdown File

Look closely at the pre-populated markdown and R code chunks in your RMD file.

Does any of the markdown syntax look familiar?

• Are any words in bold?
• Are any words in italics?
• Are any words highlighted as code?

If you are unsure, the answers are at the bottom of this page.

Code chunks

Code chunks are where your R code goes. All code chunks start and end with ``` – three backticks or graves. On your keyboard, the backticks can be found on the same key as the tilde. Graves are not the same as an apostrophe!

The initial line of a code chunk must appear as:

 ```{r chunk-name-with-no-spaces} 
# code goes here
 ```

The r part of the chunk header identifies this chunk as an R code chunk and is mandatory. Next to the {r, there is a chunk name. This name is not required for basic knitting however, it is good practice to give each chunk a unique name as it is required for more advanced knitting approaches.

```{r setup-library }
   
   library(rgdal)
   library(raster)
 
 ```

```{r code-setwd}

   # set working directory to ensure R can find the file we wish to import.
   # This will depend on your local environment.
   setwd("~/Documents/data/NEONDI-2016/") 

```

You changed the working directory to ~/Documents/data/NEONDI-2016/ (probably via setwd()). It will be restored to [directory path of current .rmd file]. See the Note section in ?knitr::knit ?knitr::knit

```{r code-setwd-stringvariable}

	# set working directory as a string variable for use in other code chunks.
	# This will depend on your local environment.
	wd <- "~/Documents/data/NEONDI-2016/"
	setwd(wd)

 ```

```{r load-dsm-raster }
   # check for the working directory
   getwd()
   # In this new chunk, the working directory has reverted to default upon kitting. 
   # Combining the working directory string variable and 
   # additional path to the file, import a DSM file.
   teak_dsm <- raster(paste0(wd, "NEONdata/D17-California/TEAK/2013/lidar/TEAK_lidarDSM.tif"))
   
 ```

Code chunk options

You can also add arguments or options to each code chunk. These arguments allow you to customize how or if you want code to be processed or appear on the output HTML document. Code chunk arguments are added on the first line of a code chunk after the name, within the curly brackets.

The example below, is a code chunk that will not be "run", or evaluated, by R. The code within the chunk will appear on the output document, however there will be no outputs from the code.

```{r intro-option, eval=FALSE}
# the code here will not be processed by R 
# but it will appear on your output document
1+2
 ```

We use eval=FALSE often when the chunk is exporting an file that we don't need to re-export but we want to document the code used to export the file.

Three common code chunk options are:

• eval = FALSE: Do not evaluate (or run) this code chunk when knitting the RMD document. The code in this chunk will still render in our knitted HTML output, however it will not be evaluated or run by R.
• echo = FALSE: Hide the code in the output. The code is evaluated when the RMD file is knit, however only the output is rendered on the output document.
• results = hide: The code chunk will be evaluated but the results of the code will not be rendered on the output document. This is useful if you are viewing the structure of a large object (e.g. outputs of a large data.frame).

Multiple code chunk options can be used for the same chunk. For more on code chunk options, read R Markdown: The Definitive Guide or the knitr documentation.

Now continue on to the next tutorial to learn how to knit this document into a HTML file.

This tutorial we will work with the knitr and rmarkdown packages within RStudio to learn how to effectively and efficiently document and publish our workflows online.

Documentation Is Important

As we read in the Reproducible Science overview, the four facets of reproducible science are:

• Documentation
• Organization,
• Automation and
• Dissemination.

This week we will learn about the R Markdown file format (and R package) which can be used with the knitr package to document and publish (disseminate) your code and code output.

View Slideshow: Share, Publish & Archive - from the Reproducible Science Curriculum

The Tools We Will Use

R Markdown

“R Markdown is an authoring format that enables easy creation of dynamic documents, presentations, and reports from R. It combines the core syntax of markdown (an easy to write plain text format) with embedded R code chunks that are run so their output can be included in the final document. R Markdown documents are fully reproducible (they can be automatically regenerated whenever underlying R code or data changes)." -- RStudio documentation.

We use markdown syntax in R Markdown (.rmd) files to document workflows and to share data processing, analysis and visualization outputs. We can also use it to create documents that combine R code, output and text.

Why R Markdown?

There are many advantages to using R Markdown in your work:

• Human readable syntax.
• Simple syntax - it can be learned quickly.
• All components of your work are clearly documented. You don't have to remember what steps, assumptions, tests were used.
• You can easily extend or refine analyses by modifying existing or adding new code blocks.
• Analysis results can be disseminated in various formats including HTML, PDF, slide shows and more.
• Code and data can be shared with a colleague to replicate the workflow.

Knitr

The knitr package for R allows us to create readable documents from R Markdown files.

In the next tutorial we will learn more about working with the R Markdown format in RStudio.

The primary goal of this tutorial is to explain how to set a working directory in R. The working directory is where your R session interacts with your hard drive. This is where you can read data that you want to use, and save new information such as derived data products, tables, maps, and figures. It is a good practice to store your information in an organized set of directories, so you will often want to change your working directory depending on what kinds of information that you need to access. This tutorial teaches how to download and unzip the data files that accompany many NEON Data Skills tutorials, and also covers the concept of file paths. You can read from top to bottom, or use the menu bar at left to navigate to your desired topic.

Set Up For NEON Data Skills Tutorials

Many NEON data tutorials utilize teaching data subsets which are hosted on the NEON Data Skills figshare repository. If a data subset is required for a tutorial it can be downloaded at the top of each tutorial in the Download Data section.

Prior to working with any data in R, we must set the working directory to the location of the data files. Setting the working directory tells R where the data files are located on the computer. If the working directory is not correctly set first, when we try to open a file we will get an error telling us that R cannot find the file.

• Wondering why we're saying directory instead of folder? See our discussion of Directory vs. Folder in the middle of this tutorial.

Download & Uncompress the Data

1) Download

First, we will download the data to a location on the computer. To download the data for this tutorial, click the blue button Download NEON Teaching Data Subset: Meteorological Data for Harvard Forest within the box at the top of this page.

Note: In other NEON Data Skills tutorials, download all data subsets in the Download Data section prior to starting the tutorial. Here, the second data subset is for those wishing to practice these skills in a Challenge activity and will be downloaded at that time.

After clicking on the Download Data button, the data will automatically download to the computer.

2) Locate .zip file

Second, we need to find the downloaded .zip file. Many browsers default to downloading to the Downloads directory on your computer. Note: You may have previously specified a specific directory (folder) for files downloaded from the internet, if so, the .zip file will download there.

3) Move to data directory

Third, we must move the data files to the location we want to work with them. We recommend moving the .zip to a dedicated data directory within the Documents directory on your computer. This data directory can then be a repository for all data subsets you use for the NEON Data Skills tutorials. Note: If you chose to store your data in a different directory (e.g., not in ~/Documents/data), modify the directions below with the appropriate file path to your data directory.

4) Unzip/uncompress

Fourth, we need to unzip/uncompress the file so that the data files can be accessed. Use your favorite tool that can unpackage/open .zip files (e.g., winzip, Archive Utility, etc). The files will now be accessible in a directory named NEON-DS-Met-Time-Series containing all the subdirectories and files that make up the dataset or the subdirectories and files will be unzipped directly into the data directory. If the latter happens, they need to be moved into a data/NEON-DS-Met-Time-Series directory.

The directory should be the same as in this screenshot (below). Note that NEON-DS-Site-Lyout-Files directory will only be in your directory if you completed the challenge above. If you did not, your directory should look the same but without that directory.

Directory vs. Folder

"Directory" and "Folder" both refer to the same thing. Folder makes a lot of sense when we think of an isolated folder as a "bin" containing many files. However, the analogy to a physical file folder falters when we start thinking about the relationship between different folders and how we tell a computer to find a specific folder. This is why the term directory is often preferred. Any directory (folder) can hold other directories and/or files. When we set the working directory, we are telling the computer the location of the directory (or folder) to start with when looking for other files or directories, or to save any output to.

Full, Base, and Relative Paths

The data downloaded and unzipped in the previous steps are located within a nested set of directories:

• primary-level/home directory: neon
  • This directory isn't obvious as we are within this directory once we log into the computer.
  • You will see your own user ID.
• secondary-level directory: neon/Documents
• tertiary-level directory: neon/Documents/data
• quaternary-level directory: neon/Documents/data/NEON-DS-Met-Time-Series
• quaternary-level directory: neon/Documents/data/NEON-DS-Site-Layout-Shapefiles

Full & Base Path

The full path is essentially the complete "directions" for how to find the desired directory or file. It always starts with the home directory or root (e.g., /Users/neon/). A full path is the base path when used to set the working directory to a specific directory. The base path for the NEON-DS-Met-Time-Series directory would be:

 /Users/neon/Documents/data/NEON-DS-Met-Time-Series 

Relative Path

A relative path is a path to a directory or file that starts from the location determined by the working directory. If our working directory is set to the data directory,

 /Users/neon/Documents/data/

we can then create a relative path for all directories and files within the data directory.

The relative path for the meanNDVI_HARV_2011.csv file would be:

 NEON-DS-Met-Time-Series/HARV/NDVI/meanNDVI_HARV_2011.csv

The R Working Directory

In R, the working directory is the directory where R starts when looking for any file to open (as directed by a file path) and where it saves any output.

Without a working directory all R scripts would need the full file path written any time we wanted to open or save a file. It is more efficient if we have a base file path set as our working directory and then all file paths written in our scripts only consist of the file path relative to that base path (a relative path).

• If you are unfamiliar with the term full path, base path, or relative path, please see the section below on Full, Base, and Relative Paths for a more detailed explanation before continuing with this tutorial.

Find a Full Path to a File in Unknown Location

If you are unsure of the path to a specific directory or file, you can find this information for a particular file/directory of interest by looking in the:

• Windows: Properties, General tab (right click on the file/directory) or in the file path bar at the top of each window (select versions).
• Mac OS X: Get Info (right clicking/control+click on the file/directory)

The file path may appear as:

Computer > Users > neon > Documents > data > NEON-DS-Met-Time-Series

Copy and paste this information to automatically reformat into the full path to the directory or file:

• Windows: C:\Users\neon\Documents\data\NEON-DS-Met-Time-Series
• Mac OS X: /Users/neon/Documents/data/NEON-DS-Met-Time-Series

Determine Current Working Directory

Once we are in the R program, we can view the current working directory using the code getwd().

# view current working directory 
getwd()
[1] "/Users/neon"

The working directory is currently set to the home directory /Users/neon. Remember, your current working directory will have a different user name and may appear different based on operating system.

This code can be used at any time to determine the current working directory.

Set the Working Directory

To set our current working directory to the location where our data are located, we can either set the working directory in the R script or use our current GUI to select the working directory.

We want to set our working directory to the data directory.

Set the Working Directory: Base Path in Script

We can set the working directory using the code setwd("PATH") where PATH is the full path to the desired directory. You can enter "PATH" as a string (as shown below), or save the file path as a string to a variable (e.g., wd <- "~/Documents/data") and then set the working directory based on that variable (e.g., setwd(wd)).

This latter method is used in many of the NEON Data Skills tutorials because of the advantages that this method provides. First, this method is extermely flexible across different compute environments and formats, including personal computers, Linux-based servers on 'the cloud' (e.g., AWS, CyVerse), and when using Rmarkdown (.Rmd) documents. Second this method allows the tutorial's user to set their working directory once as a string and then to reuse that string as needed to reference input files, or make output files. For example, some functions must have a full filepath to an input file (such as when reading in HDF5 files). Third, this method simplifies the process that NEON uses internally to create and update these tutorials. All in all, saving the working directory as a string variable makes the code more explicit and determanistic without relying on working knowledge of relative filepaths, making your code more producible and easier for an outsider to interpret.

To practice, use these techniques to set your working directory to the directory where you have the data saved, and check that you set the working directory correctly. There is no R output from setwd(). If we want to check that the working directory is correctly set we can use getwd().

Example Windows File Path

Notice the the backslashes used in Windows paths must be changed to slashes in R.

# set the working directory to `data` folder
wd <- "C:/Users/neon/Documents/data"
setwd(wd)

# check to ensure path is correct
getwd()
[1] "C:/Users/neon/Documents/data"

Example Mac OS X File Path

# set the working directory to `data` folder
wd <- "/Users/neon/Documents/data"
setwd(wd)

# check to ensure path is correct
getwd()
[1] "/Users/neon/Documents/data"

Set the Working Directory: Using RStudio GUI

You can also set the working directory using the Rstudio and/or R graphical user interface (GUI). This method is easy for beginners to learn, but it also makes your code less reproducible because it relies on a person to follow certain instructions, which is a process that introduces human error. It may also be impossible for an observer to determine where your input data are stored, which can make troubleshooting more difficult as well. Use this method when getting started, or when you will find it helpful to use a graphical user interface to navigate your files. Note that this method will run a single line setwd() command in the console when you select your working directory, so you can copy/paste that line into your script for future use!

1. Go to Session in menu bar,
2. select Select Working Directory,
3. select Choose Directory,
4. in the new window that appears, select the appropriate directory.

Set the Working Directory: Using R GUI

Windows Operating Systems:

1. Go to the File menu bar,
2. select Change dir... or Change Working Directory,
3. in the new window that appears, select the appropriate directory.

Mac Operating Systems:

1. Go to the Misc menu,
2. select Change Working Directory,
3. in the new window that appears, select the appropriate directory.

This tutorial explores how to import and plot a multiband raster in R. It also covers how to plot a three-band color image using the plotRGB() function in R.

The Basics of Imagery - About Spectral Remote Sensing Data

About Raster Bands in R

As discussed in the Intro to Raster Data tutorial, a raster can contain 1 or more bands.

To work with multiband rasters in R, we need to change how we import and plot our data in several ways.

• To import multiband raster data we will use the stack() function.
• If our multiband data are imagery that we wish to composite, we can use plotRGB() (instead of plot()) to plot a 3 band raster image.

About MultiBand Imagery

One type of multiband raster dataset that is familiar to many of us is a color image. A basic color image consists of three bands: red, green, and blue. Each band represents light reflected from the red, green or blue portions of the electromagnetic spectrum. The pixel brightness for each band, when composited creates the colors that we see in an image.

Getting Started with Multi-Band Data in R

To work with multiband raster data we will use the terra package.

# terra package to work with raster data

library(terra)



# package for downloading NEON data

library(neonUtilities)



# package for specifying color palettes

library(RColorBrewer)



# set working directory to ensure R can find the file we wish to import

wd <- "~/data/" # this will depend on your local environment environment

# be sure that the downloaded file is in this directory

setwd(wd)

In this tutorial, the multi-band data that we are working with is imagery collected using the NEON Airborne Observation Platform high resolution camera over the NEON Harvard Forest field site. Each RGB image is a 3-band raster. The same steps would apply to working with a multi-spectral image with 4 or more bands - like Landsat imagery, or even hyperspectral imagery (in geotiff format). We can plot each band of a multi-band image individually.

byTileAOP(dpID='DP3.30010.001', # rgb camera data
          site='HARV',
          year='2022',
          easting=732000,
          northing=4713500,
          check.size=FALSE, # set to TRUE or remove if you want to check the size before downloading
          savepath = wd)

## Downloading files totaling approximately 351.004249 MB 
## Downloading 1 files
## 

|
| | 0% |
|====================================================================================================================================================================| 100% ## Successfully downloaded 1 files to ~/data//DP3.30010.001

Or we can plot each bands separately as follows:

# Determine the number of bands

num_bands <- nlyr(RGB_HARV)



# Define colors to plot each

# Define color palettes for each band using RColorBrewer

colors <- list(
  brewer.pal(9, "Reds"),
  brewer.pal(9, "Greens"),
  brewer.pal(9, "Blues")
)



# Plot each band in a loop, with the specified colors

for (i in 1:num_bands) {
  plot(RGB_HARV[[i]], main=paste("Band", i), col=colors[[i]])
}

Image Raster Data Attributes

We can display some of the attributes about the raster, as shown below:

# Print dimensions

cat("Dimensions:\n")

## Dimensions:

cat("Number of rows:", nrow(RGB_HARV), "\n")

## Number of rows: 10000

cat("Number of columns:", ncol(RGB_HARV), "\n")

## Number of columns: 10000

cat("Number of layers:", nlyr(RGB_HARV), "\n")

## Number of layers: 3

# Print resolution

resolutions <- res(RGB_HARV)

cat("Resolution:\n")

## Resolution:

cat("X resolution:", resolutions[1], "\n")

## X resolution: 0.1

cat("Y resolution:", resolutions[2], "\n")

## Y resolution: 0.1

# Get the extent of the raster

rgb_extent <- ext(RGB_HARV)



# Convert the extent to a string with rounded values

extent_str <- sprintf("xmin: %d, xmax: %d, ymin: %d, ymax: %d", 
                      round(xmin(rgb_extent)), 
                      round(xmax(rgb_extent)), 
                      round(ymin(rgb_extent)), 
                      round(ymax(rgb_extent)))



# Print the extent string

cat("Extent of the raster: \n")

## Extent of the raster:

cat(extent_str, "\n")

## xmin: 732000, xmax: 733000, ymin: 4713000, ymax: 4714000

Let's take a look at the coordinate reference system, or CRS. You can use the parameters describe=TRUE to display this information more succinctly.

crs(RGB_HARV, describe=TRUE)

##                    name authority  code
## 1 WGS 84 / UTM zone 18N      EPSG 32618
##                                                                                                                                                                                                                                                          area
## 1 Between 78°W and 72°W, northern hemisphere between equator and 84°N, onshore and offshore. Bahamas. Canada - Nunavut; Ontario; Quebec. Colombia. Cuba. Ecuador. Greenland. Haiti. Jamaica. Panama. Turks and Caicos Islands. United States (USA). Venezuela
##            extent
## 1 -78, -72, 84, 0

Let's next examine the raster's minimum and maximum values. What is the range of values for each band?

# Replace Inf and -Inf with NA

values(RGB_HARV)[is.infinite(values(RGB_HARV))] <- NA



# Get min and max values for all bands

min_max_values <- minmax(RGB_HARV)



# Print the results

cat("Min and Max Values for All Bands:\n")

## Min and Max Values for All Bands:

print(min_max_values)

##     2022_HARV_7_732000_4713000_image_1 2022_HARV_7_732000_4713000_image_2 2022_HARV_7_732000_4713000_image_3
## min                                  0                                  0                                  0
## max                                255                                255                                255

This raster contains values between 0 and 255. These values represent the intensity of brightness associated with the image band. In the case of a RGB image (red, green and blue), band 1 is the red band. When we plot the red band, larger numbers (towards 255) represent pixels with more red in them (a strong red reflection). Smaller numbers (towards 0) represent pixels with less red in them (less red was reflected). To plot an RGB image, we mix red + green + blue values into one single color to create a full color image - this is the standard color image a digital camera creates.

Other Types of Multi-band Raster Data

Multi-band raster data might also contain:

1. Time series: the same variable, over the same area, over time.
2. Multi or hyperspectral imagery: image rasters that have 4 or more (multi-spectral) or more than 10-15 (hyperspectral) bands. Check out the NEON Data Skills Imaging Spectroscopy HDF5 in R tutorial to learn how to work with hyperspectral data cubes.

The true color image plotted at the beginning of this lesson looks pretty decent. We can explore whether applying a stretch to the image might improve clarity and contrast using stretch="lin" or stretch="hist".

# What does stretch do?



# Plot the linearly stretched raster

plotRGB(RGB_HARV, stretch="lin")

# Plot the histogram-stretched raster

plotRGB(RGB_HARV, stretch="hist")

In this case, the stretch doesn't enhance the contrast our image significantly given the distribution of reflectance (or brightness) values is distributed well between 0 and 255, and applying a stretch appears to introduce some artificial, almost purple-looking brightness to the image.

In this tutorial, we will review the fundamental principles, packages and metadata/raster attributes that are needed to work with raster data in R. We discuss the three core metadata elements that we need to understand to work with rasters in R: CRS, extent and resolution. We also explore missing and bad data values as stored in a raster and how R handles these elements. Finally, we introduce the GeoTiff file format.

About Raster Data

Raster or "gridded" data are stored as a grid of values which are rendered on a map as pixels. Each pixel value represents an area on the Earth's surface.

Raster Data in R

Let's first import a raster dataset into R and explore its metadata. To open rasters in R, we will use the terra package.

library(terra)



# set working directory, you can change this if desired

wd <- "~/data/" 

setwd(wd)

Download LiDAR Raster Data

We can use the neonUtilities function byTileAOP to download a single elevation tiles (DSM and DTM). You can run help(byTileAOP) to see more details on what the various inputs are. For this exercise, we'll specify the UTM Easting and Northing to be (732000, 4713500), which will download the tile with the lower left corner (732000,4713000). By default, the function will check the size total size of the download and ask you whether you wish to proceed (y/n). This file is ~8 MB, so make sure you have enough space on your local drive. You can set check.size=TRUE if you want to check the file size before downloading.

byTileAOP(dpID='DP3.30024.001', # lidar elevation

          site='HARV',

          year='2022',

          easting=732000,

          northing=4713500,

          check.size=FALSE, # set to TRUE or remove if you want to check the size before downloading

          savepath = wd)

This file will be downloaded into a nested subdirectory under the ~/data folder, inside a folder named DP3.30024.001 (the Data Product ID). The file should show up in this location: ~/data/DP3.30024.001/neon-aop-products/2022/FullSite/D01/2022_HARV_7/L3/DiscreteLidar/DSMGtif/NEON_D01_HARV_DP3_732000_4713000_DSM.tif.

Open a Raster in R

We can use terra's rast("path-to-raster-here") function to open a raster in R.

# Load raster into R

dsm_harv_file <- paste0(wd, "DP3.30024.001/neon-aop-products/2022/FullSite/D01/2022_HARV_7/L3/DiscreteLidar/DSMGtif/NEON_D01_HARV_DP3_732000_4713000_DSM.tif")

DSM_HARV <- rast(dsm_harv_file)



# View raster structure

DSM_HARV 

## class       : SpatRaster 
## dimensions  : 1000, 1000, 1  (nrow, ncol, nlyr)
## resolution  : 1, 1  (x, y)
## extent      : 732000, 733000, 4713000, 4714000  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 18N (EPSG:32618) 
## source      : NEON_D01_HARV_DP3_732000_4713000_DSM.tif 
## name        : NEON_D01_HARV_DP3_732000_4713000_DSM 
## min value   :                               317.91 
## max value   :                               433.94

# plot raster 

plot(DSM_HARV, main="Digital Surface Model - HARV")

Types of Data Stored in Raster Format

Raster data can be continuous or categorical. Continuous rasters can have a range of quantitative values. Some examples of continuous rasters include:

1. Precipitation maps.
2. Maps of tree height derived from LiDAR data.
3. Elevation values for a region.

The raster we loaded and plotted earlier was a digital surface model, or a map of the elevation for Harvard Forest derived from the NEON AOP LiDAR sensor. Elevation is represented as a continuous numeric variable in this map. The legend shows the continuous range of values in the data from around 300 to 420 meters.

Some rasters contain categorical data where each pixel represents a discrete class such as a landcover type (e.g., "forest" or "grassland") rather than a continuous value such as elevation or temperature. Some examples of classified maps include:

1. Landcover/land-use maps.
2. Tree height maps classified as short, medium, tall trees.
3. Elevation maps classified as low, medium and high elevation.

Categorical Landcover Map for the United States

Categorical Elevation Map of the NEON Harvard Forest Site

The legend of this map shows the colors representing each discrete class.

# add a color map with 5 colors

col=terrain.colors(3)

# add breaks to the colormap (4 breaks = 3 segments)

brk <- c(250,350, 380,500)



# Expand right side of clipping rect to make room for the legend

par(xpd = FALSE,mar=c(5.1, 4.1, 4.1, 4.5))

# DEM with a custom legend

plot(DSM_HARV, 
	col=col, 
	breaks=brk, 
	main="Classified Elevation Map - HARV",
	legend = FALSE
	)



# turn xpd back on to force the legend to fit next to the plot.

par(xpd = TRUE)

# add a legend - but make it appear outside of the plot

legend( 733100, 4713700,
        legend = c("High Elevation", "Middle","Low Elevation"), 
        fill = rev(col))

What is a GeoTIFF??

Raster data can come in many different formats. In this tutorial, we will use the geotiff format which has the extension .tif. A .tif file stores metadata or attributes about the file as embedded tif tags. For instance, your camera might store a tag that describes the make and model of the camera or the date the photo was taken when it saves a .tif. A GeoTIFF is a standard .tif image format with additional spatial (georeferencing) information embedded in the file as tags. These tags can include the following raster metadata:

1. A Coordinate Reference System (CRS)
2. Spatial Extent (extent)
3. Values that represent missing data (NoDataValue)
4. The resolution of the data

In this tutorial we will discuss all of these metadata tags.

More about the .tif format:

• GeoTIFF on Wikipedia
• OSGEO TIFF documentation

Coordinate Reference System

The Coordinate Reference System or CRS tells R where the raster is located in geographic space. It also tells R what method should be used to "flatten" or project the raster in geographic space.

What Makes Spatial Data Line Up On A Map?

There are many great resources that describe coordinate reference systems and projections in greater detail (read more, below). For the purposes of this activity, it is important to understand that data from the same location but saved in different projections will not line up in any GIS or other program. Thus, it's important when working with spatial data in a program like R to identify the coordinate reference system applied to the data and retain it throughout data processing and analysis.

Read More:

• A comprehensive online library of CRS information.
• QGIS Documentation - CRS Overview.
• Choosing the Right Map Projection.
• NCEAS Overview of CRS in R.

How Map Projections Can Fool the Eye

Check out this short video, from Buzzfeed, highlighting how map projections can make continents seems proportionally larger or smaller than they actually are!

View Raster Coordinate Reference System (CRS) in R

We can view the CRS string associated with our R object using thecrs() method. We can assign this string to an R object, too.

# view crs description

crs(DSM_HARV,describe=TRUE)

##                    name authority  code
## 1 WGS 84 / UTM zone 18N      EPSG 32618
##                                                                                                                                                                                                                                                          area
## 1 Between 78°W and 72°W, northern hemisphere between equator and 84°N, onshore and offshore. Bahamas. Canada - Nunavut; Ontario; Quebec. Colombia. Cuba. Ecuador. Greenland. Haiti. Jamaica. Panama. Turks and Caicos Islands. United States (USA). Venezuela
##            extent
## 1 -78, -72, 84, 0

# assign crs to an object (class) to use for reprojection and other tasks

harvCRS <- crs(DSM_HARV)

The CRS of our DSM_HARV object tells us that our data are in the UTM projection, in zone 18N.

The CRS in this case is in a char format. This means that the projection information is strung together as a series of text elements.

We'll focus on the first few components of the CRS, as described above.

• name: The projection of the dataset. Our data are in WGS84 (World Geodetic System 1984) / UTM (Universal Transverse Mercator) zone 18N. WGS84 is the datum. The UTM projection divides up the world into zones, this element tells you which zone the data are in. Harvard Forest is in Zone 18.
• authority: EPSG (European Petroleum Survey Group) - organization that maintains a geodetic parameter database with standard codes
• code: The EPSG code. For more details, see EPSG 32618.

Extent

The spatial extent is the geographic area that the raster data covers.

The spatial extent of an R spatial object represents the geographic "edge" or location that is the furthest north, south, east and west. In other words, extent represents the overall geographic coverage of the spatial object.

You can see the spatial extent using terra::ext:

ext(DSM_HARV)

## SpatExtent : 732000, 733000, 4713000, 4714000 (xmin, xmax, ymin, ymax)

Resolution

A raster has horizontal (x and y) resolution. This resolution represents the area on the ground that each pixel covers. The units for our data are in meters. Given our data resolution is 1 x 1, this means that each pixel represents a 1 x 1 meter area on the ground.

The best way to view resolution units is to look at the coordinate reference system string crs(rast,proj=TRUE). Notice our data contains: +units=m.

crs(DSM_HARV,proj=TRUE)

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

Display Raster Min and Max Values

It can be 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. However if they weren't already calculated, we can calculate them using the min() or max() functions.

# view the min and max values

min(DSM_HARV)

## class       : SpatRaster 
## dimensions  : 1000, 1000, 1  (nrow, ncol, nlyr)
## resolution  : 1, 1  (x, y)
## extent      : 732000, 733000, 4713000, 4714000  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 18N (EPSG:32618) 
## source(s)   : memory
## varname     : NEON_D01_HARV_DP3_732000_4713000_DSM 
## name        :    min 
## min value   : 317.91 
## max value   : 433.94

We can see that the elevation at our site ranges from 317.91 m to 433.94 m.

NoData Values in Rasters

Raster data often has a NoDataValue associated with it. This is a value assigned to pixels where data are missing or no data were collected.

By default the shape of a raster is always square or rectangular. So if we have a dataset that has a shape that isn't square or rectangular, some pixels at the edge of the raster will have NoDataValues. This often happens when the data were collected by an airplane which only flew over some part of a defined region.

Let's take a look at some of the RGB Camera data over HARV, this time downloading a tile at the edge of the flight box.

byTileAOP(dpID='DP3.30010.001',

          site='HARV',

          year='2022',

          easting=737500,

          northing=4701500,

          check.size=FALSE, # set to TRUE or remove if you want to check the size before downloading

          savepath = wd)

This file will be downloaded into a nested subdirectory under the ~/data folder, inside a folder named DP3.30010.001 (the Camera Data Product ID). The file should show up in this location: ~/data/DP3.30010.001/neon-aop-products/2022/FullSite/D01/2022_HARV_7/L3/Camera/Mosaic/2022_HARV_7_737000_4701000_image.tif.

In the image below, the pixels that are black have NoDataValues. The camera did not collect data in these areas.

# Use rast function to read in all bands

RGB_HARV <- 
  rast(paste0(wd,"DP3.30010.001/neon-aop-products/2022/FullSite/D01/2022_HARV_7/L3/Camera/Mosaic/2022_HARV_7_737000_4701000_image.tif"))



# Create an RGB image from the raster

par(col.axis="white",col.lab="white",tck=0)

plotRGB(RGB_HARV, r = 1, g = 2, b = 3, axes=TRUE)

In the next image, the black edges have been assigned NoDataValue. R doesn't render pixels that contain a specified NoDataValue. R assigns missing data with the NoDataValue as NA.

# reassign cells with 0,0,0 to NA

func <- function(x) {
  x[rowSums(x == 0) == 3, ] <- NA
  x}



newRGBImage <- app(RGB_HARV, func)

## 

|---------|---------|---------|---------|

par(col.axis="white",col.lab="white",tck=0)

# Create an RGB image from the raster stack

plotRGB(newRGBImage, r = 1, g = 2, b = 3, axis=TRUE)

## Warning in plot.window(...): "axis" is not a graphical parameter

## Warning in plot.xy(xy, type, ...): "axis" is not a graphical parameter

## Warning in title(...): "axis" is not a graphical parameter

NoData Value Standard

The assigned NoDataValue varies across disciplines; -9999 is a common value used in both the remote sensing field and the atmospheric fields. It is also the standard used by the National Ecological Observatory Network (NEON).

If we are lucky, our GeoTIFF file has a tag that tells us what is the NoDataValue. If we are less lucky, we can find that information in the raster's metadata. If a NoDataValue was stored in the GeoTIFF tag, when R opens up the raster, it will assign each instance of the value to NA. Values of NA will be ignored by R as demonstrated above.

Bad Data Values in Rasters

Bad data values are different from NoDataValues. Bad data values are values that fall outside of the applicable range of a dataset.

Examples of Bad Data Values:

• The normalized difference vegetation index (NDVI), which is a measure of greenness, has a valid range of -1 to 1. Any value outside of that range would be considered a "bad" value.
• Reflectance data in an image should range from 0-1 (or 0-10,000 depending upon how the data are scaled). Thus a value greater than 1 or greater than 10,000 is likely caused by an error in either data collection or processing. These erroneous values can occur, for example, in water vapor absorption bands, which contain invalid data, and are meant to be disregarded.

Find Bad Data Values

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

Create A Histogram of Raster Values

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

# view histogram of data

hist(DSM_HARV,
     main="Distribution of Digital Surface Model Values\n NEON Harvard Forest (HARV)",
     xlab="DSM Elevation Value (m)",
     ylab="Frequency",
     col="lightblue")

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.

Raster Bands

The Digital Surface Model object (DSM_HARV) that we've been working with is a single band raster. This means that there is only one dataset stored in the raster: surface elevation in meters for one time period.

A raster dataset can contain one or more bands. We can use the rast()function to import all bands from a single OR multi-band raster. We can view the number of bands in a raster using thenlyr()` function.

# view number of bands in the Lidar DSM raster

nlyr(DSM_HARV)

## [1] 1

# view number of bands in the RGB Camera raster

nlyr(RGB_HARV)

## [1] 3

As we see from the RGB camera raster, raster data can also be multi-band, meaning one raster file contains data for more than one variable or time period for each cell. By default the terra::rast() function imports all bands of a multi-band raster. You can set lyrs = 1 if you only want to read in the first layer, for example.

Jump to the fourth tutorial in this series for a tutorial on multi-band rasters: Work with Multi-band Rasters: Images in R.

View Raster File Attributes

Remember that a GeoTIFF contains a set of embedded tags that contain metadata about the raster. So far, we've explored raster metadata after importing it in R. However, we can use the describe("path-to-raster-here") function to view raster information (such as metadata) before we open a file in R. Use help(describe) to see other options for exploring the file contents.

# view metadata attributes before opening the file

describe(path.expand(dsm_harv_file),meta=TRUE)

##  [1] "AREA_OR_POINT=Area"                                                                                                                                                                                                                                                                                   
##  [2] "TIFFTAG_ARTIST=Created by the National Ecological Observatory Network (NEON)"                                                                                                                                                                                                                         
##  [3] "TIFFTAG_COPYRIGHT=The National Ecological Observatory Network is a project sponsored by the National Science Foundation and managed under cooperative agreement by Battelle. This material is based in part upon work supported by the National Science Foundation under Grant No. DBI-0752017."      
##  [4] "TIFFTAG_DATETIME=Flown on 2022080312, 2022080412, 2022081213, 2022081413"                                                                                                                                                                                                                             
##  [5] "TIFFTAG_IMAGEDESCRIPTION=Elevation LiDAR - NEON.DP3.30024 acquired at HARV by RIEGL LASER MEASUREMENT SYSTEMS Q780 2220855 as part of 2022-P3C1"                                                                                                                                                      
##  [6] "TIFFTAG_MAXSAMPLEVALUE=434"                                                                                                                                                                                                                                                                           
##  [7] "TIFFTAG_MINSAMPLEVALUE=318"                                                                                                                                                                                                                                                                           
##  [8] "TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)"                                                                                                                                                                                                                                                               
##  [9] "TIFFTAG_SOFTWARE=Tif file created with a Matlab script (write_gtiff.m) written by Tristan Goulden (tgoulden@battelleecology.org) with data processed from the following scripts: create_tiles_from_mosaic.m, combine_dtm_dsm_gtif.m, lastools_workflow.csh which implemented LAStools version 210418."
## [10] "TIFFTAG_XRESOLUTION=1"                                                                                                                                                                                                                                                                                
## [11] "TIFFTAG_YRESOLUTION=1"

Specifying options=c("stats") will show some summary statistics:

# view summary statistics before opening the file

describe(path.expand(dsm_harv_file),options=c("stats"))

##  [1] "Driver: GTiff/GeoTIFF"                                                                                                                                                                                                                                                                                  
##  [2] "Files: C:/Users/bhass/Documents/data/DP3.30024.001/neon-aop-products/2022/FullSite/D01/2022_HARV_7/L3/DiscreteLidar/DSMGtif/NEON_D01_HARV_DP3_732000_4713000_DSM.tif"                                                                                                                                   
##  [3] "       C:/Users/bhass/Documents/data/DP3.30024.001/neon-aop-products/2022/FullSite/D01/2022_HARV_7/L3/DiscreteLidar/DSMGtif/NEON_D01_HARV_DP3_732000_4713000_DSM.tif.aux.xml"                                                                                                                           
##  [4] "Size is 1000, 1000"                                                                                                                                                                                                                                                                                     
##  [5] "Coordinate System is:"                                                                                                                                                                                                                                                                                  
##  [6] "PROJCRS[\"WGS 84 / UTM zone 18N\","                                                                                                                                                                                                                                                                     
##  [7] "    BASEGEOGCRS[\"WGS 84\","                                                                                                                                                                                                                                                                            
##  [8] "        ENSEMBLE[\"World Geodetic System 1984 ensemble\","                                                                                                                                                                                                                                              
##  [9] "            MEMBER[\"World Geodetic System 1984 (Transit)\"],"                                                                                                                                                                                                                                          
## [10] "            MEMBER[\"World Geodetic System 1984 (G730)\"],"                                                                                                                                                                                                                                             
## [11] "            MEMBER[\"World Geodetic System 1984 (G873)\"],"                                                                                                                                                                                                                                             
## [12] "            MEMBER[\"World Geodetic System 1984 (G1150)\"],"                                                                                                                                                                                                                                            
## [13] "            MEMBER[\"World Geodetic System 1984 (G1674)\"],"                                                                                                                                                                                                                                            
## [14] "            MEMBER[\"World Geodetic System 1984 (G1762)\"],"                                                                                                                                                                                                                                            
## [15] "            MEMBER[\"World Geodetic System 1984 (G2139)\"],"                                                                                                                                                                                                                                            
## [16] "            ELLIPSOID[\"WGS 84\",6378137,298.257223563,"                                                                                                                                                                                                                                                
## [17] "                LENGTHUNIT[\"metre\",1]],"                                                                                                                                                                                                                                                              
## [18] "            ENSEMBLEACCURACY[2.0]],"                                                                                                                                                                                                                                                                    
## [19] "        PRIMEM[\"Greenwich\",0,"                                                                                                                                                                                                                                                                        
## [20] "            ANGLEUNIT[\"degree\",0.0174532925199433]],"                                                                                                                                                                                                                                                 
## [21] "        ID[\"EPSG\",4326]],"                                                                                                                                                                                                                                                                            
## [22] "    CONVERSION[\"UTM zone 18N\","                                                                                                                                                                                                                                                                       
## [23] "        METHOD[\"Transverse Mercator\","                                                                                                                                                                                                                                                                
## [24] "            ID[\"EPSG\",9807]],"                                                                                                                                                                                                                                                                        
## [25] "        PARAMETER[\"Latitude of natural origin\",0,"                                                                                                                                                                                                                                                    
## [26] "            ANGLEUNIT[\"degree\",0.0174532925199433],"                                                                                                                                                                                                                                                  
## [27] "            ID[\"EPSG\",8801]],"                                                                                                                                                                                                                                                                        
## [28] "        PARAMETER[\"Longitude of natural origin\",-75,"                                                                                                                                                                                                                                                 
## [29] "            ANGLEUNIT[\"degree\",0.0174532925199433],"                                                                                                                                                                                                                                                  
## [30] "            ID[\"EPSG\",8802]],"                                                                                                                                                                                                                                                                        
## [31] "        PARAMETER[\"Scale factor at natural origin\",0.9996,"                                                                                                                                                                                                                                           
## [32] "            SCALEUNIT[\"unity\",1],"                                                                                                                                                                                                                                                                    
## [33] "            ID[\"EPSG\",8805]],"                                                                                                                                                                                                                                                                        
## [34] "        PARAMETER[\"False easting\",500000,"                                                                                                                                                                                                                                                            
## [35] "            LENGTHUNIT[\"metre\",1],"                                                                                                                                                                                                                                                                   
## [36] "            ID[\"EPSG\",8806]],"                                                                                                                                                                                                                                                                        
## [37] "        PARAMETER[\"False northing\",0,"                                                                                                                                                                                                                                                                
## [38] "            LENGTHUNIT[\"metre\",1],"                                                                                                                                                                                                                                                                   
## [39] "            ID[\"EPSG\",8807]]],"                                                                                                                                                                                                                                                                       
## [40] "    CS[Cartesian,2],"                                                                                                                                                                                                                                                                                   
## [41] "        AXIS[\"(E)\",east,"                                                                                                                                                                                                                                                                             
## [42] "            ORDER[1],"                                                                                                                                                                                                                                                                                  
## [43] "            LENGTHUNIT[\"metre\",1]],"                                                                                                                                                                                                                                                                  
## [44] "        AXIS[\"(N)\",north,"                                                                                                                                                                                                                                                                            
## [45] "            ORDER[2],"                                                                                                                                                                                                                                                                                  
## [46] "            LENGTHUNIT[\"metre\",1]],"                                                                                                                                                                                                                                                                  
## [47] "    USAGE["                                                                                                                                                                                                                                                                                             
## [48] "        SCOPE[\"Navigation and medium accuracy spatial referencing.\"],"                                                                                                                                                                                                                                
## [49] "        AREA[\"Between 78°W and 72°W, northern hemisphere between equator and 84°N, onshore and offshore. Bahamas. Canada - Nunavut; Ontario; Quebec. Colombia. Cuba. Ecuador. Greenland. Haiti. Jamaica. Panama. Turks and Caicos Islands. United States (USA). Venezuela.\"],"                        
## [50] "        BBOX[0,-78,84,-72]],"                                                                                                                                                                                                                                                                           
## [51] "    ID[\"EPSG\",32618]]"                                                                                                                                                                                                                                                                                
## [52] "Data axis to CRS axis mapping: 1,2"                                                                                                                                                                                                                                                                     
## [53] "Origin = (732000.000000000000000,4714000.000000000000000)"                                                                                                                                                                                                                                              
## [54] "Pixel Size = (1.000000000000000,-1.000000000000000)"                                                                                                                                                                                                                                                    
## [55] "Metadata:"                                                                                                                                                                                                                                                                                              
## [56] "  AREA_OR_POINT=Area"                                                                                                                                                                                                                                                                                   
## [57] "  TIFFTAG_ARTIST=Created by the National Ecological Observatory Network (NEON)"                                                                                                                                                                                                                         
## [58] "  TIFFTAG_COPYRIGHT=The National Ecological Observatory Network is a project sponsored by the National Science Foundation and managed under cooperative agreement by Battelle. This material is based in part upon work supported by the National Science Foundation under Grant No. DBI-0752017."      
## [59] "  TIFFTAG_DATETIME=Flown on 2022080312, 2022080412, 2022081213, 2022081413"                                                                                                                                                                                                                             
## [60] "  TIFFTAG_IMAGEDESCRIPTION=Elevation LiDAR - NEON.DP3.30024 acquired at HARV by RIEGL LASER MEASUREMENT SYSTEMS Q780 2220855 as part of 2022-P3C1"                                                                                                                                                      
## [61] "  TIFFTAG_MAXSAMPLEVALUE=434"                                                                                                                                                                                                                                                                           
## [62] "  TIFFTAG_MINSAMPLEVALUE=318"                                                                                                                                                                                                                                                                           
## [63] "  TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)"                                                                                                                                                                                                                                                               
## [64] "  TIFFTAG_SOFTWARE=Tif file created with a Matlab script (write_gtiff.m) written by Tristan Goulden (tgoulden@battelleecology.org) with data processed from the following scripts: create_tiles_from_mosaic.m, combine_dtm_dsm_gtif.m, lastools_workflow.csh which implemented LAStools version 210418."
## [65] "  TIFFTAG_XRESOLUTION=1"                                                                                                                                                                                                                                                                                
## [66] "  TIFFTAG_YRESOLUTION=1"                                                                                                                                                                                                                                                                                
## [67] "Image Structure Metadata:"                                                                                                                                                                                                                                                                              
## [68] "  INTERLEAVE=BAND"                                                                                                                                                                                                                                                                                      
## [69] "Corner Coordinates:"                                                                                                                                                                                                                                                                                    
## [70] "Upper Left  (  732000.000, 4714000.000) ( 72d10'28.52\"W, 42d32'36.84\"N)"                                                                                                                                                                                                                              
## [71] "Lower Left  (  732000.000, 4713000.000) ( 72d10'29.98\"W, 42d32' 4.46\"N)"                                                                                                                                                                                                                              
## [72] "Upper Right (  733000.000, 4714000.000) ( 72d 9'44.73\"W, 42d32'35.75\"N)"                                                                                                                                                                                                                              
## [73] "Lower Right (  733000.000, 4713000.000) ( 72d 9'46.20\"W, 42d32' 3.37\"N)"                                                                                                                                                                                                                              
## [74] "Center      (  732500.000, 4713500.000) ( 72d10' 7.36\"W, 42d32'20.11\"N)"                                                                                                                                                                                                                              
## [75] "Band 1 Block=1000x1 Type=Float32, ColorInterp=Gray"                                                                                                                                                                                                                                                     
## [76] "  Min=317.910 Max=433.940 "                                                                                                                                                                                                                                                                             
## [77] "  Minimum=317.910, Maximum=433.940, Mean=358.584, StdDev=17.156"                                                                                                                                                                                                                                        
## [78] "  NoData Value=-9999"                                                                                                                                                                                                                                                                                   
## [79] "  Metadata:"                                                                                                                                                                                                                                                                                            
## [80] "    STATISTICS_MAXIMUM=433.94000244141"                                                                                                                                                                                                                                                                 
## [81] "    STATISTICS_MEAN=358.58371301653"                                                                                                                                                                                                                                                                    
## [82] "    STATISTICS_MINIMUM=317.91000366211"                                                                                                                                                                                                                                                                 
## [83] "    STATISTICS_STDDEV=17.156044149253"                                                                                                                                                                                                                                                                  
## [84] "    STATISTICS_VALID_PERCENT=100"

It can be useful to use describe to explore your file before reading it into R.