Mosaic rasterstacks based on minimum of certain layer?

Question

I am trying to mosaic 42 remote sensing rasterstacks (with 250 bands) based on the criterion that in overlapping areas, the pixel should be taken that has the most nadir viewing angle

Beside my rasterstacks I also have 42 rasters (so one for each stack) with the corresponding viewing angle for each pixel.

Any idea how to solve this?

I tried to include the viewing angle raster in the stack, and use something similar to

mosaic(a,b,fun=function(x)(min(x[[251]])) 

but that didn't work...

Any advice?

Answer 1

This is not a full solution, but it is too big for a comment.

The main problem that you are facing is that the mosaic function appears to operates on a per-band basis. As such, you can't use the last band to determine how it should mosaic the first band.
A dirty work-around would be to by-pass the mosaic function entirely:

1. Determine the combined total extent of the resulting mosaic - use extent() and @xmin, @xmax, @ymin and @ymax for this
2. Create temporary versions of all your rasters in the full mosaic extent using the extend() function. Note the difference between extent() and extend().
3. Create a stack of all the full-sized viewing angle rasters.
4. Determine which viewing angle is the best for each pixel using which.min() on view angle raster stack.

5. Use the new 'best angle' raster for determining which of your original rasters to use data from by utilizing a double for-loop - see code example below:

   #First we create an empty raster
   TemplateRaster <- setValues(raster(FullExtentOfMosaic,nrows=nrow(ViewAngleMosaicStack),ncols=ViewAngleMosaicStack)),rep.int(0,ncell(ViewAngleMosaicStack)))
   #Then we loop away over first the number of bands and then the number of files:
   for(i in 1:250){
   IterationRaster <- TemplateRaster
   for(j in 1:42){
   Driver <- BestViewAngleRaster==j
   IterationRaster <- IterationRaster+raster(ListOfExtendedHyperspectralRaster[j],band=i)*Driver
   }
   outputname <- paste("Hyperspectral_Band_",i,".tif",sep="")
   writeRaster(IterationRaster,filename=outputname,format="GTiff",overwrite=T)
   }

Basic idea of the code in step 5 is to write 250 single band rasters each containing the value as determined by the raster which was made using the RasterStack of view angles and which.min(). These 250 single band rasters can then be stacked and written into one big hyperspectral raster.

It should be noted that I haven't fully tested the above code on a 250-band raster, so I am uncertain on the required processing time and I have no doubt that there are more efficient methods available.