For example i tried this but it doesn't work


it returns

Error in compareRaster(e1, e2, extent = FALSE, rowcol = FALSE, crs = TRUE,  : 
different origin
  • You are attempting to add the rasters. You need to mosaic them using the raster packages "mosaic()" function. Apr 6 '15 at 18:58
  • at what point of the above code should i use the mosaic function?Checked the help but encountered difficulties. Elaborate
    – gsa
    Apr 6 '15 at 19:17
  • You should not be forcing the extent and resolution using "res" and "extent". Just rasterize, using your reference raster "rast_a", which I assume represents the entire extent and desired resolution, and then mosaic the results. Apr 6 '15 at 20:03
  • i cant really grasp what you are saying since without the "res", "rast_a" doesn't have the resolution info. Also, i have three shps that i have to combine after i turn them into raster and i don't know how to do it with mosaic. An example would be greatly appreciated
    – gsa
    Apr 6 '15 at 21:00
  • Generally, the OP has to provide a reproducible example, not the opposite.
    – user32309
    Apr 7 '15 at 3:20

Here's a simple reproducible example. I give 2 approaches, one using mosaic/merge and another that just does the initial rasterize at the total extent of all shapefiles combined. The results are the same.

Added based on Jeffrey Evans comment: You also need to consider what to do when you have overlapping polygons. If you want to apply a function (e.g. sum, mean, min, max) to the input rasters to get the output use mosaic. If you want the output to be based on the order of input rasters, use merge.

# Example polygons
p1 <- rbind(c(-180,-20), c(-140,55), c(10, 0), c(-140,-60), c(-180,-20))
hole <- rbind(c(-150,-20), c(-100,-10), c(-110,20), c(-150,-20))
p1 <- SpatialPolygons(list(Polygons(list(Polygon(p1), Polygon(hole, hole=TRUE)), 1)))
p2 <- SpatialPolygons(list(Polygons(list(Polygon(rbind(c(-10,0), c(140,60), c(160,0), c(140,-55), c(-10,0)))), 2)))
p3 <- SpatialPolygons(list(Polygons(list(Polygon(rbind(c(-125,0), c(0,60), c(40,5), c(15,-45), c(-125,0)))), 3)))
# Add some fake data
p1 <- SpatialPolygonsDataFrame(p1, data=data.frame(val=1), match.ID=F)
p2 <- SpatialPolygonsDataFrame(p2, data=data.frame(val=2), match.ID=F)
p3 <- SpatialPolygonsDataFrame(p3, data=data.frame(val=3), match.ID=F)

# Method 1: rasterize at full extent of all shapefile
e <- union(union(extent(p1), extent(p2)), extent(p3))
r <- raster(e, res=1)
r1a <- rasterize(p1, r, field='val')
r2a <- rasterize(p2, r, field='val')
r3a <- rasterize(p3, r, field='val')
rsum <- sum(r1a, r2a, r3a, na.rm=T)

plot(e, axes=F, xlab='', ylab='')
plot(p1, add=T, col='red')
plot(p2, add=T, col='blue')
plot(p3, add=T, col='green')

# Method 2: rasterize at extents of individual shapefiles
r1b <- rasterize(p1, raster(extent(p1), res=1), field='val')
r2b <- rasterize(p2, raster(extent(p2), res=1), field='val')
r3b <- rasterize(p3, raster(extent(p3), res=1), field='val')
# Mosaic together, summing values that overlap
# Change the fun parameter to mean, min, or max as necessary
rMosaic <- mosaic(r1b, r2b, r3b, fun=sum, na.rm=T)

# Check that they're the same
cellStats( ((rMosaic - rsum) != 0), sum)

# Method 3: if you want outcome of an overlap to be determined by order instead of a function
# In the case of an overlap, the value is assigned based on order of arguments in merge
rMerge <- merge(r1a, r2a, r3a)

Using mosiac and summing input rasters:
Using merge:

  • Nice job. One more option is the merge function in raster. This is the same as mosaic but does not require a function and orders overlap values on order of input. For the most part you do not want to sum rasters to mosaic them. In your example, overlap values represent the sum of the overlap which, is not always desirable. The difference is a max value of 6 vs a mean of 3. For mosaic operations you often want the mean, min or max. Contrary to the OP's example, I am not confident that they actually want the sum of the overlap. Apr 7 '15 at 22:50
  • Thanks, I wasn't really clear what the difference between mosaic and merge was. I've edited the original post to incorporate both these functions. Thanks!
    – Matt SM
    Apr 7 '15 at 23:17
  • @MattSM i need clarifications, why you need the rbind?how many vectors you have inside the binds and why those numbers. Too complicated.Lots of questions.i cant implement it on my own problem. Thanks
    – gsa
    Apr 9 '15 at 15:16
  • @gsa I'm just generating some made up data to create a reproducible example. The bit of code you posted isn't very useful because you don't provide any information on what shp_1, shp_2, and shp_3 are. You'll want to start with the line e <- union(union(extent(p1), extent(p2)), extent(p3)) and just replace p1, p2, and p3 with your shapefiles (shp_1, shp_2, shp_3)
    – Matt SM
    Apr 9 '15 at 21:48
  • okay,now its more clear. the "val" field what is it supposed to contain? and is it obligatory to add a "field" parameter when we rasterize? May i guess that in this position we put a "field that is common between all the shapefiles?
    – gsa
    Apr 10 '15 at 9: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.