How to extract data from PRISM raster for a given spatial point in R

I have a data set consisting of points with lat and long (My whole data set consists of nearly 12500 points). I want to get the daily tmin, tmax, and ppt values for these points from the PRISM rasters for the years 1981-2015.

Suppose I have only 5 data points and I want to get daily tmin and tmax for 1991-01-01 to 1991-01-31. This is what I’ve done in R. First load the required libraries

library(tidyr)
library(stringr)
library(prism)
library(raster)
library(magrittr)

and download PRISM data.

# SEE: https://github.com/ropensci/prism
# Get PRISM data
options(prism.path = "~/prismtmp")
get_prism_dailys(type="tmin", minDate = "1991-01-01", maxDate = "1991-01-31", keepZip = FALSE)
get_prism_dailys(type="tmax", minDate = "1991-01-01", maxDate = "1991-01-31", keepZip = FALSE)

After downloading we can plot PRISM data easily. For example for 1991-01-01 the tmax looks like as follows.

plot of chunk unnamed-chunk-3

We can stack up the raster files together to get a raster stack. With lots of raster files.

mystack <- ls_prism_data() %>%  prism_stack(.)  # Stack files

Now we convert our point into a spatial points data frame. We should be careful about projections. First extract projection from raster stack and then convert points to spatial points data frame.

# Get proj from raster stack
mycrs <- mystack@crs@projargs

# My points
mypoints <- data.frame(id = c(1, 2, 3, 4, 5),
                       lat = c(37.78392, 39.4987278, 39.604634, 38.7315512, 38.151546),
                       long = c(-95.6069, -95.244729, -101.5161912, -95.4363938, -97.3255883)
)

# Convert points to spatial points data frame
coordinates(mypoints) <- c('long', 'lat')
proj4string(mypoints) <- CRS(mycrs)

After adding projection to our spatial points data frame, now we can extract values from raster stack for our points with the following code.

# Extract data from raster
data <- data.frame(coordinates(mypoints), mypoints$id, extract(mystack, mypoints))

# Reshape data
data <- data %>%  gather(date, value, 4:ncol(data))

# Split date into type, year, month, and day
data$date <- gsub("PRISM_", "", data$date) %>% gsub("_stable_4kmD1_", "", .) %>% gsub("_bil", "", .)
data <- separate(data, "date", c("type", "Year", "Month", "Day"), sep = c(4, 8, 10))

# Reshape data again
data <- data %>%  spread(type, value)

# Order data
data <- data[order(data$mypoints.id),]

The data looks like this in the end.

# Print data
head(data)
##        long      lat mypoints.id Year Month Day  tmax   tmin
## 63 -95.6069 37.78392           1 1991    01  01 -3.53 -14.85
## 64 -95.6069 37.78392           1 1991    01  02  3.68  -8.95
## 65 -95.6069 37.78392           1 1991    01  03 -4.81 -13.82
## 66 -95.6069 37.78392           1 1991    01  04 -4.71 -13.96
## 67 -95.6069 37.78392           1 1991    01  05  0.87  -7.59
## 68 -95.6069 37.78392           1 1991    01  06  1.50  -5.90
tail(data)
##         long      lat mypoints.id Year Month Day  tmax   tmin
## 57 -97.32559 38.15155           5 1991    01  26 -2.22 -17.44
## 58 -97.32559 38.15155           5 1991    01  27 -0.29 -15.23
## 59 -97.32559 38.15155           5 1991    01  28  4.26  -8.01
## 60 -97.32559 38.15155           5 1991    01  29  7.65  -9.43
## 61 -97.32559 38.15155           5 1991    01  30 -7.93 -16.08
## 62 -97.32559 38.15155           5 1991    01  31  4.63 -15.02