Why does the Extent of NIR Band 8 Landsat image not match other bands, even when cropping to a specified extent?

Question

I'm trying to load and then stack four bands 2,3,4 and 8 RGB and NIR from Landsat images. I have downloaded and extracted them and they are in .tif format.

# 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)

But then I am getting this error message

Error in compareRaster(x) : different extent

When I look at the extent of the cropped bands 2,3 and 4 all have the same extent: e.g.

extent(studysiteB2)

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

But the extent of B8 cropped image is .5 more on the xmin / xmax

extent(studysiteB8)

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

This is the same on the original images as well. Any ideas on why this is happening and how to solve it?

Here is a DropBox link to the images: https://www.dropbox.com/sh/fk2acf07lc9437l/AAD-DCopQxFFIwHRWAfew2PMa?dl=0

Answer 1

A shell-script one-liner tells me the basic geometry of the TIFs:

$ for f in *TIF ; do  echo $f; gdalinfo $f | grep Size; done

Which outputs this:

LC08_L1TP_015053_20160413_20170223_01_T1_B2.TIF
Size is 7551, 7711
Pixel Size = (30.000000000000000,-30.000000000000000)
LC08_L1TP_015053_20160413_20170223_01_T1_B3.TIF
Size is 7551, 7711
Pixel Size = (30.000000000000000,-30.000000000000000)
LC08_L1TP_015053_20160413_20170223_01_T1_B4.TIF
Size is 7551, 7711
Pixel Size = (30.000000000000000,-30.000000000000000)
LC08_L1TP_015053_20160413_20170223_01_T1_B8.TIF
Size is 15101, 15421
Pixel Size = (15.000000000000000,-15.000000000000000)

Which tells me the third TIF (B8) is a higher resolution image, and so it looks like your mask is getting a slightly larger extent when cropped because of the finer pixels.

Resample B8 to the same resolution as the other bands first?

Answer 2

I don't know how to solve the issue directly but I found a different way round it, following instructions in "Remote Sensing and GIS for Ecologists" by Wegman, Leutner and Dech. This also removed the need to load the large images. Before I had loaded the bands separately rather than via the meta data and this, for some reason, produced the small differences.

# Using the meta data obtained from the Landsat download:
meta <- readMeta("./Data/LC08_L1TP_015053_20160413_20170223_01_T1_MTL.txt")
summary(meta)
# Then using the stackMeta function to Stack all the bands together
site2016 <- stackMeta(meta)
site2016

I changed the extent to be a simple box rather than a .shp file, but it works either way.

# now to clip this to extent I want:
ex <- extent(c(785313, 841487, 1145281, 1199888))
site_crop <- crop(site2016, ex)
plot(site_crop)