0

I have 12 GeoTIFFs to stack in R and have run in to a problem. When I try to process the images as a stack in R it fails as their extent differs.

The rasters originate from Google Earth Engine. They are extracts of processed Landsat 8 scenes.

The AOI is on the boundary between two Landsat paths. Due to cloud/scene availability, some images derive from one path, and some the other.

To understand the problem (and as I am relatively new to R) I mapped the images in QGIS. Although they are aligned some have the CRS EPSG:32653, and some EPSG:32654 (corresponding to the different paths) which tilts one vs the other:

Mis-matched image projections

The [complete] code in R is works at all other study sites, and works with another set of rasters at this site, but fails with the following error:

# Bungala

# Set directory for geoTIFFs
pre_tif_Dir <- "..."
post_tif_Dir <- "..." 

# Create filename list
pre_fileNames <- list.files(pre_tif_Dir,full.names=TRUE,pattern="tif$")
post_fileNames <- list.files(post_tif_Dir,full.names=TRUE,pattern="tif$")

# Import mask shapefile (masked buffers)
mask <- shapefile(".../Bungala_maskedBuffers.shp")

## Create a stack (files don't yet have same crs so can't use brick)
pre_stk <- stack(as.list(pre_fileNames))
post_stk <- stack(as.list(post_fileNames))

pre_stk generates the error:

Error in compareRaster(rasters) : different extent

My code notes from late last night are then:

# Extent differs for some scenes (months 1,3,5,0,12 are slightly rotated due to differing crs)
# Tried to crop/correct in QGIS.  Didn't work.

# raster <- projectRaster(raster, crs=crs(mask)) ?  Do each image individually?

I have tried reprojecting and cropping, but this hasn't worked - either the image shifts many miles away, or the tilt remains as the CRS issue is unchanged. I don't know how to do this in R, so had to do it in QGIS this morning.

How do I correct the CRS to maintain alignment, but create a stack for processing using R?

2
  • 2
    reprojecting shouldn't shift things wrongly. Of the two you've shown, projecting one to 32653 should work. can you do that and show us the extents you get before and after? or can you point us to the data so we can try this?
    – Spacedman
    Mar 19, 2022 at 11:33
  • Sure - the data are at livelancsac-my.sharepoint.com/:f:/g/personal/… Mar 19, 2022 at 14:31

1 Answer 1

1

This works for me. First get all file names:

tifs <- list.files(".","*.tif")

then get the first one as a raster, and its coordinate system. We're going to transform the others to this projection, and crop/extend to this raster's extent:

r1 <- raster(tifs[1])
cr <- projection(r1)

Next loop over the rest, and if the projection isn't what we want, use projectRaster using to to project it to the grid and CRS of the first raster. Then use crop and extend to cut the area of interest. I'm not sure if both of these are essential but...

rlist <- lapply(tifs[-1],function(rfile){
    r = raster(rfile)
    if(projection(r) != cr){
        r = projectRaster(r, to = r1)
    }
    r = crop(extend(r, r1), r1)
    
    return(r)
}
)

Then put the first one back and make a stack:

s <- stack(c(r1, rlist))

enter image description here

I've not got the shapefile you might be masking to but with crop and extend this should all work. I've not seen any problems with raster transforms.

2
  • Thanks for this. I understand the principle but my lack of R knowledge (especially syntax) really slows me down. :) Mar 19, 2022 at 16:55
  • I had to set the working directory to the one containing the rasters, but the mask seems to work fine (as you suggested) and with the stack working I can complete the analysis I needed to. Thanks for the help. Mar 23, 2022 at 20:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.