Introduction

In this exercise, we will download and extract future precipitation data from the MACA dataset and calculate changes in future precipitation relative to historic precipitation for the Big Sandy watershed.

Specifically, we analyze monthly precipitation for the US for a historic period (2000-2004) and future period (2096-2099) for the Relative Concentration Pathway (RCP) 8.5. RCP 8.5 is characterized by increasing greenhouse gas emissions over time, representative of scenarios in the literature that lead to high greenhouse gas concentration levels (Riahi et al. 2007). For more information on RCPs, see https://skepticalscience.com/rcp.php.

Future climate in this exercise is based on a single GCM, the Norwegian Earth System Model - NorESM model. NorESM is a global, coupled model system for the physical climate system, which can be run with various degrees of interactions with bio-geo-chemical processes in the earth system (https://folk.uib.no/ngfhd/EarthClim/index.htm#no).

The NorESM is just one GCM but the extraction and analysis approach presented in this exercise is applicable to any of teh 20 GCMs contained in the MACA dataset.

MACA

MACA stands for Multivariate Adaptive Constructed Analogs, MACA is a statistical method for downscaling Global Climate Models (GCMs) from their native coarse resolution to a higher spatial resolution that captures reflects observed patterns of daily near-surface meteorology and simulated changes in GCMs experiments. This method has been shown to be slightly preferable to direct daily interpolated bias correction in regions of complex terrain due to its use of a historical library of observations and multivariate approach.

MACA includes two different downscaled datasets covering CONUS-plus. Both use a common set of 20 CMIP5 GCMs from models that provided daily output of requisite variables for historical (1950-2005) and future experiments under RCP4.5 and RCP8.5.

These data are available at different spatial and temporal scales for the United States and includes maximum temperature, minimum temperature, maximum and minimum relative humidity, precipitation accumulation, downward surface shortwave radiation, wind-velocity, and specific humidity.

You can learn more about MACA at http://www.climatologylab.org/maca.html

Load package(s)

For this exercise, we will use the packages called below.

library(raster)
## Loading required package: sp
library(lubridate)
## 
## Attaching package: 'lubridate'
## The following objects are masked from 'package:raster':
## 
##     intersect, union
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union
library(sf)
## Linking to GEOS 3.5.1, GDAL 2.2.2, PROJ 4.9.2
library(tidyverse)
## ── Attaching packages ────────────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.3.0     ✓ purrr   0.3.4
## ✓ tibble  3.0.1     ✓ dplyr   0.8.5
## ✓ tidyr   1.0.2     ✓ stringr 1.4.0
## ✓ readr   1.3.1     ✓ forcats 0.5.0
## ── Conflicts ───────────────────────────────────────────────── tidyverse_conflicts() ──
## x lubridate::as.difftime() masks base::as.difftime()
## x lubridate::date()        masks base::date()
## x tidyr::extract()         masks raster::extract()
## x dplyr::filter()          masks stats::filter()
## x lubridate::intersect()   masks raster::intersect(), base::intersect()
## x dplyr::lag()             masks stats::lag()
## x dplyr::select()          masks raster::select()
## x lubridate::setdiff()     masks base::setdiff()
## x lubridate::union()       masks raster::union(), base::union()
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

In addition, we will need to install and load the ‘ncd4’ package, which is used to interface with netcdf files, a common file structure for storing climate data.

Install ‘ncd4’ using ‘Install’ udner the ‘Packages’ tab in the File Browser window and load

library(ncdf4)

Preloaded datafiles

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

  1. watershed boundary polygon - “rockville_guage_wshed.shp.shp”
  2. MACA monthly precipitation data:
  • macav2metdata_pr_NorESM1-M_r1i1p1_historical_2000_2004_CONUS_monthly - historic (2000-2004)
  • macav2metdata_pr_NorESM1-M_r1i1p1_rcp85_2096_2099_CONUS_monthly - future (2096-2099) RCP 8.5

Download and Analyze MACA climate data

Data were dowloaded from the MACA website using the instructions in “Download_Instructions_Maca_Climate_Data_5-15-2020.docx” file

Define projection, load shapes, and rasters

# Define a projection and load watershed boundary shapefile

prj<-"+proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0"

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

# create a shapefile object you want to extract from the netcdfs, so it have the same extent.
shp<-shapefile(BigSandy) 

# reproject shape
shp_rp<-spTransform(shp,CRS(prj))

Get a list of netcdf files and check the name of the variable in the netcdf file

# obtain names of files to use, copying the name of the file into R
pr_historic<-"macav2metdata_pr_NorESM1-M_r1i1p1_historical_2000_2004_CONUS_monthly (1).nc"
pr_future<-"macav2metdata_pr_NorESM1-M_r1i1p1_rcp85_2096_2099_CONUS_monthly.nc"

# list names in dataset (use historic as example)
raster(pr_historic)
## [1] "vobjtovarid4: error #F: I could not find the requsted var (or dimvar) in the file!"
## [1] "var (or dimvar) name: crs"
## [1] "file name: /cloud/project/macav2metdata_pr_NorESM1-M_r1i1p1_historical_2000_2004_CONUS_monthly (1).nc"
## class      : RasterLayer 
## band       : 1  (of  60  bands)
## dimensions : 13, 25, 325  (nrow, ncol, ncell)
## resolution : 0.0416654, 0.04166603  (x, y)
## extent     : 279.998, 281.0396, 39.45868, 40.00034  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
## source     : /cloud/project/macav2metdata_pr_NorESM1-M_r1i1p1_historical_2000_2004_CONUS_monthly (1).nc 
## names      : Monthly.Precipitation.Amount 
## z-value    : 2000-01-15 
## zvar       : precipitation
# zvar is precipitation, the name of the variable that you want to extract in the netcdf files.

Define function to extract data from MACA files

#function to extract the netcdf: in this case the precipitation
extract_fun<-function(netcdf){ 
  big_raster<-stack(netcdf,varname="precipitation", quick=TRUE) # gets 365 days into stacked layers
  big_raster@extent[1]<-big_raster@extent[1]-360 # for unknown reason, longitutde is incorrect, this   should fix it. 
  big_raster@extent[2]<-big_raster@extent[2]-360
  crs(big_raster)<-prj #reproject the raster 
  sum_raster<-sum(big_raster)/5 # Calculates annual precipitation by adding all the monthly dimensions
  return(sum_raster)# the return value of the function
}

Apply the function and save rasters

# Create raster by applying functions: 
pr_2000_2004<-pr_historic %>% extract_fun() 
## [1] "vobjtovarid4: error #F: I could not find the requsted var (or dimvar) in the file!"
## [1] "var (or dimvar) name: crs"
## [1] "file name: /cloud/project/macav2metdata_pr_NorESM1-M_r1i1p1_historical_2000_2004_CONUS_monthly (1).nc"
pr_2096_2099<-pr_future %>% extract_fun() 
## [1] "vobjtovarid4: error #F: I could not find the requsted var (or dimvar) in the file!"
## [1] "var (or dimvar) name: crs"
## [1] "file name: /cloud/project/macav2metdata_pr_NorESM1-M_r1i1p1_rcp85_2096_2099_CONUS_monthly.nc"
# Create raster as *.tif
writeRaster(x = pr_2000_2004,filename = "mean_precip_historic.tif",overwrite=T)
writeRaster(x=pr_2096_2099,filename="mean_precip_future.tif",overwrite=T)

Read saved rasters, extract precipitation using shapefile, and spatially average data

raster_mean<-raster("mean_precip_historic.tif") 

# Return spatially-averaged mean annual precipitation for historic period  
Big_sandy_historic<-raster::extract(pr_2000_2004,shp_rp,fun=mean) 

# Return spatially-averaged mean annual precipitation for future period  
Big_sandy_future<-raster::extract(pr_2096_2099,shp_rp,fun=mean)

Calculate change in future precipitation relative to historic period and plot

difference<-(Big_sandy_future-Big_sandy_historic)/Big_sandy_historic *100
difference
##           [,1]
## [1,] -13.72069
#The future P in 2096-2099 according to the Norwegian model for RCP 8.5 results in 13.82% less annual precipitation for the Big Sandy. 
brks<-c(700,800,900,1000,1100,1200,1300,1400,1500,1600)

par(mfrow=c(2,1), mar=c(5,5,3,0))
plot(pr_2000_2004, col=blues9, breaks=brks, main="Historic Precipitation",xlab="Lon",ylab="Lat")
plot(shp_rp,add=T)
plot(pr_2096_2099,col=blues9, breaks=brks, main="Future Precipitation", xlab="Lon",ylab="Lat")
plot(shp_rp,add=T)