Spatial analysis can be extremely helpful when it comes to determining what is happening in a region. In the case of our work we use it to analyze what is happening in rivers and watersheds in Southern Virginia. After metric calculations have been done to determine flow we can create images that show us regions with high or low flow and how that can change over time.

This tutorial is meant to step you through accessing shapefiles and other data that is used to create maps and perform spatial analysis of a desired study area. It also explains how to adjust the colors, shape and size of the images you are producing.

Download Libraries

R has libraries that let us do alot of amazing things. We have to turn these libraries on though for them to work. To create images in R we will need the following libraries. If you do not have these libraries already installed in your RStudio simply do install.packages(' '). Be sure that the package name is in quotes and in the parenthesis.

library(stringr)
library(rapportools)
library(tidyr)
library(questionr)
library(rgdal)
library(raster)
library(randomcoloR)

Determining Map Projection

Before creating a map you have to determine what projection you want it to be in so that all of your components will be projected the same. This requires a Proj4 code.Search the projection you want (in this case WGS 1984 was used becuase it is compatible with VA Hydro) and choose the link for ’spatialreference.org.’Once on the website click on the Proj4 link and copy and paste the provided line in quotes after its associated variable name; like so.

Projection<- '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs'

Accessing Local Data

Now that we know our projection let’s find the folder that you have stored this tutorial in so that we can get our data. Odds are it is in your downloads folder, and the filepath looks something like this: “C:\\Users\\Username\\Downloads\\Arc_GIS”.

Within that folder are two shapefiles and one csv file. We are going to open these files in RStudio using the readOGR() and read.csv() commands. Go ahead and open a new script in R.

Copy Entire File Path

First we’ll open the necessary shapefiles one at a time. This means that we are going to paste in the entire path name where the shapefile is stored. Copy the following line into your script, replacing the filepath with yours. Remember to include either double backslashes (\\) or single forward slashes (/) within the path.

Let’s start with state boundaries. So the way this works is you have to paste the entire file path to where the shapefile is located followed by a comma and the actual shapefile name. We then have to set the projection so that the data is in a form that we can work with. To do that we use spTransform()in the parenthesis we put the shapefile variable that we just read in followed by CRS= the projection.

State <-readOGR('C:\\Users\\Person\\Downloads\\Arc_GIS\\states', 'state')  
State<- spTransform(State, CRS=Projection)

Now we just need to repeat this process for the River Segments shapefile and metrics csv.

RivSeg<-readOGR ('C:\\Users\\Person\\Downloads\\Arc_GIS\\BaseMap', 'RiverSegs')            
RivSeg<-spTransform(RivSeg, CRS=Projection) 
Metrics<- read.csv('C:\\Users\\Person\\Downloads\\Arc_GIS\\Percent Error All Metrics.csv')


Prepping Data so Shapefile and Metrics Can be Combined

First we will specify which metric we want to map spatially. We will then need to create an empty column for the desired data that will be applied to the river segments. The metrics will be added to the shapefile for the river segments as seprate columns and the desired metric will be given its own column. To do so we set up the code just like we are working with any normal vector but we have to include an @data between the variable name and and column name that we want to reference.

DesiredMetric<- 'AugLowFlow'
RivSeg@data$RiverSeg<-as.character(RivSeg@data$RiverSeg)
RivSeg@data$Metric<-NA

Metrics[,paste0("DesiredMetric")] <- Metrics[DesiredMetric]


Combining Shapefile and Metric

To do this we will need to create a loop that will look at every river segment and determine if it has an associated metric. The if statement in the loop is asking if the river segment ID is in the metrics file, if so make it true and go inside the if statement, if not skip the if statement and progress in the loop.

for (i in 1:length(RivSeg@data$RiverSeg)){
  if (RivSeg@data$RiverSeg[i]%in%Metrics$river_seg){
    RivSeg@data$Metric[i]<- Metrics$DesiredMetric[Metrics$river_seg==RivSeg@data$RiverSeg[i]]
  }
}



Creating the Map

Now for the fun part, we have all of the data we need and it has been formatted properly so now we make the map. What we are going to do first is create the path that we want the map saved to (it is much eaiser to just go ahead and save the map as a png than to try and export it later).

folder_location<- 'C:\\Users\\Person\\Downloads\\Arc_GIS'
dir.create(paste0(folder_location,"\\OUTPUT\\SpatialComparisons"), showWarnings = FALSE);
png(filename=paste0(folder_location,"\\OUTPUT\\SpatialComparisons\\", DesiredMetric ,"- Error.png"), 
    width=1400, height=950, units="px")

Now we will pick the colors and groupings that we will put our river segments in. The way this works, is R puts the metrics in caetgories based on number so we have to set those categories and then assign a color to them as follows;

RivSeg@data$color<- cut(RivSeg@data$Metric,c(-Inf,-50,-20,20,50,Inf), labels=c('red', 'palegreen', 'green4', 'yellowgreen', 'orange'))

What this has done is it says our first category or bin is anything from -infinity to -50, anything that falls in that category will be red. The next bin is -50 to -20 with anything in that group being palegreen, and so on. Now that we have our colors we will cut our data to only show the riversegments that have been put in a bin. 

SouthernRivers<- RivSeg[!is.na (RivSeg@data$color),]

Everything should be set to go now so it is time to plot! The order in which you plot things is the order that they are layered. So we plot the ‘SouthernRivers’ twice, once to get the size of the plot we want and then to plot the rivers over the state boundaries.

title<-paste0((DesiredMetric),' Percent Error')
plot(SouthernRivers, col=paste0(SouthernRivers@data$color))
plot(State, add=TRUE, col='gray')
lines(State, col='white')
plot(SouthernRivers, col=paste0(SouthernRivers@data$color), add=TRUE)
legend("bottomright", legend=c('< -50', '-50 to -20', '-20 to 20', '20 to 50', '> 50'), col=c('red', 'palegreen', 'green4', 'yellowgreen', 'orange'), lty=0, pch=15, pt.cex=2.5, bty='n', y.intersp=0.75, x.intersp=0.3, cex=1.5, lwd=1)
legend("topleft", legend=c(title), lty=0, pt.cex=3, bty='n', y.intersp=0.75, x.intersp=0.3, cex=2, lwd=1)

dev.off()

Here is a list of a few of the things you can do to modify this plot (Go ahead and try some of these on your own)

Plot Command What it does Ways to change it
lwd=.9 changes line weight change number
cex=1 changes point size or font size plot: change number, legend: change scale
legend(’’) changes legend position specify a position or coordinates
lty=1 changes line type in legend change number for different types
type=‘l’ specifies line type p, l, b, c, o, h, n
col=‘blue’ changes line color type colors() in console for full list
pch=15 creates box for legend change number to change shape