R/runMrt.R

In MODIS: Acquisition and Processing of MODIS Products

#' Run MODIS Reprojection Tool
#' 
#' @description 
#' Specifying input parameters, this function gets MODIS grid data from the 
#' archive (HTTP or local) and processes it with the MODIS Reprojection Tool 
#' (MRT). At any point, you are highly encouraged to consult the MRT User's 
#' Manual for further information.
#' 
#' @param product,collection,begin,end,extent,tileH,tileV,SDSstring,job See 
#'   [runGdal()] and functions linked therein.
#' @param datum The output datum used for datum conversion as `character`, 
#'   defaults to `"NODATUM"`. Supported datums are `"NAD27"`, `"NAD83"`, 
#'   `"WGS66"`, `"WGS72"` and `"WGS84"`, see MRT User's Manual, p. 7-8.
#' @param zone Output zone number as `integer`, relevant only for UTM 
#'   projections (i.e., `outProj = "UTM"`). Valid values are `-60` to `+60`.
#' @param projPara Output projection parameters as `character` string, see 
#'   Details. Ignored if 'outProj' is one of `c("SIN", "GEO")`. If not specified
#'   and using another target projection, the default settings for `"GEO"` are 
#'   assumed.
#' @param mosaic A `logical` that toggles mosaicking on (default) or off. One 
#'   example where `mosaic = FALSE` makes sense is for large spatial extents 
#'   because maximum supported HDF4 file size is 2GB. If crossed, mosaicking 
#'   will fail.
#' @param anonym A `logical`, defaults to `TRUE`. If `FALSE`, the job name is 
#'   appended to the root filename.
#' @param ... Additional arguments passed to [MODISoptions()], see also Details 
#'   for some MRT specific settings.
#' 
#' @return 
#' A `list` of output file names summarized by product and date, see also Value 
#' in [runGdal()].
#' 
#' @details 
#' Please note that in contrast to [runGdal()], MRT's `resample` function does 
#' not offer an 'overwrite' option, meaning that existing files will be 
#' overwritten (see also MRT User's Manual, p. 59).
#' Further arguments that require particular attention when operating MRT are 
#' summarized in the following list:
#'  
#' **dataFormat**:
#' Output file formats include:
#' 
#' * `"raw binary"`: `.hdr` and `.dat`
#' * `"HDF-EOS"`: `.hdf`
#' * `"GeoTiff"`: `.tif` (default)
#' 
#' Any other format specified through [MODISoptions()] or 'dataFormat' is 
#' ignored and set to `"GeoTiff"`.
#' 
#' **outProj**:
#' MRT uses calls to the General Cartographic Transformation Package (GCTP) and 
#' as such allows projection to the following mapping grids:
#' 
#' * Albers Equal Area (`"AEA"`)
#' * Equirectangular (`"ER"`)
#' * Geographic (`"GEO"`)
#' * Hammer (`"HAM"`)
#' * Integerized Sinusoidal (`"ISIN"`)
#' * Interrupted Goode Homolosine (`"IGH"`)
#' * Lambert Azimuthal (`"LA"`)
#' * Lambert Conformal Conic (`"LCC"`)
#' * Mercator (`"MERCAT"`)
#' * Molleweide (`"MOL"`)
#' * Polar Stereographic (`"PS"`)
#' * Sinusoidal (`"SIN"`)
#' * Transverse Mercator (`"TM"`)
#' * Universal Transverse Mercator (`"UTM"`)
#' 
#' See also 'References' and MRT User's Manual, pp. 6 and 29.
#' 
#' **projPara**:
#' Output projection parameters are autodetected for 
#' `outProj \%in\% c("SIN", "GEO")`:
#' 
#' * `"SIN"`: `"6371007.18 0.00 0.00 0.00 0.00 0.00 0.00 0.00 86400.00 0.00 0.00 0.00 0.00 0.00 0.00"`
#' * `"GEO"`: `"0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0"`
#' 
#' For detailed information on defining parameters for other target projections,
#' please refer to 'Appendix C: Projection Parameters' in the MRT User's Manual,
#' p. 65-66. 
#' 
#' @author 
#' Matteo Mattiuzzi, Forrest Stevens and Florian Detsch
#' 
#' @seealso 
#' [MODISoptions()], [runGdal()].
#' 
#' @source 
#' The MRT software has been [retired](https://lpdaac.usgs.gov/news/downloadable-modis-reprojection-tool-mrt-and-mrtswath-tool-have-been-retired), 
#' and is hence no longer officially available for download through LP DAAC.
#' 
#' @references 
#' Dwyer J, Schmidt G (2006) The MODIS Reprojection Tool, 162-177,
#' \doi{10.1007/978-3-540-37294-3_9}. In: Qu JJ, Gao W, Kafatos M, Murphy RE, 
#' Salomonson VV (eds) Earth Science Satellite Remote Sensing. Springer: Berlin, 
#' Heidelberg.
#' 
#' Elassal AA (1989) General Cartographic Transformation Package (GCTP), Version 
#' II. NOAA Technical Report NOS124 CGS9. NOAA: Rockville, MD, USA. Available 
#' online \href{https://www.cmascenter.org/ioapi/documentation/all_versions/html/GCTP.pdf}{here}
#' (2018-09-13).
#' 
#' @examples 
#' \dontrun{
#' geo = runMrt(product="MOD11A1", extent="austria", begin="2010001", end="2010002", SDSstring="101",
#'              job="ExampleGEOdelme", outProj="GEO")
#' sin = runMrt(product="MOD11A1", extent="austria", begin="2010001", end="2010002", SDSstring="101",
#'              job="ExampleSINdelme", outProj="SIN")
#' utm = runMrt(product="MOD11A1", extent="austria", begin="2010001", end="2010002", SDSstring="101",
#'              job="ExampleUTMdelme", outProj="UTM", zone = 33)
#' }
#' 
#' @export
runMrt = function(
  product
  , collection = NULL
  , begin = NULL
  , end = NULL
  , extent = NULL
  , tileH = NULL
  , tileV = NULL
  , SDSstring = NULL
  , job = NULL
  , datum = c("NODATUM", "NAD27", "NAD83", "WGS66", "WGS72", "WGS84")
  , zone = NULL
  , projPara = NULL
  , mosaic = TRUE
  , anonym = TRUE
  , ...
) {

  datum = match.arg(
    datum
  )
  
  dots = list(...)
    opts = do.call(combineOptions, dots)

    if (!opts$mrtOk)
    {
        stop("MRT path not set or MRT not installed on your system!")
    }
    
    product     <- getProduct(product, quiet = TRUE, collection = collection)
    tLimits          <- transDate(begin=begin,end=end)
    
    opts$localArcPath <- setPath(opts$localArcPath)
    opts$outDirPath   <- setPath(opts$outDirPath)
    
    if(!tolower(opts$dataFormat) %in% c('raw binary', 'hdf-eos', 'hdf4image','gtiff', 'geotiff'))
    {
        stop('dataFormat=\'',opts$dataFormat ,'\' is not supported by MRT (only \'raw binary\', \'HDF-EOS\' or \'GeoTiff\')')
    } 
    ext = checkMrtWriteDriver(opts$dataFormat)
    
    if (!inherits(extent, "MODISextent")) {
      extent = if (product@TYPE[1] == "Tile" |
                   (all(!is.null(extent) | !is.null(tileH) & !is.null(tileV)) &
                    product@TYPE[1]=="CMG")) {
        getTile(x = extent, tileH = tileH, tileV = tileV
                , outProj = "asIn"
                , pixelSize = "asIn")
      } else {
        NULL
      }
    }
    
    opts$resamplingType <- checkResamplingType(opts$resamplingType,tool="mrt",quiet=TRUE)
    opts$outProj        <- checkOutProj(opts$outProj,tool="mrt",quiet=TRUE)
    
    if (opts$outProj[1]=="asIn")
    {
      opts$outProj <- list(short="SIN",long="Sinusoidal")
    }
    cat("Output projection:", opts$outProj$long,"\n")    
    if (opts$outProj$short=="UTM")
    {
      zone = checkUTMZone(zone)
      
        if (is.null(zone)) 
        {
            cat("No UTM zone specified using MRT autodetection.\n")            
        } else 
        {
            cat("Using UTM zone:", zone,"\n")
        }
    }

    cat("Output pixel size:", opts$pixelSize,"\n")
    cat("Resampling method:", opts$resamplingType,"\n")
 
    if (!toupper(datum) %in% c("NAD27", "NAD83", "WGS66", "WGS72", "WGS84", "NODATUM"))
    {
      stop('Output "datum" must be one of: "NAD27", "NAD83", "WGS66", "WGS72", "WGS84" or "NODATUM"')
    } else {
      cat("Datum conversion:", datum, "\n")
    }
    
    if (is.null(projPara)) 
    {
        if(opts$outProj$short=="SIN") # maybe we should add other
        {
            projPara <- "6371007.18 0.00 0.00 0.00 0.00 0.00 0.00 0.00 86400.00 0.00 0.00 0.00 0.00 0.00 0.00"
        } else 
        {
            cat("No output projection parameters specified, using 'GEO' default settings!\n")
            projPara <- "0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0"
        }
    } else 
    {
        cat("Output projection parameters specified!\nUsing:\n ",projPara,"\n")
    }

    lst_product <- vector("list", length(product@PRODUCT))
    for (z in 1:length(product@PRODUCT))
    {
      
      # if (product@TYPE[z]=="CMG") 
      # {
      #   tileID="GLOBAL"
      #   ntiles=1 
      # } else 
      # {
      #   extent <- getTile(x=extent,tileH=tileH,tileV=tileV)
      #   ntiles    <- length(extent@tile)
      # }

        todo <- paste(product@PRODUCT[z],".",product@CCC[[product@PRODUCT[z]]],sep="")    
    
        lst_todo <- vector("list", length(todo))
        for(u in 1:length(todo))
        {
            if (is.null(job))
            {
                job <- paste(todo[u],"_",format(Sys.time(), "%Y%m%d%H%M%S"),sep="")    
                cat("No 'job' name specified, generated (date/time based)):",job,"\n")
            }
            outDir <- gsub("//", "/", file.path(opts$outDirPath,job,fsep="/"))
            dir.create(outDir, showWarnings = FALSE)

            ######################## along platform (TerraAqua)
            ftpdirs <- list()    ##  FRS: Fix provided by Ahmadou Dicko
            ftpdirs[[1]] <- as.Date(getStruc(product=strsplit(todo[u],"\\.")[[1]][1],collection=strsplit(todo[u],"\\.")[[1]][2],begin=tLimits$begin,end=tLimits$end,server=opts$MODISserverOrder[1])$dates)
            
            prodname <- strsplit(todo[u],"\\.")[[1]][1] 
            coll     <- strsplit(todo[u],"\\.")[[1]][2]
            
            avDates <- ftpdirs[[1]]
            avDates <- avDates[avDates!=FALSE]
            avDates <- avDates[!is.na(avDates)]        
            
            sel     <- as.Date(avDates)
            
            st = correctStartDate(tLimits$begin, sel, prodname, quiet = opts$quiet)
            us      <- sel >= st & sel <= tLimits$end
            
            if (sum(us,na.rm=TRUE)>0)
            {
                avDates <- avDates[us]

            ######################### along begin -> end date
                lst_ofile <- as.list(rep(NA, length(avDates)))
                for (l in 1:length(avDates))
                { 
                  files = unlist(
                    do.call(getHdf, append(list(product = product@PRODUCT[z]
                                                , collection = strsplit(todo[u], "\\.")[[1]][2]
                                                , begin = avDates[l], end = avDates[l]
                                                , tileH = extent@tileH, tileV = extent@tileV)
                                      , opts))
                  )

                    if (length(files)!=0)
                    {
                        mos <- mosaic
        
                        if (mos)
                        {
                            # if not all files available switch "off" mosaicking and process single files. Problematic in areas with tiles outside land!
                            if (sum(file.exists(files)) < length(extent@tile))
                            {
                                mos <- FALSE
                            } else {
                                mos <- TRUE
                            }
            
                        } else { 
                            mos <-  FALSE
                        }
            
                        if (mos)
                        {
                            v <- 1
                        } else {
                            v <- seq_along(files)
                        }
            
                        for (q in v)
                        {
                            w <- options("warn")
                            options(warn=-1)
                            if (is.null(SDSstring))
                            {
                                SDSstring <- rep(1,length(getSds(HdfName=files[q],method="mrt")$SDSnames))
                            }    
                
                            SDSstringIntern <- getSds(HdfName=files[q],SDSstring=SDSstring,method="mrt")
                            options(warn=w$warn)
                            
                            if (!opts$quiet && u == 1 && l == 1)
                            {
                                cat("\n#############################\nExtracing SDS:",SDSstringIntern$SDSnames,"#############################\n",sep="\n")
                            }
                            
                            if (mos)
                            {
                                TmpMosNam <- paste("TmpMosaic",makeRandomString(),".hdf",sep="")
                                ### in subset
                                paraname <- file.path(outDir,"/MRTgMosaic.prm",fsep="/") # create mosaic prm file
                                on.exit({
                                  try(unlink(paraname), silent = TRUE)
                                  try(unlink(paste(outDir,TmpMosNam,sep="/")), silent = TRUE)
                                })
                                
                                filename = file(paraname, open="wt")
                                on.exit(try(close(filename), silent = TRUE))
                                
                                write(paste("\"",files,"\"",sep='',collapse=' '), filename)
                                close(filename)
                        
                                # run mosaic
                                if (.Platform$OS=="unix")
                                {
                                    system(paste("mrtmosaic -i ",paraname," -o ",outDir,"/",TmpMosNam," -s '",SDSstringIntern$SDSstring,"'" ,sep=""))
                                } else 
                                {
                                    shell(paste("mrtmosaic -i \"",paraname,"\" -o \"", normalizePath(outDir) ,"\\",TmpMosNam,"\" -s \"",SDSstringIntern$SDSstring,"\"" ,sep=""))
                                }
                                unlink(paraname)
                                Sys.sleep(opts$wait) # without wait the skript can break here. "wait" is a try but it seams to work!!!
                            }
            
                            basenam <- strsplit(files[q],"/")[[1]]
                            basenam <- basenam[length(basenam)]
        
                            if (mos)
                            {
                                basenam <- paste(strsplit(basenam,"\\.")[[1]][c(1,2,4)],collapse=".")
                            } else {
                                basenam <- paste(strsplit(basenam,"\\.")[[1]][c(1,2,3,4)],collapse=".")    
                            }
        
                            if (!anonym)
                            {   
                                basenam <- paste(basenam,job,sep=".")
                            }
    
                            #### Write prm File
                            paraname <- paste(outDir,"/MRTgResample.prm",sep="")
                            on.exit({
                              try(unlink(paraname), silent = TRUE)
                            })
                            
                            filename = file(paraname, open="wt")
                            on.exit(try(close(filename), silent = TRUE))

                            if (mos)
                            {
                                write(paste('INPUT_FILENAME = "',outDir,"/",TmpMosNam,'"',sep=''), filename)
                            } else 
                            {
                                write(paste('SPECTRAL_SUBSET = ( ',SDSstringIntern$SDSstring,' )',sep=''), filename)
                                write(paste('INPUT_FILENAME = "',files[q],'"',sep=''), filename)
                            }
    
                            write('SPATIAL_SUBSET_TYPE = INPUT_LAT_LONG',filename)
    
                            if (!is.null(extent@extent))
                            {
                                write(paste('SPATIAL_SUBSET_UL_CORNER = (',extent@extent@ymax,' ',extent@extent@xmin,')',sep=''),filename)
                                write(paste('SPATIAL_SUBSET_LR_CORNER = (',extent@extent@ymin,' ',extent@extent@xmax,')',sep=''),filename)
                            }
                            if (opts$pixelSize!="asIn")
                            {
                                write(paste('OUTPUT_PIXEL_SIZE = ',opts$pixelSize,sep=''),filename) 
                            }    
                            write(paste('OUTPUT_FILENAME = ',outDir,"/",basenam,ext,sep=''),filename) 
                            write(paste('RESAMPLING_TYPE = ',opts$resamplingType,sep=''),filename)
                            
                            write(paste('OUTPUT_PROJECTION_TYPE = ',opts$outProj$short,sep=''),filename)
                            
                            if (opts$outProj$short=="UTM" && !is.null(zone))
                            {
                                write(paste('UTM_ZONE = ',zone,sep=''),filename)
                            }
        
                            if (!is.null(projPara))
                            {
                                write(paste('OUTPUT_PROJECTION_PARAMETERS = ( ',projPara,' )',sep=''),filename)
                            }
                            if (!is.null(datum))
                            {
                                write(paste('DATUM = ', datum,sep=''),filename)
                            }
                            close(filename)
    
                            if (.Platform$OS=="unix")
                            {
                                system(paste("resample -p ",paraname,sep=""))
                            } else {
                                shell(paste("resample -p \"",paraname,"\" -o \"",outDir,"/",basenam, ext,sep=""))
                            }
                            unlink(paraname)
    
                            if (mos)
                            {
                                unlink(paste(outDir,TmpMosNam,sep="/"))
                            }
                        }
                      
                      xtr = gsub(" ", "_", SDSstringIntern$SDSnames)
                      lst_ofile[[l]] = sapply(xtr, function(pttrn) {
                        list.files(outDir, full.names = TRUE, ignore.case = TRUE
                                   , pattern = paste(basenam, pttrn, ext, sep = ".*"))
                      })
                      lst_ofile[[l]] = as.character(lst_ofile[[l]])
                      
                    } else {
                      warning(paste0("No file found for date: ",avDates[l]))
                      lst_ofile[[l]] <- NA
                    }
                } # l, avDates
                
              names(lst_ofile) <- avDates
              lst_todo[[u]] <- lst_ofile
                
            } else {
              warning(paste("No", product@PRODUCT, "files found for the period from"
                            , tLimits$begin, "to", paste0(tLimits$end, ".")))
              lst_todo[[u]] <- NA
            }
        } # u   
        
        lst_product[[z]] <- lst_todo[[1]]
    } # z

    names(lst_product) <- paste(product@PRODUCT, product@CCC, sep = ".")
    return(lst_product)
}