Extract MODIS phenology dates (as DOY) from the MCD12Q2 HDF files and export as tiff.

extract_MODIS_phenology_layer.r

# calculate adjusted phenology dates
# for greenness onset and maximum

extract_MODIS_phenology_layer <- function(value="increase",tiles_file="mytiles.txt"){
  
  # load required libraries
  require(raster)
  require(MODIS) # to import the SDS layers
  
  # list all hdf files, and extract unique years
  all_hdf_files = list.files(".",pattern=paste("^.*\\.hdf$",sep=""),recursive=TRUE)
  years = unique(substr(sapply(strsplit(all_hdf_files,split="MCD12Q2.A"),"[[",2),1,4))
  
  # correction factor needed to convert to DOY
  correction_factor = (as.numeric(years) - 2000)*365
  
  # list files and tiles
  #all_files = strsplit(list.files(".",pattern="*.hdf",recursive=TRUE),split='\\.')
  tiles = unlist(read.table(tiles_file))
  
  # for a given tile, for every year extract the layer of interest
  # convert to DOY, remove 32767 values, write to disk
  
  # for all tiles calculate
  for (j in tiles){
    
    # print some feedback
    print(j)
    
    # list all files for a given tile
    hdf_files = list.files(".",pattern=paste("^.*\\.",j,"\\..*\\.hdf$",sep=""),recursive=TRUE)
    
    # extract the SDS layer from the file headers
    if(value == "increase"){
      SDSlayers = unlist(lapply(hdf_files,function(x)getSds(x)$SDS4gdal[1]))
    }
    
    if (value == "maximum"){
      SDSlayers = unlist(lapply(hdf_files,function(x)getSds(x)$SDS4gdal[2]))
    }
    
    if (value == "decrease"){
      SDSlayers = unlist(lapply(hdf_files,function(x)getSds(x)$SDS4gdal[3]))
    }
    
    if (value == "minimum"){
      SDSlayers = unlist(lapply(hdf_files,function(x)getSds(x)$SDS4gdal[4]))
    }
    
    if (value == "amplitude"){
      SDSlayers_min = unlist(lapply(hdf_files,function(x)getSds(x)$SDS4gdal[5]))
      SDSlayers_max = unlist(lapply(hdf_files,function(x)getSds(x)$SDS4gdal[6]))
    }
    
    if (value == "area"){
      SDSlayers = unlist(lapply(hdf_files,function(x)getSds(x)$SDS4gdal[7]))
    }
    
    # check if the number of years equals the number of SDSlayers
    # if not skip to the next tile (need to download missing data)
    if (length(years) != length(SDSlayers)){
      print("skipping, missing years")
      next
    }
    
    # one by one open the SDSlayer - apply the correction factor
    # needed by the MODIS phenology product and save as a new file
    # [compressed geotiff]
    
    # process amplitude values
    if (value == "amplitude"){
      for (i in 1:length(SDSlayers)){
        print(SDSlayers_min[i])
        r_min = raster(SDSlayers_min[i],b=1)  
        r_max = raster(SDSlayers_max[i],b=1)  
        r = (r_max - r_min) * 0.0001 # this is the multiplier for EVI
        if (i == 1){  
          assign("c",r,envir=.GlobalEnv) # create composite
        } else {
          c = addLayer(c,r) # add to composite
        }
      }
      
      # write to file
      filename = paste("./",value,".",j,".tif",sep="")
      writeRaster(c,filename,overwrite=TRUE,options=c("COMPRESS=DEFLATE"))
      
      # clean up
      removeTmpFiles()
      rm(c) # remove stacked file
    }
    
    # process area under the curve values
    if (value == "area"){
      for (i in 1:length(SDSlayers)){
        print(SDSlayers[i])
        r = raster(SDSlayers[i],b=1)  
        r = r * 0.01 # this is the multiplier for EVI
        if (i == 1){  
          assign("c",r,envir=.GlobalEnv) # create composite
        } else {
          c = addLayer(c,r) # add to composite
        }
      }
      
      # write to file
      filename = paste("./",value,".",j,".tif",sep="")
      writeRaster(c,filename,overwrite=TRUE,options=c("COMPRESS=DEFLATE"))
      
      # clean up
      removeTmpFiles()
      rm(c) # remove stacked file  
      
    } 
    
    # process phenology values
    if (value == "increase" || value == "maximum" || value == "decrease" || value == "minimum"){
      for (i in 1:length(SDSlayers)){
        r = raster(SDSlayers[i],b=1)
        r = r - correction_factor[i]
        if (i == 1){
          assign("c",r,envir=.GlobalEnv)
        } else {
          c = addLayer(c,r)
        }
      }
      
      # write to file
      filename = paste("./",value,".",j,".tif",sep="")
      writeRaster(c,filename,overwrite=TRUE,options=c("COMPRESS=DEFLATE"))
        
      # clean up temporary scratch files
      # remove c raster stack (composite of all years of a given tile)
      removeTmpFiles()
      rm(c) # remove stacked file
    }
    
  } # end looping over tiles

}