0

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

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?

  • thanks for the answer @Spacedman, it made me go back and start again from the beginning. Solution below. Cheers – TomCLewis Nov 21 '18 at 8:45
0

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)
  • Where do the readMeta and stackMeta functions come from? – Spacedman Nov 21 '18 at 9:07
  • They are from the RStoolbox package – TomCLewis Nov 21 '18 at 11:42

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.