[R] R getting slower until it breaks...

[Bastien Ferland-Raymond]

Hello R-users,

I'm currently facing a pretty hard problem which I'm hopping you'll be able to help me with.  I'm using R to create images.  That alone is not the problem, the problem is that I'm using R to create 168 000 images...  My code (which is given below) use different package (raster and rgdal) to import a image (size 20gig) and divide it into 168 000 pictures that are 100 pixel x 100 pixel.  The code works fine for making the images, but if I ask it to run all 168 000, it always breaks around 15 000.  

It starts with the code being able to make around 2 pictures per second, but then it slows down and after around 2000 pictures it's only 1 picture per second.  Later on it's getting closer to 1 pictures every 3 seconds etc.  until it bugs.  I have no error message, only Windows that tells me that "R encounter a problem and most be close..."  Initially I though it was a Windows problem, that I couldn't put too many file into a folder and it was slowering it down.  Then I divided my batch process into smaller (5000 files) folder but it didn't help, still breaks at 15 000.  I also try to do gc() after each 5000 pictures to save memory but it didn't help either.  I removed every loops from the code because I thought it was the problem, but it was just faster at bugging... After the bug, I need to restart the computer if I want to go back to the initial speed.

I'm pretty much running out of options.  It's there limitation in R as the number of files it can create in one session?  Is it a windows problem?  Is there better way to clear the memory than gc()? Any thought on that?

I'm using R 2.11.1, win XP, my hard drive is NTSF, computer: intel core2 duo E6750 32 bit with 2 gig of Ram.

Here is my code, but I doubt it would help much with my problem:

########
# It made of 4 functions (sorry, it's french):

##########################################################################
##########################################################################
###  Ensemble des fonctions pour faire les images NDVI rouge et verte  ###
##########################################################################
######  Bastien Ferland-Raymond, 5 oct 2010  #############################
##########################################################################

########
## Simplement rouler le script au complet
########
### Library nécessaire:
library(raster)
library(rgdal)
library(shapefiles)

#############################################################################
## Fonction 1  -  NDVI a partir de coordonnee Pixel et largeur #####  
 calculate_NDVI<- function(Type, object, VALUE) { 
   redorgreen <- ifelse(Type=="red",2,3)
   list1 <- unstack(object)
   rast1 <- list1[[1]]
   rast2 <- list1[[redorgreen]]
   NAvalue(rast1)<- -99999
   NAvalue(rast2)<- -99999
   cells1 <- getValuesBlock(rast1,row=VALUE[[2]],nrow=VALUE[[3]],col=VALUE[[1]],ncol=VALUE[[3]])
   cells2 <- getValuesBlock(rast2,row=VALUE[[2]],nrow=VALUE[[3]],col=VALUE[[1]],ncol=VALUE[[3]])
   cells1[is.na(cells1)]<-0;
   cells2[is.na(cells2)]<-0;
   calculNDVI <-(cells1 - cells2) / (cells1 + cells2)
   NDVImatrix <- matrix(calculNDVI,nrow=VALUE[[3]],ncol=VALUE[[3]], byrow=TRUE)
   NDVImatrix <- NDVImatrix + 1
   NDVImatrix <- NDVImatrix * (255/2)
   return(NDVImatrix)
   }
#################################################################################
## Fonction 1b  -  Faire le tiff
 make.tiff<- function(NV=newValues,TT=Type,img=imgRaster,nom){
 pixelNDVIMatrix <- calculate_NDVI(TT,img,NV[c(1,2,3)])
 newRaster <- raster(pixelNDVIMatrix)
 NAvalue(newRaster)<-999999
 nnom<-nom[NV[4]]
 writeRaster(newRaster, filename=nnom,datatype="INT1U",format="GTiff",overwrite=FALSE)
 aaa<-2
}
  
#################################################################################
## Fonction 2  -  Creation de fonction convertissant les coordonnee metrique en coordonnee pixels #####  
 latlong_to_pixels<- function(Coord, facteur, meterWidth=NULL) {   #Coord doit être c(x,y)
  newX <- Coord[1] / facteur
  newY <- Coord[2] / facteur
  if(!is.null(meterWidth)){
   newWidth <- meterWidth / facteur
   return(c(newX,newY,newWidth))
  }
  return(c(newX,newY))
 }

#############################################################################
####  Fonction 3  -  Fonction principal   #####
 make.NDVI.photo<- function(tableDesPlacettes,Type,newImagesDirectory,textAndImgDirectory="U:\\kNN_Valcartier\\Photo aerienne"){
 lastWD<- getwd()
setwd(textAndImgDirectory)
 imgRaster<- stack(imageAssociee)
 x1<- tableDesPlacettes[,2] - xmin(imgRaster) - (tailleFenetres/2)     # The image origin for calculation is in the top left corner
 y1<- ymax(imgRaster) - tableDesPlacettes[,3] - (tailleFenetres/2)
coo <- cbind(x1,y1)
 newValues<- t(apply(coo,1,latlong_to_pixels,facteurMetreParPixel,tailleFenetres))
 newImgName<- paste(newImagesDirectory,substr(Type,1,3),"_","GC",tableDesPlacettes[,1],".tif",sep="")
apply(cbind(newValues,1:length(newImgName)),1,make.tiff,Type,imgRaster,nom=newImgName)
setwd(lastWD)
}


###########################
## Executing fonctions:
#############################


###  loader les données brutes de fenêtres
fichier.fenetre.brute<-read.dbf("U:\\kNN_Valcartier\\Fenêtres 30x30 1 octobre\\168700_30m_centroid.dbf", header=T) 
###  Sélectionner les fenêtres complètes
fenetre.complete<-round(fichier.fenetre.brute[[1]][,2],1)==900
###  Sortir les centroides pour extraction
centro.tout.900<-fichier.fenetre.brute[[1]][fenetre.complete,c(1,5,6)]
#rm(fichier.fenetre.brute) ; gc()

## données nécessaires pour la fonction
 imageAssociee<- "mosaique_all_v1.img"  # nom de l'image
 facteurMetreParPixel<- 0.3         # combien de metre vaut un pixel
 tailleFenetres= 30            # en metre

 start.time<-Sys.time();start.time
make.NDVI.photo(centro.tout.900[19137:24136,],"red","BFR\\NDVI_red_fenetre\\batch 3\\")
 stop.time<-Sys.time()
 time.run<-stop.time-start.time
 alarm()
 time.run 
gc()
 start.time<-Sys.time();start.time
make.NDVI.photo(centro.tout.900[24137:29136,],"red","BFR\\NDVI_red_fenetre\\batch 4\\")
 stop.time<-Sys.time()
 time.run<-stop.time-start.time
 alarm()
 time.run 
gc()

#############
    

Voilà, 

Thanks!

Bastien Ferland-Raymond