5

I am writing a script for the processing of different hyperspectral datacubes. The datacubes contain flight lines that I've recorded with a UAV. In total there are 10 flight lines that need to be processed and mosaiced. The files are in .bsq format and contain 101 bands:

    class       : RasterBrick 
dimensions  : 1293, 749, 968457, 101  (nrow, ncol, ncell, nlayers)
resolution  : 0.25, 0.25  (x, y)
extent      : 651041.9, 651229.1, 5687300, 5687623  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=31 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
data source : \HyperspectralDatacube_FlightLine02_PixelSize=0.25m_DSM=Photogrammetric.bsq 
names       : X450.00000.Nanometers, X455.00000.Nanometers, X460.00000.Nanometers, X465.00000.Nanometers, X470.00000.Nanometers, ... 

I load all 10 hyperspectral rasters using lapply and brick to load all 101 layers:

# Load hyperspectral datacubes
files <- list.files(path=".", pattern = ".bsq$")
r <- lapply(files, brick)

Then I want to calculate vegetation indices for these hyperspectral rasters and mosaic the results.

# Create NDVI function
NDVI <- function(R,NIR) 
{
  (NIR-R)/(NIR+R)
}

# Calculate NDVI
ndvi1 <- NDVI(raster(r[[1]], layer=41),(raster(r[[1]], layer=61)))
ndvi2 <- NDVI(raster(r[[2]], layer=41),(raster(r[[2]], layer=61)))
.......

# Mosaic rasters
mndvi <- mosaic(ndvi1,ndvi2, .... fun=mean, tolerance=0.05)

How can I create a loop for the NDVI calculations? Instead of writing the code line by line like the example above. I'm fairly new to R and I dont quite understand yet how I can write this loop with the variable rasters and layers for each calculation. The results need to be mosaiced like the example above.

2

For the raster generation just look up loop constructs in R.

The problematic bit is the last line, since you want to avoid specifying all the rasters manually. You can add them to a list in the same loop and use that in the mosaic call. Look at this answer to see how do.call can be used if your R version doesn't support a list parameter to mosaic anymore. It also has loop examples.

  • Thank you for your help! The do.call function seems very helpful, as indeed i could not mosaic lists. I still however cannot get the loop to work. I've now tried this for looping the NDVI function on the raster layers: ndvi <- lapply(r, (NDVI(raster(r[layer=41]),(raster(r[layer=61]))))) I do however get this error: Error in .local(x, ...) : list has no "x" – cf2 Nov 2 '15 at 11:26
  • Almost a new question, but as always, check if the subexpressions work first and then go higher. My R is rusty, but isn't r[layer=41] a whole column? – lynxlynxlynx Nov 2 '15 at 18:14

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.