Seperate 8-Bit MODIS (MOD10A2) snow cover data in R (MultiCore)

with No Comments

Description

I recently had to seperate 8-BIT MODIS data into individual days. The MODIS snow cover product consists of one file for eight days. The files are therefore in a 8-Bit format where every Bit represents one day and where 0=nosnow and 1=snow. Have a look at the MODIS snow cover description. First I wanted to have a look at the data, so I wanted to batch process all my files (690 8Bit files) to output an image of the snow extent for every single day. I then wanted to save the one raster file for every day to be able to easily further process the data. To achieve this I wrote a R script, which utilises various additional libraries. Due to the fact that this takes a ton of processing time (>10h with 1 x 3.00GHz) for so many files when using only 1 core (R standard), I have written the script to utilise as many cores as possible (<1h with 32 x 3.00GHz). Feel free to comment, or add functionality. I have not had time to thoroughly test this script, but it seems to be working for 8-Bit MODIS snow extent data.

The Code

#================================================================================================================#
#           .--.              .--.                                                                               #
#          : (\ ". _......_ ." /) :       Author: Manuel Bär @ GeoNET                                            #
#           '.    `        `    .'        Date: 09.04.2015                                                       #
#            /'   _        _   `\         Version: 1.1                                                           #
#           /     0}      {0     \        Features: MutliCore Enabled!                                           #
#          |       /      \       |                                                                              #
#          |     /'        `\     |       Description: Iterates through a folder containing 8-Bit MODIS snow     #
#           \   | .  .==.  . |   /        cover files. Takes each file, splits the file into single 8 days and   #
#            '._ \.' \__/ './ _.'         according to options saves daily TIF-files, daily images of extent or  #
#            /  ``'._-''-_.'``  \         both. Also cleans up overlapping days.                                 #
#                    `--`   jgs                                                                                  #
#================================================================================================================#


    #install needed packages and libraries
    install.packages("fields")
    install.packages("raster")
    install.packages("rgdal")     

    library(raster)
    library(rgdal)
    library(fields)
    library(parallel)
        
    # ######################:-  OPTIONS  :-######################
    # Define various options
    
    # --------- Directories and files -------------
    # Define directory where original data is stored 
    dataDir <- "c:/test_Data/"
    # Define the filename pattern using REGEX
    fileName <- "[a-zA-Z0-9_]_lyr2$"
    # Define output directory of files
    outDirFiles <- "c:/dailySnowRasterTest"
    # Define output directory of images
    outDirImages <- "c:/dailySnowImagesTest" 

    # --------- Additional options -------------
    # Specify number of Bits (only tested with 8BIT MODIS files where each bit is a day
    noBits <- 8
    # List contaiing days of each year of Data being processed
    daysInYears <- list(
           "2000" = 366 #2000-2001
          ,"2001" = 365 #2001-2002
          ,"2002" = 365 #2002-2003
          ,"2003" = 365 #2003-2004
          ,"2004" = 366 #2004-2005
          ,"2005" = 365 #2005-2006
          ,"2006" = 365 #2006-2007
          ,"2007" = 365 #2007-2008
          ,"2008" = 366 #2008-2009
          ,"2009" = 365 #2009-2010
          ,"2010" = 365 #2010-2011
          ,"2011" = 365 #2011-2012
          ,"2012" = 366 #2012-2013
          ,"2013" = 365 #2013-2014
          ,"2014" = 365 #2014-2015
    
    )
    
    # Should daily files be created? (0=no;1=yes)
    createFiles <- 1
    # Should daily images be created? (0=no;1=yes)
    createImages <- 1
    
        
    # #########################################
    # Load Files into array
    
    fileList <- list.files(path = dataDir, pattern = fileName)
    
    
    # #########################################
    # Function to convert INT to BINARY with optional input of number of BITS
    
    number2binary = function(number, noBits) {
        binary_vector = rev(as.numeric(intToBits(number)))
        if(missing(noBits)) {
            return(binary_vector)
        } else {
            binary_vector[-(1:(length(binary_vector) - noBits))]
        }
    }
        
    # #########################################
    # THE MAGIC IS REAL!
    # Function that opens a file, converts 8day-files to 
    # daily files and saves in format

    raster2files <- function(i) {
        # Load files into raster
        r <- readGDAL(paste(dataDir,fileList[i],sep=""))
        rRas <- raster(r)
        nCols = ncol(rRas)
        nRows = nrow(rRas)
        rTemp <- raster(ncols=nCols, nrows=nRows)
        
        # #########################################
        # Save PNG
        
        for(j in 1:noBits){  
            
            fun <- function(x){number2binary(x,noBits)[(noBits+1)-j]} 
            rTemp <- calc(rRas,fun)
            fileYear <- substr(fileList[i],9,12)
			daysInYear <- as.numeric(daysInYears[toString(fileYear)][1])
            fileNumber <- as.numeric(substr(fileList[i],13,15))+(j-1)
			if(fileNumber <= daysInYear){			
				if(createFiles==1){
					fileName <- paste(fileYear,"_",fileNumber,".tif",sep="")
					setwd(outDirFiles)
					writeRaster(rTemp, filename=fileName, overwrite=TRUE)
				}
				if(createImages==1){
					fileName <- paste(fileYear,"_",fileNumber,".png",sep="")
					setwd(outDirImages)
					png(fileName)
					image.plot((rTemp),zlim=c(0,1),col=bpy.colors(2), main=paste("Year: ",fileYear," ; Day: ",fileNumber, sep=""), axis.args=list(at=c(0,1), labels=c(0,1)))    
					dev.off()
				}
			}
        }        
    }
    
    # #########################################
    # Setup cluster for multicore processing
    
    numCores <- detectCores()            
    cl <- makePSOCKcluster(numCores) 
    setDefaultCluster(cl)
    
    # #########################################
    # Give the cluster workers the needed libraries

    clusterEvalQ(NULL, library(rgdal))
    clusterEvalQ(NULL, library(raster))
    clusterEvalQ(NULL, library(fields))
    clusterExport(NULL, c('number2binary','raster2files','dataDir','outDirFiles','outDirImages','createFiles','createImages','daysInYears','fileList','noBits'))
    
    # #########################################
    # Distibute processing to the different workers 
    inputs <- 1:length(fileList)
    parLapply(cl, inputs, raster2files)  
    
    # #########################################
    # Close clusters again
    stopCluster(cl)  
    
    
#================================================================================================================#
# Improvements and ToDo: 
# - More options to specify number of cores?
# - Cluster in network?
# 
# --------------------------------------------------------------------------------
#
# - If you have any additional comments or would like to contribute to the code, 
#   please contact me: mbaer[at]geonet[dot]ch
#  
#================================================================================================================#