1

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
  • The raster layers you are trying to perform the focal analysis on are not defined in your example, so it's hard to say what might be going on. What it the result when you run compareRaster(raster1, raster2, raster3, res=T, orig=T, showwarning=T) – Jacob F Sep 21 '18 at 14:39
2

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
  • Thank you for your response, I will attempt this code when I am in a position to do so and follow up. One clarification regarding your comment that the above code created the same result twice, the difference between the code is if nrow = 3 or ncol = 3. They produce different results. – Wincow Sep 21 '18 at 17:30
  • of course, I have removed my comment and fixed my solution – Robert Hijmans Sep 21 '18 at 17:44
  • Apologies for the late reply, I was away from my workstation for a few days. This code worked like a charm. If you don't mind clarifying how one part of the code works I'd really appreciate it. Specifically: what does the "t" in rbind(0,t(m),0) do? Thanks again for your help – Wincow Sep 25 '18 at 15:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.