Search code examples
rrasterr-sp

Extracting values from one of multiple PRISM rasters in R


I have a set of 100 lat/longs and 100 (sometimes repeating) years. For each of the lat/longs, I would like to find the annual temperature (using PRISM data) in the corresponding year.

I've created for loops and functions to iterate through the process, but they keep getting stuck. I know there must be a simpler way of doing this (preferably avoiding a for loop!).

As an example, here are 4 lat/longs and years, and my (laughable) attempt to loop over them.

library(prism)
library(raster)

locs <- data.frame( lat = c(46.30101, 42.65503, 44.38075, 43.90637), lon = c(-91.764380  -86.201983, -88.951511, -91.081340, -87.896017))

years <- c(1989,1954,2010,1954)

coordinates(locs) <- c('lat', 'lon')

temps <- NULL

for(i in 1:length(years)) {
  tryCatch({dir.create(paste0(getwd(),"/",years[i]))}, error=function(e){}) # skip making a new directory for any years that already exist
  options(prism.path = paste0(getwd(),"/",years[i]))
  get_prism_annual(type = "tmean", years = as.numeric(years[i])) # Get the data
  climate_data <- prism_stack(ls_prism_data()) # Stack it
  climate_crs <- climate_data@crs@projargs # Get the projection
  proj4string(occ.latlong) <- CRS(climate_crs) # Project the climate data's CRS onto the coordinates
  
  temps <- rbind(temps, extract(climate_data, locs[i,]))
}

This loop throws NAs for repeat years (1954 above). Is there an easier way to do this?!


Solution

  • Because you are downloading a lot of data, I would split the loop into two parts, and first do the downloads. In that way it is also easier to debug your code. I also simplified the code were possible.

    years <- c(1989, 1954, 2010, 1954)
    
    library(prism)
    for (y in unique(years)) {
        dir.create(y, FALSE)
        options(prism.path = y)
        get_prism_annual(type = "tmean", years = y) 
    }
    

    Now extract for the locations.

    locs <- data.frame( lat = c(46.30101, 42.65503, 44.38075, 43.90637), lon = c(-91.764380  -86.201983, -88.951511, -91.081340, -87.896017))
    # lon/lat order!    
    locs <- locs[,2:1]
    

    Keep things simple and put the results in a list first

    temps <- list()
    for(i in 1:length(years)) {
        options(prism.path = years[i])
        climate_data <- prism_stack(ls_prism_data()) 
        temps[[i]] <- extract(climate_data, locs)
    }
    

    Combine the results

    x <- do.call(cbind, temps)
    

    Note that it can be rather inefficient to extract the same values (same years) multiple times, rather than reusing the extracted values.