# loading the bands
B2 <- raster("./Data/LC08_L1TP_015053_20160413_20170223_01_T1_B2.tif")
B3 <- raster("./Data/LC08_L1TP_015053_20160413_20170223_01_T1_B3.tif")
B4 <- raster("./Data/LC08_L1TP_015053_20160413_20170223_01_T1_B4.tif")
B8 <- raster("./Data/LC08_L1TP_015053_20160413_20170223_01_T1_B8.tif")

# then I want to clip the rasters to the study site
# first I have to load the raster of the Costa Rica provinces
shape <- readOGR("./Data/gadm36_CRI_1.shp")

# Then add that projection to the shapefile
shape.proj <- spTransform(shape, CRS(proj4string(B2)))
proj4string(shape.proj)

# Combining the two provinces that I want : Heredia and Alajuela
p1 <- shape.proj[shape.proj$NAME_1=='Heredia', ]
p2 <- shape.proj[shape.proj$NAME_1=='Alajuela', ]
p2$NAME_1 <- NULL
x <- bind(p1, p2)

# Crop raster by extent of shapefile
studysiteB2 <- mask(crop(B2, x), x)
studysiteB3 <- mask(crop(B3, x), x)
studysiteB4 <- mask(crop(B4, x), x)
studysiteB8 <- mask(crop(B8, x), x)

# now stack the Bands together
image <- stack(studysiteB2,studysiteB3,studysiteB4,studysiteB8)

Error in compareRaster(x) : different extent

extent(studysiteB2)

class       : Extent
xmin        : 710385
xmax        : 860325
ymin        : 1087635
ymax        : 1225725

extent(studysiteB8)

class       : Extent
xmin        : 710392.5
xmax        : 860332.5
ymin        : 1087628
ymax        : 1225733