Subtracting rasters after gradient field calculation in R fails?

Question

I am attempting to write a code in R that imports 3 rasters with the same extent/projection, calculates the horizontal and vertical gradient fields (rise over run in cells), then subtract the outputs from each other (for example, after horizontal gradient fields are calculated for each of the three rasters, the horizontal gradient output for raster 1 is subtracted from the horizontal gradient output for raster 2). I’m using this code on multiple "triplets" of rasters. Sometimes this code works for some triplets, but not for others. This would suggest the problem is with the input rasters, and not the code, but when I inspect the rasters that cause the code to fail, they appear no different than other rasters that the code runs through successfully.

For rasters when the code does not work, I first get the warning message:

Warning message:
In .rasterFromRasterFile(grdfile, band = band, objecttype, ...) :
size of values file does not match the number of cells (given the 
data type)

Then, when I try to subtract the vertical rasters from each other (this problem only happens for subtracting vertical gradient field rasters, not horizontal), I get this error message:

Error in setValues(r, methods::callGeneric(getValues(x))) : 
length(values) is not equal to ncell(x), or to 1

It’s important to note I only get this error when trying to find the absolute value of the differences; if I just subtract the two rasters I get the warning, but not the error message.

Below is the code I have written for this.

library(raster)

#Load rasters
raster1<- raster(*input path to raster 1*)
raster2<- raster(*input path to raster 2*)
raster3<- raster(*input path to raster 3*)

#Calculate vertical (i) and horizontal (j) gradient fields of raster 
#files. Warning message occurs for calculation of vertical gradient
#fields, but not horizontal. 

zone1_igrad1<- focal(zone1_raster11, matrix(c(-1/2,0,1/2), ncol = 3))
zone1_jgrad1<- focal(zone1_raster11, matrix(c(-1/2,0,1/2), nrow = 3))
zone1_igrad2<- focal(zone1_raster21, matrix(c(-1/2,0,1/2), ncol = 3))
zone1_jgrad2<- focal(zone1_raster21, matrix(c(-1/2,0,1/2), nrow = 3))
zone1_igrad3<- focal(zone1_raster31, matrix(c(-1/2,0,1/2), ncol = 3))
zone1_jgrad3<- focal(zone1_raster31, matrix(c(-1/2,0,1/2), nrow = 3))

#Subtract Horizontal Gradient Field Outputs. 
zone1_hori1<- abs(zone1_jgrad1-zone1_jgrad2)
zone1_hori2<- abs(zone1_jgrad1-zone1_jgrad3)
zone1_hori3<- abs(zone1_jgrad2-zone1_jgrad3)
zone1_hori_field<- zone1_hori1 + zone1_hori2 + zone1_hori3

#Subtract Vertical Gradient Field Outputs. (Error message occurs 
#here)
zone1_vert1<- abs(zone1_igrad1-zone1_igrad2)
zone1_vert2<- abs(zone1_igrad1-zone1_igrad3)
zone1_vert3<- abs(zone1_igrad2-zone1_igrad3)
zone1_vert_field<- zone1_vert1 + zone1_vert2 + zone1_vert3

Answer 1

It looks like this is a bug that occurs when you use a one column/row weights matrix and when the results are written to disk (if the dataset is large, to a temp file). I will fix that, but for now I think you can work around by making a 3x3 matrix (with added zeros)

library(raster)
# force writing to disk, even with small datasets
rasterOptions(todisk=TRUE)

f <- system.file("external/rlogo.grd", package="raster")
raster1<- raster(f,1)
raster2<- raster(f,2)
raster3<- raster(f,3)

m <- matrix(c(-1/2,0,1/2))

# fix the problem: 
m1 <- cbind(0, m, 0)
m2 <- rbind(0, t(m), 0)


zone1_igrad1 <- focal(raster1, m1)
zone1_jgrad1 <- focal(raster1, m2)
zone1_igrad2 <- focal(raster2, m1)
zone1_jgrad2 <- focal(raster2, m2)
zone1_igrad3 <- focal(raster3, m1)
zone1_jgrad3 <- focal(raster3, m2)

zone1_hori1 <- abs(zone1_jgrad1-zone1_jgrad2)
zone1_hori2 <- abs(zone1_jgrad1-zone1_jgrad3)
zone1_hori3 <- abs(zone1_jgrad2-zone1_jgrad3)
zone1_hori_field <- zone1_hori1 + zone1_hori2 + zone1_hori3    
zone1_vert1 <- abs(zone1_igrad1-zone1_igrad2)
zone1_vert2 <- abs(zone1_igrad1-zone1_igrad3)
zone1_vert3 <- abs(zone1_igrad2-zone1_igrad3)
zone1_vert_field <- zone1_vert1 + zone1_vert2 + zone1_vert3