1

I need to apply r.neighbours or Neighbourhood analysis (QGIS plugin LecoS (or equivalent) to 3000 rasters and define a neighbourhood of 5 x 5 and calculate the mean.

I suppose I need to create a script with a loop that goes through the raster files and applies the same analysis settings. Unfortunately I only know how to program in R. So I'm asking if anyone can show me how to create a script that does this.

I use QGIS and I would prefer to do this using the QGIS Processing Framework. But GRASS GIS or SAGA GIS are also OK for me.

Thanks in advance


EDIT: This how I did it on Ubuntu 12.04. Thanks @Curlew.

I created GRASS Location/Mapstet and batch imported the 3000 rasters using "Common formats import (r.in.gdal)" to the mapset (I recommend deactivating the option "add to layer tree"). I used one of the rasters for setting the region (g.region).

Then I opened R, set the workspace to the folder where the rasters are stored and run this:

library(raster)
library(spgrass6)
initGRASS(gisBase = "/usr/lib/grass64", home = tempdir(), gisDbase = "/path/to/gis/data/directory", location = "Location", mapset = "Mapset", override = TRUE) #Initate the previously created Location/Mapset

ls <- lapply(list.files(".","*.tif"),FUN=raster) # Import all raster to the list "ls"

for(i in seq_along(ls)){
execGRASS("r.neighbors", input=names(ls[[i]]),size=5,output="teste",flags=c("overwrite"))
execGRASS("r.out.gdal", input="teste",output=paste(paste("/path/to/output/folder",sep="",names(ls[[i]])),sep="",".tif"))
}

This solves my problem, but I'm still curious about how the equivalent Python code for QGIS Processing looks like.

1 Answer 1

2

Just quick, as i don't have the time to code sth. for you. You know that you can call GRASS and also SAGA modules in R?
For GRASS: spgrass6
For SAGA: RSAGA

Then you just read in all the files you need and build a loop to calculate the needed statistics.

library(raster)
library(RSAGA)
library(spgrass6)
ls <- lapply(list.files(".","*.tif"),FUN=raster)

for(i in ls){
  res <- doGRASS("r.neighbors") # like that
  res <- RSAGA::rsaga.filter.simple() # or that
  # you could write the result to a new file or use assign to save it as variable
}

However i would suggest to learn a bit python.

3
  • Hi. Thanks. One question: I see that I can do stuff like ths: "execGRASS("r.neighbors", input=example,size=5,output="exp",flags=c("overwrite"))" to GRASS objects in the mapset. So, what is "ls <- lapply(list.files(".","*.tif"),FUN=raster)" for? Commented Jan 31, 2014 at 12:17
  • ahh, it was generally to import all *.tif files in the current folder into R. Figured you might need sth. like that with 3000 raster layers. You can modify the command to load in the files as srgd for SAGA or into GRASS instead as well...
    – Curlew
    Commented Jan 31, 2014 at 12:51
  • I see. It is handy. Is it possible to convert a "RasterLayer" to GRASS Raster and upload it to the Mapset? writeRAST6(), according to the help file, only accepts "SpatialGridDataFrame". Commented Jan 31, 2014 at 13:40

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.