9.3 PRISM with prism | R as GIS for Economists

9.3 PRISM with `prism`

9.3.1 Basics

PRISM dataset provide model-based estimates of precipitation, tmax, and tmin for the U.S. at the 4km by 4km spatial resolution. Here, we use get_prism_dailys() from the prism package (Hart and Bell 2015) to download daily data. Here is its general syntax:

#--- NOT RUN ---#
get_prism_dailys(
  type = variable type,
  minDate = starting date as character,
  maxDate = ending date as character,
  keepZip = TRUE or FALSE
)

The variables types you can select from is “ppt” (precipitation), “tmean” (mean temperature), “tmin” (minimum temperature), and “tmax” (maximum temperature). For minDate and maxDate, the dates must be specified in a specific format of “YYYY-MM-DD.” keepZip = FALSE does not keep the zipped folders of the downloaded files as the name suggests.

Before you download PRISM data using the function, it is recommended that you set the path to folder in which the downloaded PRISM will be stored using options(prism.path = "path"). For example, the following set the path to “Data/PRISM/” relative to the current working directory.

options(prism.path = "Data/PRISM/")

The following code downloads daily precipitation data from January 1, 2000 to Jan 10, 2000.

#--- NOT RUN ---#
get_prism_dailys(
  type = "ppt",
  minDate = "2000-01-01",
  maxDate = "2000-01-10",
  keepZip = FALSE
)

When you download data using the above code, you will notice that it creates one folder for one day. For example, for precipitation data for “2000-01-01,” you can get the path to the downloaded file as follows:

var_type <- "ppt" # variable type
dates_prism_txt <- str_remove_all("2000-01-01", "-") # date without dashes

#--- folder name ---#
folder_name <- paste0("PRISM_", var_type, "_stable_4kmD2_", dates_prism_txt, "_bil") 

#--- file name of the downloaded data inside the above folder ---#
file_name <- paste0("PRISM_", var_type, "_stable_4kmD2_", dates_prism_txt, "_bil.bil") 

#--- path to the file relative to the designated data folder (here, it's "Data/PRISM/") ---#
(
file_path <- paste0("Data/PRISM/", folder_name, "/", file_name)
)

[1] "Data/PRISM/PRISM_ppt_stable_4kmD2_20000101_bil/PRISM_ppt_stable_4kmD2_20000101_bil.bil"

We can then easily read the data using terra::rast() or stars::read_stars() if you prefer the stars way.

#--- as SpatRaster ---#
(
prism_2000_01_01_sr <- rast(file_path) 
)

class       : SpatRaster 
dimensions  : 621, 1405, 1  (nrow, ncol, nlyr)
resolution  : 0.04166667, 0.04166667  (x, y)
extent      : -125.0208, -66.47917, 24.0625, 49.9375  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat NAD83 
source      : PRISM_ppt_stable_4kmD2_20000101_bil.bil 
name        : PRISM_ppt_stable_4kmD2_20000101_bil 
min value   :                                   0 
max value   :                              49.848

#--- as stars ---#
(
prism_2000_01_01_stars <- read_stars(file_path) 
)

stars object with 2 dimensions and 1 attribute
attribute(s):
                                   Min. 1st Qu. Median      Mean 3rd Qu.   Max.
PRISM_ppt_stable_4kmD2_2000010...     0       0      0 0.4952114       0 49.848
                                     NA's
PRISM_ppt_stable_4kmD2_2000010...  390874
dimension(s):
  from   to   offset      delta refsys point values x/y
x    1 1405 -125.021  0.0416667  NAD83    NA   NULL [x]
y    1  621  49.9375 -0.0416667  NAD83    NA   NULL [y]

Here is a quick visualization of the data (Figure 9.2):

plot(prism_2000_01_01_stars)

Figure 9.2: PRISM precipitation data for January 1, 2000

As you can see, the data covers the entire contiguous U.S.

9.3.2 Download daily PRISM data for many years and build your own datasets

Here, an example of how to create your own sets of PRISM datasets is presented. Creating such datasets and have them locally can be useful if you expect to use the data for many different projects in the future.

Suppose we are interested in saving daily PRISM precipitation data by year-month from 1980 to 2018. We will write a loop that loops over all the year-month combinations. Before writing a loop, let’s work on the code for a particular year-month combination: 1990-December. However, we will write codes in a way that can be easily translated into a looped operations later. Specifically, we will define the following variables and use them as if they are variables to be looped over in a loop.

#--- month to work on ---#  
temp_month <- 12

#--- year to work on ---#  
temp_year <- 1990

We first need to set the path to the folder in which daily PRISM files will be downloaded.

#--- set your own path ---#
options(prism.path = "Data/PRISM/")

We then set the start and end dates for get_prism_dailys().

#--- starting date of the working month-year ---#
(
start_date <- dmy(paste0("1/", temp_month, "/", temp_year))
)

[1] "1990-12-01"

#--- ending date: add a month and then go back 1 day ---#
(
end_date <- start_date %m+% months(1) - 1
)

[1] "1990-12-31"

We now download PRISM data for the year-month we are working on.

#--- download daily PRISM data for the working month-year ---#
get_prism_dailys(
  type = "ppt",
  minDate = as.character(start_date),
  maxDate = as.character(end_date),
  keepZip = FALSE
)

Once all the data are downloaded, we will read and import them onto R. To do so, we will need the path to all the downloaded files.

#--- list of dates of the working month-year ---#
dates_ls <- seq(start_date, end_date, "days") 

#--- remove dashes ---#
dates_prism_txt <- str_remove_all(dates_ls, "-")

#--- folder names ---#
folder_name <- paste0("PRISM_", var_type, "_stable_4kmD2_", dates_prism_txt, "_bil") 
#--- the file name of the downloaded data ---#
file_name <- paste0("PRISM_", var_type, "_stable_4kmD2_", dates_prism_txt, "_bil.bil") 
#--- complete path to the downloaded files ---#
(
file_path <- paste0("Data/PRISM/", folder_name, "/", file_name) 
)

 [1] "Data/PRISM/PRISM_ppt_stable_4kmD2_19901201_bil/PRISM_ppt_stable_4kmD2_19901201_bil.bil"
 [2] "Data/PRISM/PRISM_ppt_stable_4kmD2_19901202_bil/PRISM_ppt_stable_4kmD2_19901202_bil.bil"
 [3] "Data/PRISM/PRISM_ppt_stable_4kmD2_19901203_bil/PRISM_ppt_stable_4kmD2_19901203_bil.bil"
 [4] "Data/PRISM/PRISM_ppt_stable_4kmD2_19901204_bil/PRISM_ppt_stable_4kmD2_19901204_bil.bil"
 [5] "Data/PRISM/PRISM_ppt_stable_4kmD2_19901205_bil/PRISM_ppt_stable_4kmD2_19901205_bil.bil"
 [6] "Data/PRISM/PRISM_ppt_stable_4kmD2_19901206_bil/PRISM_ppt_stable_4kmD2_19901206_bil.bil"
 [7] "Data/PRISM/PRISM_ppt_stable_4kmD2_19901207_bil/PRISM_ppt_stable_4kmD2_19901207_bil.bil"
 [8] "Data/PRISM/PRISM_ppt_stable_4kmD2_19901208_bil/PRISM_ppt_stable_4kmD2_19901208_bil.bil"
 [9] "Data/PRISM/PRISM_ppt_stable_4kmD2_19901209_bil/PRISM_ppt_stable_4kmD2_19901209_bil.bil"
[10] "Data/PRISM/PRISM_ppt_stable_4kmD2_19901210_bil/PRISM_ppt_stable_4kmD2_19901210_bil.bil"
[11] "Data/PRISM/PRISM_ppt_stable_4kmD2_19901211_bil/PRISM_ppt_stable_4kmD2_19901211_bil.bil"
[12] "Data/PRISM/PRISM_ppt_stable_4kmD2_19901212_bil/PRISM_ppt_stable_4kmD2_19901212_bil.bil"
[13] "Data/PRISM/PRISM_ppt_stable_4kmD2_19901213_bil/PRISM_ppt_stable_4kmD2_19901213_bil.bil"
[14] "Data/PRISM/PRISM_ppt_stable_4kmD2_19901214_bil/PRISM_ppt_stable_4kmD2_19901214_bil.bil"
[15] "Data/PRISM/PRISM_ppt_stable_4kmD2_19901215_bil/PRISM_ppt_stable_4kmD2_19901215_bil.bil"
[16] "Data/PRISM/PRISM_ppt_stable_4kmD2_19901216_bil/PRISM_ppt_stable_4kmD2_19901216_bil.bil"
[17] "Data/PRISM/PRISM_ppt_stable_4kmD2_19901217_bil/PRISM_ppt_stable_4kmD2_19901217_bil.bil"
[18] "Data/PRISM/PRISM_ppt_stable_4kmD2_19901218_bil/PRISM_ppt_stable_4kmD2_19901218_bil.bil"
[19] "Data/PRISM/PRISM_ppt_stable_4kmD2_19901219_bil/PRISM_ppt_stable_4kmD2_19901219_bil.bil"
[20] "Data/PRISM/PRISM_ppt_stable_4kmD2_19901220_bil/PRISM_ppt_stable_4kmD2_19901220_bil.bil"
[21] "Data/PRISM/PRISM_ppt_stable_4kmD2_19901221_bil/PRISM_ppt_stable_4kmD2_19901221_bil.bil"
[22] "Data/PRISM/PRISM_ppt_stable_4kmD2_19901222_bil/PRISM_ppt_stable_4kmD2_19901222_bil.bil"
[23] "Data/PRISM/PRISM_ppt_stable_4kmD2_19901223_bil/PRISM_ppt_stable_4kmD2_19901223_bil.bil"
[24] "Data/PRISM/PRISM_ppt_stable_4kmD2_19901224_bil/PRISM_ppt_stable_4kmD2_19901224_bil.bil"
[25] "Data/PRISM/PRISM_ppt_stable_4kmD2_19901225_bil/PRISM_ppt_stable_4kmD2_19901225_bil.bil"
[26] "Data/PRISM/PRISM_ppt_stable_4kmD2_19901226_bil/PRISM_ppt_stable_4kmD2_19901226_bil.bil"
[27] "Data/PRISM/PRISM_ppt_stable_4kmD2_19901227_bil/PRISM_ppt_stable_4kmD2_19901227_bil.bil"
[28] "Data/PRISM/PRISM_ppt_stable_4kmD2_19901228_bil/PRISM_ppt_stable_4kmD2_19901228_bil.bil"
[29] "Data/PRISM/PRISM_ppt_stable_4kmD2_19901229_bil/PRISM_ppt_stable_4kmD2_19901229_bil.bil"
[30] "Data/PRISM/PRISM_ppt_stable_4kmD2_19901230_bil/PRISM_ppt_stable_4kmD2_19901230_bil.bil"
[31] "Data/PRISM/PRISM_ppt_stable_4kmD2_19901231_bil/PRISM_ppt_stable_4kmD2_19901231_bil.bil"

We now read them as a stars object, set the third dimension as date using Dates object class, and then save it as an R dataset (This will ensure that date dimensions is kept. See Chapter 7.4).

(
#--- combine all the PRISM files as stars ---#
temp_stars <- read_stars(file_path, along = 3)
  #--- set the third dimension as data ---# 
  st_set_dimensions("band", values = dates_ls, name = "date")
)

#--- save the stars as an rds file ---#
saveRDS(
  temp_stars, 
  paste0("Data/PRISM/PRISM_", var_type, "_y", temp_year, "_m", temp_month, ".rds")
)

You could alternatively read the files into a SpatRaster object and save it data as a GeoTIFF file.

(
#--- combine all the PRISM files as a RasterStack ---#
temp_stars <- terra::rast(file_path) 
)

#--- save as a multi-band GeoTIFF file ---#
writeRaster(temp_stars, paste0("Data/PRISM/PRISM_", var_type, "_y", temp_year, "_m", temp_month, ".tif"), overwrite = T)

Note that this option of course does not have date as the third dimension. Moreover, the RDS file above takes up only 14 Mb, while the tif file occupies 108 Mb.

Finally, if you would like, you can delete all the individual PRISM files:

#--- delete all the downloaded files ---#
unlink(paste0("Data/PRISM/", folder_name), recursive = TRUE)

Okay, now that we know what to do with a particular year-month combination, we can easily write a loop to go over all the year-month combinations for the period of interest. Since all the processes we observed above for a single year-month combination is embarrassingly parallel, it is easy to parallelize using future.apply::future_lapply() or parallel::mclapply() (Linux/Mac users only). Here we use future_lapply(). Let’s first get the number of logical cores.

library(parallel)
num_cores <- detectCores() 

plan(multiprocess, workers = num_cores)

The following function goes through all the steps we saw above for a single year-month combination.

#--- define a function to download and save PRISM data stacked by month ---#
get_save_prism <- function(i, var_type) {

  print(paste0("working on ", i))

  temp_month <- month_year_data[i, month] # working month
  temp_year <- month_year_data[i, year] # working year
  
  #--- starting date of the working month-year ---#
  start_date <- dmy(paste0("1/", temp_month, "/", temp_year))
  #--- end date ---#
  end_date <- start_date %m+% months(1) - 1

  
  #--- download daily PRISM data for the working month-year ---#
  get_prism_dailys(
    type = var_type,
    minDate = as.character(start_date),
    maxDate = as.character(end_date),
    keepZip = FALSE
  ) 

  #--- list of dates of the working month-year ---#
  dates_ls <- seq(start_date, end_date, "days") 

  #--- remove dashes ---#
  dates_prism_txt <- str_remove_all(dates_ls, "-")

  #--- folder names ---#
  folder_name <- paste0("PRISM_", var_type, "_stable_4kmD2_", dates_prism_txt, "_bil") 
  #--- the file name of the downloaded data ---#
  file_name <- paste0("PRISM_", var_type, "_stable_4kmD2_", dates_prism_txt, "_bil.bil") 
  #--- complete path to the downloaded files ---#
  file_path <- paste0("Data/PRISM/", folder_name, "/", file_name)

  #--- combine all the PRISM files as a RasterStack ---#
  temp_stars <- stack(file_path) %>% 
    #--- convert to stars ---#
    st_as_stars() %>% 
    #--- set the third dimension as data ---# 
    st_set_dimensions("band", values = dates_ls, name = "date")

  #--- save the stars as an rds file ---#
  saveRDS(
    temp_stars, 
    paste0("Data/PRISM/PRISM_", var_type, "_y", temp_year, "_m", temp_month, ".rds")
  )

  #--- delete all the downloaded files ---#
  unlink(paste0("Data/PRISM/", folder_name), recursive = TRUE)
}

We then create a data.frame of all the year-month combinations:

(
#--- create a set of year-month combinations to loop over ---#
month_year_data <- expand.grid(month  = 1:12, year = 1990:2018) %>% 
  data.table()
)

     month year
  1:     1 1990
  2:     2 1990
  3:     3 1990
  4:     4 1990
  5:     5 1990
 ---           
344:     8 2018
345:     9 2018
346:    10 2018
347:    11 2018
348:    12 2018

We now do parallelized loop over all the year-month combinations (by looping over the rows of the month_year_data):

#--- run the above code in parallel ---#
future_lapply(
  1:nrow(month_year_data), 
  function (x) get_save_prism(x, "ppt")
)

That’s it. Of course, you can do the same thing for tmax by this:

#--- run the above code in parallel ---#
future_lapply(
  1:nrow(month_year_data), 
  function (x) get_save_prism(x, "tmax")
)

Now that you have PRISM datasets, you can extract values from the raster layers for vector data for your analysis, which is covered extensively in Chapters 5 and 7 (for stars objects).

If you want to save the data by year (each file would be about 168 Mb). You could do this.

#--- define a function to download and save PRISM data stacked by year ---#
get_save_prism_y <- function(temp_year, var_type) {

  print(paste0("working on ", temp_year))
  
  #--- starting date of the working month-year ---#
  start_date <- dmy(paste0("1/1/", temp_year))
  #--- end date ---#
  end_date <- dmy(paste0("1/1/", temp_year + 1)) - 1

  #--- download daily PRISM data for the working month-year ---#
  get_prism_dailys(
    type = var_type,
    minDate = as.character(start_date),
    maxDate = as.character(end_date),
    keepZip = FALSE
  ) 

  #--- list of dates of the working month-year ---#
  dates_ls <- seq(start_date, end_date, "days") 

  #--- remove dashes ---#
  dates_prism_txt <- str_remove_all(dates_ls, "-")

  #--- folder names ---#
  folder_name <- paste0("PRISM_", var_type, "_stable_4kmD2_", dates_prism_txt, "_bil") 
  #--- the file name of the downloaded data ---#
  file_name <- paste0("PRISM_", var_type, "_stable_4kmD2_", dates_prism_txt, "_bil.bil") 
  #--- complete path to the downloaded files ---#
  file_path <- paste0("Data/PRISM/", folder_name, "/", file_name)

  #--- combine all the PRISM files as a RasterStack ---#
  temp_stars <- stack(file_path) %>% 
    #--- convert to stars ---#
    st_as_stars() %>% 
    #--- set the third dimension as data ---# 
    st_set_dimensions("band", values = dates_ls, name = "date")

  #--- save the stars as an rds file ---#
  saveRDS(
    temp_stars, 
    paste0("Data/PRISM/PRISM_", var_type, "_y", temp_year, ".rds")
  )

  #--- delete all the downloaded files ---#
  unlink(paste0("Data/PRISM/", folder_name), recursive = TRUE)
} 

#--- run the above code in parallel ---#
future_lapply(
  1990:2018, 
  function (x) get_save_prism_y(x, "tmax")
)

References

Hart, Edmund M., and Kendon Bell. 2015. Prism: Download Data from the Oregon Prism Project. https://doi.org/10.5281/zenodo.33663.