R: Creating a new raster/shape file from calculations done using variables from a shape file and a raster

Question

I am currently trying to create a new raster or shape file based on a conditional calculation that needs to be done over every value in a shape value based on a value in a raster file. I don't usually work with raster and shape files, so I am pretty out of my element here. I'm asking this in general terms, but here is the data I am using so hopefully it will give a better understanding of what I am trying to accomplish:

 rast_norm <- ftp://prism.nacse.org/normals_4km/tmean/PRISM_tmean_30yr_normal_4kmM2_04_bil.zip
 shp_probs <- ftp://ftp.cpc.ncep.noaa.gov/GIS/us_tempprcpfcst/seastemp_201603.zip

The main objective is to take the probability associated with each point (latitude and longitude) in shp_probs and multiply it by the value that corresponds to the same latitude and longitude in rast_norm, along with some other calculations afterward. If I had two data.tables, I could do something like the following:

 dt1 <- data.table(col1 = c(0:3), col2 = c(1:4)*11, factor1 = sqrt(c(285:288))

 # # Output # #
 # col1 col2  factor1
 #    0   11 16.88194
 #    1   22 16.91153
 #    2   33 16.94107
 #    3   44 16.97056

 dt2 <- data.table(col1 = c(0:3), col2 = c(1:4)*11, factor2 = abs(sin(c(1:4))))

 # # Output # #
 # col1 col2   factor1
 #    0   11 0.8414710
 #    1   22 0.9092974
 #    2   33 0.1411200
 #    3   44 0.7568025

 dt3 <- merge(dt1, dt2, by = c("col1", "col2"))
 dt3$factor1 <- dt3$factor1 * dt3$factor2
 dt3$factor2 <- NULL

 # # Output # #
 # col1 col2   factor1
 #    0   11 14.205665
 #    1   22 15.377615
 #    2   33  2.390725
 #    3   44 12.843364     

Easy-peasy using data tables. But I am at a loss trying to do this with a Raster and a SpatialPolygonsDataFrame. Here's what I have so far to read in and clean up the files:

 # Importing the "rast_norm" file, the first listed above with a link
 rast_norm <- "/my/file/path/PRISM_tmean_30yr_normal_4kmM2_04_bil.zip"
 zipdirec <- "/my/zip/directory"
 unzip(rast_norm, exdir = zipdirec)

 # Get the correct file from the file list
 rast_norm <- list.files(zipdirec, full.names = TRUE, pattern = ".bil")
 rast_norm <- rast_norm[!grepl("\\.xml", rast_norm)]

 # Convert to raster
 rast_norm <- raster(rast_norm)

Plotting rast_norm on its own gives this map.

 # Importing the "shp_probs" file, the second listed above with a link
 shp_probs <- "/my/file/path/seastemp_201603.zip"
 zipdirec <- "/my/zip/directory"
 unzip(shp_probs, exdir = zipdirec, overwrite = TRUE)

 # Get the correct file from the list of file names and find the layer name
 layer_name <- list.files(zipdirec, pattern = "lead14")
 layer_name <- layer_name[grepl(".shp", layer_name)]
 layer_name <- layer_name[!grepl("\\.xml", layer_name)]
 layer_name <- do.call("rbind", strsplit(layer_name, "\\.shp"))[,1]
 layer_name <- unique(layer_name)

 # Use the layer name to read in the shape file
 shp_probs <- readOGR(shp_probs, layer = layer_name)
 names_levels <- paste0(shp_probs$Cat, shp_probs$Prob)
 names_levels <- gsub("Below", "-", names_levels)
 names_levels <- gsub("Above", "+", names_levels)
 names_levels <- as.integer(names_levels)
 shp_probs@data$id <- names_levels
 shp_probs <- as(shp_probs, "SpatialPolygons")

 # Create a data frame of values to use in conjunction with the existing id's
 weights <- data.table(id = shp_probs$id, weight = shp_probs$id)
 weights$weight <- c(.80, .80, .10, .10, .10, .10, .10, .10, .80, .10, .10, .10, .10, .10)
 shp_probs <- SpatialPolygonsDataFrame(otlk_sp, weights, match.ID = FALSE)

Plotting shp_probs on its own gives this map.

I now want to take the probabilities that are associated with the shp_probs file and multiply it by the amounts of rainfall associated with the rast_norm file and multiply again by the weight associated with the probability in the shp_probs file.

I really don't know what to do and any help would be very much appreciated. How do I extract all of the corresponding data points for matching latitudes and longitudes? I think if I knwo that, I will know what to do.

Thank you, in advance.

Answer 1

Assuming that you want to perform this calculation for each grid cell of your raster, you can do something like this:

1. Download/read data, and add weight column. Note that here I've just used random weights, since your example seems to assign 14 weights to 7 polygons. Also, I'm not sure what purpose your id column serves, so I've skipped that part.

   library(raster)
   library(rgdal)
   
   download.file('ftp://prism.nacse.org/normals_4km/tmean/PRISM_tmean_30yr_normal_4kmM2_04_bil.zip', 
                 fr <- tempfile(), mode='wb')
   download.file('ftp://ftp.cpc.ncep.noaa.gov/GIS/us_tempprcpfcst/seastemp_201603.zip', 
                 fs <- tempfile(), mode='wb')
   
   unzip(fr, exdir=tempdir())
   unzip(fs, exdir=tempdir())
   
   r <- raster(file.path(tempdir(), 'PRISM_tmean_30yr_normal_4kmM2_04_bil.bil'))
   s <- readOGR(tempdir(), 'lead14_Apr_temp')
   s$weight <- runif(length(s))
   

2. Perform spatial overlay of the coordinates of the raster cells and the polygons. (Alternatively, you could use raster::rasterize twice to convert the Prob and id fields to rasters, and then multiplied the three rasters.)

   xy <- SpatialPoints(coordinates(r), proj4string=crs(r))
   o <- over(xy, s)
   

3. Create a new raster with the same extent/dimensions as the original raster, and assign the appropriate values to its cells.

   r2 <- raster(r)
   r2[] <- r[] * o$Prob * o$weight
   

With these random data, the result looks something like this: