runGdal <- function(product, collection=NULL, begin=NULL,end=NULL,
extent=NULL, tileH=NULL, tileV=NULL,
buffer=0, SDSstring=NULL, job=NULL,
checkIntegrity=TRUE, wait=0.5, quiet=FALSE,
exclList=NULL, resampList=NULL, nodataList=NULL,
scriptPath=NULL, ...)
{
opts <- combineOptions(...)
# debug:
# opts <- MODIS:::combineOptions()
# product="MOD11A1"; collection=NULL; begin='2013.06.01'; end='2013.06.05'; tileH=NULL; tileV=NULL
# buffer=0.04; SDSstring=NULL; job=NULL; checkIntegrity=TRUE; wait=0.5; quiet=FALSE; scriptPath=NULL
# exclList=list("Fpar_1km"="gt 100", "Lai_1km"="gt 100", "FparLai_QC"="255", "FparExtra_QC"="255", "FparStdDev_1km"="gt 100", "LaiStdDev_1km"="gt 100")
# resampList=list("Fpar_1km"="bilinear", "Lai_1km"="bilinear", "FparLai_QC"="mode", "FparExtra_QC"="mode", "FparStdDev_1km"="bilinear", "LaiStdDev_1km"="bilinear"), ...)
# extent <- raster::raster("~/d1/CO_UpperRioGrande/DOMAIN/geogrid_tmp.tif")
if(!opts$gdalOk)
{
stop("GDAL not installed or configured, read in '?MODISoptions' for help")
}
# absolutly needed
product <- getProduct(product,quiet=TRUE)
# check for optional python script path
if (is.null(scriptPath)) scriptPath = opts$gdalPath
# optional and if missing it is added here:
product$CCC <- getCollection(product,collection=collection)
tLimits <- transDate(begin=begin,end=end)
dataFormat <- toupper(opts$dataFormat)
if (dataFormat == 'RAW BINARY')
{
stop('in argument dataFormat=\'raw binary\', format not supported by GDAL (it is MRT specific) type: \'options("MODIS_gdalOutDriver")\' (column \'name\') to list available inputs')
}
if(dataFormat == 'HDF-EOS')
{
dataFormat <- "HDF4IMAGE"
} else if(dataFormat == 'GEOTIFF')
{
dataFormat <- "GTIFF"
}
if(is.null(opts$gdalOutDriver))
{
opts$gdalOutDriver <- gdalWriteDriver()
options("MODIS_gdalOutDriver"=opts$gdalOutDriver) # save for current session
}
if(dataFormat %in% toupper(opts$gdalOutDriver$name))
{
dataFormat <- grep(opts$gdalOutDriver$name, pattern=paste("^",dataFormat,"$",sep=""),ignore.case = TRUE,value=TRUE)
of <- paste0(" -of ",dataFormat)
extension <- MODIS:::getExtension(dataFormat)
} else
{
stop('in argument dataFormat=\'',opts$dataFormat,'\', format not supported by GDAL type: \'gdalWriteDriver()\' (column \'name\') to list available inputs')
}
#### settings with messages
# output pixel size in output proj units (default is "asIn", but there are 2 chances of changing this argument: pixelSize, and if extent comes from a Raster* object.
if (product$TYPE[1]=="Tile" | (all(!is.null(extent) | !is.null(tileH) & !is.null(tileV)) & product$TYPE[1]=="CMG"))
{
extent <- getTile(extent=extent, tileH=tileH, tileV=tileV, buffer=buffer)
} else
{
extent <- NULL
}
#### outProj
t_srs <- NULL
cat("########################\n")
if(!is.null(extent$target$outProj))
{
outProj <- MODIS:::checkOutProj(extent$target$outProj,tool="GDAL")
cat("outProj = ",outProj ," (Specified by raster*/spatial* object)\n")
} else
{
outProj <- MODIS:::checkOutProj(opts$outProj,tool="GDAL")
cat("outProj = ",outProj,"\n")
}
if (outProj == "asIn")
{
if (product$SENSOR[1]=="MODIS")
{
if (product$TYPE[1]=="Tile")
{
outProj <- "+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs"
} else
{
outProj <- "+proj=longlat +ellps=clrk66 +no_defs" # CMG proj
}
} else if (product$SENSOR[1]=="SRTM")
{
outProj <- "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"
}
}
t_srs <- paste0(' -t_srs ',shQuote(outProj))
#### pixelSize
if(!is.null(extent$target$pixelSize))
{
pixelSize <- extent$target$pixelSize
cat("pixelSize = ",pixelSize ," (Specified by raster* object)\n")
} else
{
pixelSize <- opts$pixelSize
cat("pixelSize = ",pixelSize,"\n")
}
tr <- NULL
if (pixelSize[1]!="asIn")
{
if (length(pixelSize)==1)
{
tr <- paste(" -tr",pixelSize,pixelSize)
} else
{
tr <- paste0(" -tr ", paste0(pixelSize,collapse=" "))
}
}
#### resamplingType
opts$resamplingType <- MODIS:::checkResamplingType(opts$resamplingType, tool="gdal")
cat("resamplingType = ", opts$resamplingType,"\n")
rt <- paste0(" -r ",opts$resamplingType)
#### inProj (s_srs)
if (product$SENSOR[1]=="MODIS")
{
if (product$TYPE[1]=="Tile")
{
s_srs <- paste0(' -s_srs ',shQuote("+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs"))
} else
{
s_srs <- paste0(' -s_srs ',shQuote("+proj=longlat +ellps=clrk66 +no_defs"))
}
} else if (product$SENSOR[1]=="SRTM")
{
s_srs <- paste0(' -s_srs ',shQuote("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"))
}
#### te (target @extent)
te <- NULL # if extent comes from tileV/H
if (!is.null(extent$target$extent)) # all extents but not tileV/H
{
if (is.null(extent$target$outProj)) # map or list extents (always LatLon)
{
rx <- raster(extent$target$extent,crs="+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")
rx <- projectExtent(rx,outProj)
rx <- extent(rx)
} else
{
rx <- extent$target$extent
}
te <- paste(" -te", rx@xmin, rx@ymin, rx@xmax, rx@ymax)
}
if (is.null(extent$target))
{
if(!is.null(extent$extent))
{
rx <- raster(extent$extent,crs="+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")
rx <- projectExtent(rx,outProj)
rx <- extent(rx)
te <- paste(" -te", rx@xmin, rx@ymin, rx@xmax, rx@ymax)
}
}
#### generate non-obligatory GDAL arguments
# GeoTiff BLOCKYSIZE and compression. See: http://www.gdal.org/frmt_gtiff.html
if(is.null(opts$blockSize))
{
bs <- NULL
} else
{
opts$blockSize <- as.integer(opts$blockSize)
bs <- paste0(" -co BLOCKYSIZE=",opts$blockSize)
}
# compress output data
if(is.null(opts$compression))
{
cp <- " -co compress=lzw -co predictor=2"
} else if (isTRUE(opts$compression))
{
cp <- " -co compress=lzw -co predictor=2"
} else
{
cp <- NULL
}
####
if (quiet)
{
q <- " -q"
} else
{
q <- NULL
}
# Setup exclusion range qualifiers
compStr1 <- list("none"="!=", "gt"="<=", "lt"=">=")
compStr2 <- list("none"="==", "gt"=">", "lt"="<")
for (z in seq_along(product$PRODUCT))
{ # z=1
todo <- paste(product$PRODUCT[z],".",product$CCC[z],sep="")
if(z==1)
{
if (is.null(job))
{
job <- paste0(todo[1],"_",format(Sys.time(), "%Y%m%d%H%M%S"))
cat("Output directory = ",paste0(normalizePath(opts$outDirPath,"/",mustWork=FALSE),"/",job)," (no 'job' name specified, generated (date/time based))\n")
} else
{
cat("Output Directory = ",paste0(normalizePath(opts$outDirPath,"/",mustWork=FALSE),"/",job),"\n")
}
cat("########################\n")
outDir <- file.path(opts$outDirPath,job,fsep="/")
dir.create(outDir,showWarnings=FALSE,recursive=TRUE)
}
for(u in seq_along(todo))
{ # u=1
ftpdirs <- list()
server <- ifelse (product$SOURCE[z]=="NSIDC", "NSIDC", opts$MODISserverOrder[1])
ftpdirs[[1]] <- as.Date(MODIS:::getStruc(product=strsplit(todo[u],"\\.")[[1]][1],collection=strsplit(todo[u],"\\.")[[1]][2],begin=tLimits$begin,end=tLimits$end,server=server)$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)
us <- sel >= tLimits$begin & sel <= tLimits$end
if (sum(us,na.rm=TRUE)>0)
{
avDates <- avDates[us]
for (l in seq_along(avDates))
{ # l=1
files <- unlist(
getHdf(product=prodname, collection=coll, begin=avDates[l], end=avDates[l],
tileH=extent$tileH, tileV=extent$tileV, checkIntegrity=checkIntegrity,
stubbornness=opts$stubbornness, MODISserverOrder=opts$MODISserverOrder)
)
files <- files[basename(files)!="NA"] # is not a true NA so it need to be like that na not !is.na()
if(length(files)>0)
{
w <- getOption("warn")
options("warn"= -1)
SDS <- list()
for (zz in seq_along(files))
{ # get all SDS names for one chunk
SDS[[zz]] <- MODIS:::getSds(HdfName=files[zz], SDSstring=SDSstring, method="GDAL")
}
options("warn"= w)
if (!exists("NAS"))
{
NAS <- MODIS:::getNa(SDS[[1]]$SDS4gdal)
}
for (i in seq_along(SDS[[1]]$SDSnames))
{ # i=1
outname <- paste0(paste0(strsplit(basename(files[1]),"\\.")[[1]][1:2],collapse="."),
".", gsub(SDS[[1]]$SDSnames[i],pattern=" ",replacement="_"), extension)
gdalSDS <- sapply(SDS,function(x){x$SDS4gdal[i]}) # get names of layer 'o' of all files (SDS)
naID <- which( names(NAS) == SDS[[1]]$SDSnames[i] )
if(length(naID)>0)
{
srcnodata <- paste0(" -srcnodata ",NAS[[naID]])
dstnodata <- paste0(" -dstnodata ",NAS[[naID]])
} else {
srcnodata <- NULL
dstnodata <- NULL
}
# Figure out exclusion ranges if provided
exclID<-NULL
exclStr<-NULL
if (!is.null(exclList)) {
exclID <- which( names(exclList) == SDS[[1]]$SDSnames[i] )
}
if(length(exclID)>0) {
exclStr <- exclList[[SDS[[1]]$SDSnames[i]]]
if (length(unlist(strsplit(exclStr, split=" "))) == 1) {
exclQual <- 'none'
exclNum <- unlist(strsplit(exclStr, split=" "))[1]
} else if (length(unlist(strsplit(exclStr, split=" "))) == 2) {
exclQual <- unlist(strsplit(exclStr, split=" "))[1]
exclNum <- unlist(strsplit(exclStr, split=" "))[2]
} else {
stop("Invalid exclList entry.")
return()
}
}
# Figure out resampling method if provided
resampID<-NULL
resampStr<-NULL
if( SDS[[1]]$SDSnames[i] %in% names(resampList) ) {
resampStr <- resampList[[SDS[[1]]$SDSnames[i]]]
} else {
resampStr <- opts$resamplingType
}
rt <- paste0(" -r ", resampStr)
if(length(grep(todo,pattern="M.D13C2\\.005"))>0)
{
if(i==1)
{
cat("\n###############\nM.D13C2.005 is likely to have a problem in metadata extent information, it is corrected on the fly\n###############\n")
}
ranpat <- MODIS:::makeRandomString(length=21)
randomName <- paste0(outDir,"/deleteMe_",ranpat,".tif")
#on.exit(unlink(list.files(path=outDir,pattern=ranpat,full.names=TRUE),recursive=TRUE))
for(ix in seq_along(gdalSDS))
{
cmd1 <- paste0(opts$gdalPath,"gdal_translate -a_nodata ",NAS[[naID]]," '",gdalSDS[ix],"' '",randomName[ix],"'")
cmd2 <- paste0(scriptPath,"gdal_edit.py -a_ullr -180 90 180 -90 '",randomName[ix],"'")
cmd3 <- paste0(opts$gdalPath, )
if (.Platform$OS=="unix")
{
system(cmd1,intern=TRUE)
system(cmd2,intern=TRUE)
} else
{
shell(cmd1,intern=TRUE)
shell(cmd2,intern=TRUE)
}
}
gdalSDS <- randomName
}
if (.Platform$OS=="unix")
{
ifile <- paste0(gdalSDS,collapse="' '")
ofile <- paste0(outDir, '/', outname)
# New intermediate step to convert excluded value ranges to "nodata" val before running main conversion.
# This should aid using interpolation methods other than near and mode.
if (length(exclID)>0) {
outList<-c()
for (m in 1:length(gdalSDS)) {
ranpat2 <- MODIS:::makeRandomString(length=21)
randomName <- paste0(outDir,"/deleteMe_",ranpat2,".tif")
#on.exit(unlink(list.files(path=outDir,pattern=ranpat2,full.names=TRUE),recursive=TRUE))
if (length(naID)==0) { #missing/non-existent nodata value case
stop(paste0("Missing no data value for variable ", SDS[[1]]$SDSnames[i], ". Remove exclusion values for this variable and try again."))
}
else {
cmd_pre <- paste0(scriptPath,"gdal_calc.py ",
"-A ", "'", gdalSDS[m], "'",
" --NoDataValue=", as.character(NAS[[naID]]),
" --outfile=", "'", randomName, "'",
" --calc='A*(A", compStr1[[exclQual]],
as.character(exclNum), ")+(",
as.character(NAS[[naID]]), ")*(A", compStr2[[exclQual]],
as.character(exclNum), ")'", " --overwrite")
}
if (!quiet) {print(cmd_pre)}
system(cmd_pre)
outList <- c(outList, randomName)
}
ifile <- paste0(outList, collapse="' '")
}
# Original
cmd <- paste0(opts$gdalPath,
"gdalwarp",
s_srs,
t_srs,
of,
te,
tr,
cp,
bs,
rt,
q,
srcnodata,
dstnodata,
" -overwrite",
" -multi",
" \'", ifile,"\'",
" ",
ofile
)
cmd <- gsub(x=cmd,pattern="\"",replacement="'")
if (!quiet) {print(cmd)}
system(cmd)
} else # windows
{
if (length(exclID)>0) print("Exclusion ranges not implemented in Windows. Continuing with default single value.")
cmd <- paste0(opts$gdalPath,"gdalwarp")
# ifile <- paste(shortPathName(gdalSDS),collapse='\" \"',sep=' ')
# ofile <- shortPathName(paste0(normalizePath(outDir), '\\', outname))
ofile <- paste0(outDir, '/', outname)
ifile <- paste0(gdalSDS,collapse='" "')
# GDAL < 1.8.0 doesn't support ' -overwrite'
if(file.exists(ofile))
{
invisible(file.remove(ofile))
}
shell(
paste(cmd,
s_srs,
t_srs,
of,
te,
tr,
cp,
bs,
rt,
q,
srcnodata,
dstnodata,
' -multi',
' \"', ifile,'\"',
' \"', ofile,'\"',
sep = '')
)
}
if(length(grep(todo,pattern="M.D13C2\\.005"))>0)
{
unlink(list.files(path=outDir,pattern=ranpat,full.names=TRUE),recursive=TRUE)
}
if (length(exclID)>0) {
unlink(outList, recursive=TRUE)
#unlink(list.files(path=outDir,pattern=ranpat2,full.names=TRUE),recursive=TRUE)
#unlink(list.files(path=outDir,pattern=glob2rx("deleteMe_*.tif"),full.names=TRUE),recursive=TRUE)
}
} # END end for i (SDSnames)
} else # end if(length(files)>0)
{
warning(paste0("No file found for date: ",avDates[l]))
}
} # end for l (dates)
} # end if (sum(us,na.rm=TRUE)>0)
} # end for u (collection)
} #end for z (product)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.