Reputation: 79
I am writing a script that will take any three rasters, and crop them to the minimum possible extent. All three rasters will have the same resolution and projection. However, cropping to the minimum extent does not change the extents for the three rasters. I've tried setExtent and the same thing happens. If anyone can give suggestions I would really appreciate it. Here is sample code:
library(raster)
#Projection of all three rasters
newproj<- "+proj=utm +zone=4 +datum=WGS84 +units=m +no_defs +ellps=WGS84
+towgs84=0,0,0"
#Create three rasters with varying extents
raster1p<- raster(crs = newproj)
extent(raster1p)<- c(531247, 691847, 7856684, 7987884)
res(raster1p)<- c(100, 100)
values(raster1p)<- NA
raster2p<- raster(crs = newproj)
extent(raster2p)<- c(533550.8, 646550.8, 7881307, 7973807)
res(raster2p)<- c(100, 100)
values(raster2p)<- NA
raster3p<- raster(crs = newproj)
extent(raster3p)<- c(525739, 689839, 7857305, 7996505)
res(raster3p)<- c(100, 100)
values(raster3p)<- NA
#Find minimum extent
xmin1<- c(xmin(extent(raster1p)), xmin(extent(raster2p)), xmin(extent(raster3p)))
xmax1<- c(xmax(extent(raster1p)), xmax(extent(raster2p)), xmax(extent(raster3p)))
ymin1<- c(ymin(extent(raster1p)), ymin(extent(raster2p)), ymin(extent(raster3p)))
ymax1<- c(ymax(extent(raster1p)), ymax(extent(raster2p)), ymax(extent(raster3p)))
xmin_new<- min(xmin1)
xmax_new<- min(xmax1)
ymin_new<- min(ymin1)
ymax_new<- min(ymax1)
newextent=c(xmin_new, xmax_new, ymin_new, ymax_new)
#Crop rasters to minimum extent
crop(raster1p, newextent)
crop(raster2p, newextent)
crop(raster3p, newextent)
#Compare extents
extent_check<- c(extent(raster1p), extent(raster2p), extent(raster3p))
However, when I look at the extent_check to see if the extents now match, I see that the extents have not changed at all:
> extent_check
[[1]]
class : Extent
xmin : 531247
xmax : 691847
ymin : 7856684
ymax : 7987884
[[2]]
class : Extent
xmin : 533550.8
xmax : 646550.8
ymin : 7881307
ymax : 7973807
[[3]]
class : Extent
xmin : 525739
xmax : 689839
ymin : 7857305
ymax : 7996505
Any idea what I could be doing wrong? Thank you for your time
Upvotes: 1
Views: 2409
Reputation: 47146
I think it is not so much about doing something wrong, but rater a misconception (although there is a mistake in your code).
Example data
library(raster)
prj <- "+proj=utm +zone=4 +datum=WGS84"
r1 <- raster(res=100, ext=extent(c(531247, 691847, 7856684, 7987884)), crs=prj, vals=NA)
r2 <- raster(res=100, ext=extent(c(533550.8, 646550.8, 7881307, 7973807)), crs=prj, vals=NA)
r3 <- raster(res=100, ext=extent(c(525739, 689839, 7857305, 7996505)), crs=prj, vals=NA)
Find the "minimum extent"
e <- intersect(intersect(extent(r1), extent(r2)), extent(r3))
Note that the result is different from yours because you use
xmin_new <- min(xmin1)
and ymin_new <- min(ymin1)
Where it should be
xmin_new <- max(xmin1)
and ymin_new <- max(ymin1)
Now crop
r1e <- crop(r1, e)
r2e <- crop(r2, e)
r3e <- crop(r3, e)
Inspect the resulting extents
t(sapply(c(r1e, r2e, r3e), function(i) as.vector(extent(i))))
# [,1] [,2] [,3] [,4]
#[1,] 533547.0 646547.0 7881284 7973784
#[2,] 533550.8 646550.8 7881307 7973807
#[3,] 533539.0 646539.0 7881305 7973805
They are not exactly the same, because that is not possible because the rasters do not align. Their "origins" are different
t(sapply(c(r1e, r2e, r3e), origin))
# [,1] [,2]
#[1,] 47.0 -16
#[2,] -49.2 7
#[3,] 39.0 5
To make them align, you would need to do something like this
r1e <- crop(r1, e)
r2e <- resample(r2, r1e)
r3e <- resample(r3, r1e)
t(sapply(c(r1e, r2e, r3e), function(i) as.vector(extent(i))))
# [,1] [,2] [,3] [,4]
#[1,] 533547 646547 7881284 7973784
#[2,] 533547 646547 7881284 7973784
#[3,] 533547 646547 7881284 7973784
Upvotes: 1