Introduction

In this exercise, we will download and extract Prism precipitation data and use a ‘function’ to extract and analyze long-term annual and monthly precipitation ‘Normals’ for the Big Sandy watershed.

Prism

The PRISM Climate Group, based at Oregon State University, gathers climate observations from a wide range of monitoring networks, applies quality control measures, and develops spatial climate datasets to reveal short- and long-term climate patterns. The resulting datasets incorporate a variety of modeling techniques and are available at multiple spatial/temporal resolutions, covering the period from 1895 to the present.

These data are available at 4-km and 880-m spatial resolution for the United States. While access to some Prism data requires a fee, other data, such as those used in this course, are publicaly available. Explore Prism at http://www.prism.oregonstate.edu/.

Climate ‘Normals’

What are ‘Normals’? From NOAA (https://www.ncdc.noaa.gov/news/defining-climate-normals-new-ways)- Scientists traditionally define a Climate Normal as an average over a recent 30-year period. Our most recent installment covers the period from 1981 to 2010.

Why 30 years? Close to a century ago, the International Meteorological Organization—now known as the World Metrological Organization—instructed member nations to calculate Climate Normals using 30-year periods, beginning with 1901–1930. Also, a general rule in statistics says that you need at least 30 numbers to get a reliable estimate of their mean or average. So, our scientists have traditionally defined Normals as averages over 30 years simply because that is the accepted convention—not because a 30-year average is the only logical or “right” way to define a Climate Normal.

Functions in R

What is a function? A function is a set of statements organized together to perform a specific task. R has a large number of in-built functions and the user can create their own functions (https://www.tutorialspoint.com/r/r_functions.htm). In R, a function is an object so the R interpreter is able to pass control to the function, along with arguments that may be necessary for the function to accomplish the actions. The function in turn performs its task and returns control to the interpreter as well as any result which may be stored in other objects.

Load package(s)

For this exercise, we will use the ‘raster’ and rgdal’ packages

library(raster)   
## Loading required package: sp
library(rgdal)
## rgdal: version: 1.4-8, (SVN revision 845)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.2.2, released 2017/09/15
##  Path to GDAL shared files: /usr/share/gdal/2.2
##  GDAL binary built with GEOS: TRUE 
##  Loaded PROJ.4 runtime: Rel. 4.9.2, 08 September 2015, [PJ_VERSION: 492]
##  Path to PROJ.4 shared files: (autodetected)
##  Linking to sp version: 1.4-1

Preloaded datafiles

The following files have been uploaded to the Cloud working directory:

  1. watershed boundary polygon - “rockville_guage_wshed.shp.shp”
  2. Prism climate data - Both long-term annual (denoted with ‘annual’ in file name) and monthly (denoted by months ‘1-12’ in file name) precipitation normals have been preloaded into the Cloud file window.

Download and Analyzing Prism climate data

Prism climate data can be downloaded direclty from the Prism webpage. For this course, we are intrested in Climate Normals which have been downloaded as *.bil files. Both long-term annual (denoed with ‘annual’ in file name) and monthly (denoted by months ‘1-12’ in file name) precipitation normals have been preloaded into the Cloud file window.

Part 1: Watershed-level ANNUAL extraction

The Prism dataset includes gridded 4-km raster data for the entire US. In this section, we will extact long-term annual precipitation data using the Big Sandy watershed boundary.

Step 1. Import raster from file directory and plot

rasnam<-"PRISM_ppt_30yr_normal_4kmM2_annual_bil.bil" 
ras<-raster(rasnam) # identify as raster
prj<-toString(crs(ras)) #get raster projection
plot(ras)

Step 2. Create a function to extract the mean annual precipitation

mean_annual_precipitation_function<-function(name){ # the input should be the name of the shapefile
  shp<-shapefile(name) #convert to shapefile object
  shp_rp<-spTransform(shp,CRS(prj)) #reproject shape
  pavgyr<-extract(ras,shp_rp,fun=mean) #get spatially averaged 30-year mean annual rainfall
  return(pavgyr)
}

Step 3. Extract mean annual precipitation for the Big Sandy watershed

# Call watershed boundary into work session
BigSandy<-"rockville_guage_wshed.shp" 

# Pass data contained in watershed boundary through function
BigSandy_precip<-mean_annual_precipitation_function(BigSandy) 

# Print results of function 
BigSandy_precip
##          [,1]
## [1,] 1257.708
# Export results to .csv
write.csv(BigSandy_precip,"BigSandy_annual_average_P.csv")

Part II: Watershed-level MONTHLY extraction

Step 1. Call watershed boundary shapefile, determine raster cells inside of shapefile, and get geographic coordinates

# Shapefile name in file directory
shpfn<-"rockville_guage_wshed.shp" 
shp<-shapefile(shpfn)

# Reproject shapefile based on raster projection
shp_rp<-spTransform(shp,CRS(prj)) 

grs<-matrix(unlist(cellFromPolygon(ras,shp_rp)),byrow=TRUE) 
ng<-length(grs) #determine number of cells in basin

# Determine raster cells inside shapefile
# Get the geographic coordinates in cells inside basin
geocoords<-xyFromCell(ras,grs) 

# Determine the resolution of the dataset
res<-geocoords[2,1]-geocoords[1,1] 

# Generate empty matrix to store gridded monthly data
pbgmn<-matrix(data=NA,nrow=ng,ncol=12) 

# Generate empty matrix to store spatially averaged monthly data
pavgmn<-matrix(data=NA,nrow=12,ncol=1)

Step 2. Using a loop, analyze gridded monthly data

for (m in 1:12){ #start month cycle
  
  if (m<10){ #conditional to build month to character
    
    mnchar<-paste("0",m,sep="") #add a 0 to months below October
    
  } else { #conditional for months larger than October
    
    mnchar<-toString(m) #months from October to December convert directly to character
    
  } #end of conditional to write month as character
  
  # Generates name of monthly rasters
    crasnam<-paste("PRISM_ppt_30yr_normal_4kmM2_",mnchar,"_bil.bil",sep="") 
  cras<-raster(crasnam) #open monthly raster
  pbgmn[,m]<-extract(cras,geocoords) #get monthly values by grid
  pavgmn[m,1]<-mean(pbgmn[,m]) #calculate average basin precipitation by month
  
} #end of month cycle

Step 3. Export extracted data as .csv

# Monthly average precipitation grid
write.csv(pbgmn,file="Bigsandy_monthly_precip_grid.csv") 

# Monthly average precipitation 
write.csv(pavgmn,file="Bigsandy_monthly_precip_average.csv") #save local cell numbers