0

I am attempting to merge many rasters into a single layer based on year. I have a more complicated version of the merge answer found here Merge rasters with different origins in R and I basically have the same problem as the original question as that post.

The difference for me though is that I have set up a script that works in some cases but in others I get an error. For example the rasters for year 1985 work fine but in 1986 I get the following error

Error in compareRaster(x, extent = FALSE, rowcol = FALSE, orig = TRUE,  : 
  different origin

I have narrowed down the culprit in 1986 to one raster and if removed the script runs fine but I cannot see what is different about this raster as compared to others in the dataset.

Here's the best reprex I could come up with. In this example when I try to merge r1 and r2 I get the above error (r1 is the reproduced raster that I have problems in 1986). When I merge r2 and r3 I get a completely different error

Error in .local(.Object, ...) : 

but I didn't have the patience to figure that out and we'll just pretend that it works for the sake of my question.

r1 <- new("RasterLayer", file = new(".RasterFile", name = "Fagin", datanotation = "FLT4S", byteorder = "little", nodatavalue = -Inf, NAchanged = FALSE, nbands = 1L, bandorder = "BIL", offset = 0L, toptobottom = TRUE, blockrows = 256L, blockcols = 256L, driver = "gdal", open = FALSE), data = new(".SingleLayerData", values = logical(0), offset = 0, gain = 1, inmemory = FALSE, fromdisk = TRUE, isfactor = FALSE, attributes = list(), haveminmax = FALSE, min = Inf, max = -Inf, band = 1L, unit = "", names = "Fagan"), legend = new(".RasterLegend", type = character(0), values = logical(0), color = logical(0), names = logical(0), colortable = logical(0)), title = character(0), extent = new("Extent", xmin = -110.785720561, xmax = -110.785181572, ymin = 31.840515751, ymax = 31.840785246), rotated = FALSE, rotation = new(".Rotation", geotrans = numeric(0), transfun = function () NULL), ncols = 2L, nrows = 1L, crs = new("CRS", projargs = "+proj=longlat +datum=WGS84 +no_defs"), history = list(), z = list())

r2 <- new("RasterLayer", file = new(".RasterFile", name = "E2", datanotation = "FLT4S", byteorder = "little", nodatavalue = -Inf, NAchanged = FALSE, nbands = 1L, bandorder = "BIL", offset = 0L, toptobottom = TRUE, blockrows = 256L, blockcols = 256L, driver = "gdal", open = FALSE), data = new(".SingleLayerData", values = logical(0), offset = 0, gain = 1, inmemory = FALSE, fromdisk = TRUE, isfactor = FALSE, attributes = list(), haveminmax = FALSE, min = Inf, max = -Inf, band = 1L, unit = "", names = "E2"), legend = new(".RasterLegend", type = character(0), values = logical(0), color = logical(0), names = logical(0), colortable = logical(0)), title = character(0), extent = new("Extent", xmin = -107.282829942, xmax = -107.279865502, ymin = 34.325255827, ymax = 34.327411784), rotated = FALSE, rotation = new(".Rotation", geotrans = numeric(0), transfun = function () NULL), ncols = 11L, nrows = 8L, crs = new("CRS", projargs = "+proj=longlat +datum=WGS84 +no_defs"), history = list(), z = list())

r3 <- new("RasterLayer", file = new(".RasterFile", name = "Blank", datanotation = "FLT4S", byteorder = "little", nodatavalue = -Inf, NAchanged = FALSE, nbands = 1L, bandorder = "BIL", offset = 0L, toptobottom = TRUE, blockrows = 256L, blockcols = 256L, driver = "gdal", open = FALSE), data = new(".SingleLayerData", values = logical(0), offset = 0, gain = 1, inmemory = FALSE, fromdisk = TRUE, isfactor = FALSE, attributes = list(), haveminmax = FALSE, min = Inf, max = -Inf, band = 1L, unit = "", names = "Blank"), legend = new(".RasterLegend", type = character(0), values = logical(0), color = logical(0), names = logical(0), colortable = logical(0)), title = character(0), extent = new("Extent", xmin = -105.624899254, xmax = -105.621395824, ymin = 35.682161064, ymax = 35.684317021), rotated = FALSE, rotation = new(".Rotation", geotrans = numeric(0), transfun = function () NULL), ncols = 13L, nrows = 8L, crs = new("CRS", projargs = "+proj=longlat +datum=WGS84 +no_defs"), history = list(), z = list())

l1 <- list(r1, r2)
l2 <- list(r2, r3)

tmp <- do.call(raster::merge, l1)
tmp2 <- do.call(raster::merge, l2)

So, what is it about the r1 raster that messes up my merge when I have no problem in other cases with similar data?

1 Answer 1

1

With RasterLayer x you can use as.character(x) to create representations like below that can be used for examples.

library(raster)
r1 <- raster(ncols=2, nrows=1, xmn=-110.785720561, xmx=-110.785181572, ymn=31.840515751, ymx=31.840785246, crs='+proj=longlat +datum=WGS84 +no_defs')
r2 <- raster(ncols=11, nrows=8, xmn=-107.282829942, xmx=-107.279865502, ymn=34.325255827, ymx=34.327411784, crs='+proj=longlat +datum=WGS84 +no_defs')
r3 <- raster(ncols=13, nrows=8, xmn=-105.624899254, xmx=-105.621395824, ymn=35.682161064, ymx=35.684317021, crs='+proj=longlat +datum=WGS84 +no_defs')

m <- merge(r1, r2)
#Error in compareRaster(x, extent = FALSE, rowcol = FALSE, orig = TRUE,  : 
#  different origin

But this works

m <- merge(r1, r2, tolerance=0.15)

This is what is happening under the hood:

compareRaster(r1, r2, extent=FALSE, rowcol=FALSE, orig=TRUE, res=TRUE, tolerance=0.1)
#Error in compareRaster(r1, r2, extent = FALSE, rowcol = FALSE, orig = TRUE,  : 
#   different origin
compareRaster(r2, r3, extent=FALSE, rowcol=FALSE, orig=TRUE, res=TRUE, tolerance=0.05)
# [1] TRUE

compareRaster does:

o1 <- abs(origin(r1))
o2 <- abs(origin(r2))

dif <- o1 - o2
minres <- min(res(r1))
tol <- 0.05
all.equal(dif, c(0,0), tolerance=tol, scale=minres)
#[1] "Mean scaled difference: 0.1171468"

I am not saying this is how it should be, but this is what causes the behavior.

(and when I see raster data like this, I always wonder what caused this mess in the first place. In most of such cases the workflow can be fixed upstream, by assuring that all rasters have an origin at zero).

3
  • So I considered adjusting the tolerance to correct the issue but I am not sure how to implement this when using do.call(raster::merge, l1). I have about 40 to 200 rasters per year with 35 years of data so I don't really want to go through it one at a time.
    – Mike D
    Commented May 20, 2021 at 3:17
  • Also, out of curiosity why do you think the rasters are a mess? They were created based on a polygons overlaid with Landsat data. I just used dput() to get a reproducible example because that's what I've used before. Seemed to create a very detailed reproduction just in case it was something else causing the problem.
    – Mike D
    Commented May 20, 2021 at 3:37
  • 1
    You can do l1$tolerance = 0.15. I think they are messy, because they do not perfectly align, don't have an origin of zero. Probably because they were projected without using a template. Commented May 20, 2021 at 4:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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