[GIS] Crop raster based on another raster R

Question

library(raster)

How can one clip one raster based on another raster?

RasA:

class       : RasterLayer 
dimensions  : 459, 533, 244647  (nrow, ncol, ncell)
resolution  : 10000, 10000  (x, y)
extent      : 3685000, 9015000, 655000, 5245000  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=lcc +lat_1=49 +lat_2=77 +lat_0=63.390675 +lon_0=-91.86666666666666 +x_0=6200000 +y_0=3000000 +ellps=GRS80 +datum=NAD83 +units=m +no_defs +towgs84=0,0,0 
data source : in memory
names       : Total.Precipitation 
values      : 0, 11.85  (min, max)

The resolution of RasA is 10km*10km

sample data for RasA an be found here. Shapefile included

RasB:

class       : RasterBrick 
dimensions  : 416, 885, 368160, 2  (nrow, ncol, ncell, nlayers)
resolution  : 0.1, 0.1  (x, y)
extent      : -141.1, -52.6, 41.6, 83.2  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0 
data source : in memory
names       : layer.1, layer.2 
min values  :       0,       0 
max values  :    22.6,    27.8 

The resolution of RasB is 0.1deg*0.1deg

Question:

1) crop RasA to the extent of RasB

Using @mdsummer's suggestion I solved a previous projection issue using the following method:

prr=projectRaster(RasA,RasB,method = "ngb")
range(values(prr) , na.rm = T )

But trying to change the resolution of RasA during (prr=projectRaster(RasA,RasB,method = "ngb",res=10000) and after ( res(ppr)=10000) reprojection appears not to have any effect (see below). pprcontains the same dimensions and values as RasA but the resolution is copied from RasB.

class       : RasterLayer 
dimensions  : 416, 885, 368160  (nrow, ncol, ncell)
resolution  : 0.1, 0.1  (x, y)
extent      : -141.1, -52.6, 41.6, 83.2  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0 
data source : in memory
names       : Total.Precipitation 
values      : 0, 11.85  (min, max)

Answer 1

Doing this through projection is quite weird, as you are trying to match extensions with rasters of different resolutions.

Here I suggest two options: (1) resample the clipping raster to the desired final resolution first and then projectRaster with that layer, or (2) just crop the raster to keep the original resolution (10000x10000 in your case).

Here is an example with downloadable layers:

library(raster)
library(rgdal)
library(ggplot2)
library(mapview)
library(gdalUtils)

# load the data
r1 <- raster('./data/MDT25.tif')
r2 <- raster('./data/LIDAR.tif')

# change projection if you have to (in my case r2 crs)
r2p <- projectRaster(r2, crs=crs(r1))

# resample to the clipping layer and mask r2 with it 
r3 <- crop(r1,r2p)

# see original
r1; r2; r3; mapview(r1)+r2+r3

Results:

> r1; r2; r3; mapview(r1)+r2+r3
class      : RasterLayer 
dimensions : 609, 580, 353220  (nrow, ncol, ncell)
resolution : 25, 25  (x, y)
extent     : 636512.5, 651012.5, 4803362, 4818588  (xmin, xmax, ymin, ymax)
crs        : +proj=utm +zone=29 +ellps=GRS80 +units=m +no_defs 
source     : /media/cesarkero/D1/GDrive/Proyectos/Stackexchange/GIS Stack Exchange/20210127_Crop raster based on another raster R/data/MDT25.tif 
names      : MDT25 
values     : 11.324, 770.305  (min, max)

class      : RasterLayer 
band       : 1  (of  3  bands)
dimensions : 1768, 2713, 4796584  (nrow, ncol, ncell)
resolution : 6.306866e-05, 4.403295e-05  (x, y)
extent     : -7.310935, -7.13983, 43.40091, 43.47876  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
source     : /media/cesarkero/D1/GDrive/Proyectos/Stackexchange/GIS Stack Exchange/20210127_Crop raster based on another raster R/data/LIDAR.tif 
names      : LIDAR 
values     : 0, 255  (min, max)

class      : RasterLayer 
dimensions : 359, 564, 202476  (nrow, ncol, ncell)
resolution : 25, 25  (x, y)
extent     : 636562.5, 650662.5, 4806712, 4815688  (xmin, xmax, ymin, ymax)
crs        : +proj=utm +zone=29 +ellps=GRS80 +units=m +no_defs 
source     : memory
names      : MDT25 
values     : 13.023, 748.862  (min, max)