#' Get MODIS data and process to match geogrid
#'
#' \code{GetMODIS} downloads, mosaics, and resamples MODIS data to match input
#' geogrid.
#'
#' \code{GetMODIS} reads a geogrid file and parameters on MODIS product and data
#' range and downloads MODIS tiles where they do not exist (using the R
#' \href{http://r-forge.r-project.org/projects/modis/}{MODIS} package), mosaics
#' where necessary, clips, reprojects, and resamples (using nearest neighbor) to
#' match the geogrid. Results in a set of TIF files per the
#' \href{http://r-forge.r-project.org/projects/modis/}{MODIS} package
#' specifications.
#'
#' Please see documentation on the R
#' \href{http://r-forge.r-project.org/projects/modis/}{MODIS} package for
#' details on required installs and workspace setup. This tool builds off of the
#' \href{http://r-forge.r-project.org/projects/modis/}{MODIS}
#' \code{\link[MODIS]{runGdal}} tool, so follows
#' \href{http://r-forge.r-project.org/projects/modis/}{MODIS} file directory
#' structure (see \code{\link[MODIS]{MODISoptions}}). This tool requires a local
#' \href{http://www.gdal.org/}{GDAL} installation in addition to the required R
#' packages.
#'
#' NOTE: This tool currently only works for geogrid files in Lambert Conformal
#' Conic projection OR a GDAL-friendly georeferenced TIF file.
#'
#' @param extent The pathname to the geogrid file (i.e., geo_em.d01.nc)
#' OR the path to a georeferenced TIF file of the domain OR the path to
#' a shapefile of the desired extent OR a standard country name that can
#' be identified in the mapdata package.
#' @param prodName The MODIS product name to download/process. Run the
#' \href{http://r-forge.r-project.org/projects/modis/}{MODIS} package
#' getProducts() for a complete list of supported products.
#' @param outDir Directory name to store processed TIF files. This is the
#' equivalent to the
#' \href{http://r-forge.r-project.org/projects/modis/}{MODIS} package's "job"
#' name. This is a directory name only and NOT a full path. The directory will
#' be created in the preset
#' \href{http://r-forge.r-project.org/projects/modis/}{MODIS} package
#' outDirPath.
#' @param begin Date string for the start date to download/process MODIS tiles.
#' The date string should follow
#' \href{http://r-forge.r-project.org/projects/modis/}{MODIS} package
#' convention (e.g., "2011.06.01").
#' @param end Date string for the end date to download/process MODIS tiles. The
#' date string should follow
#' \href{http://r-forge.r-project.org/projects/modis/}{MODIS} package
#' convention (e.g., "2011.06.01").
#' @param collection From MODIS::runGdal help: "Default is to download most
#' recent collection version. See: ?getCollection"
#' @param buffer From MODIS::runGdal help: "Numeric. Buffer [in units of the
#' 'outProj'] around the specified extent. See: ?getTile" (DEFAULT=0.04, which
#' seems to cover MODIS tile corners sufficiently)
#' @param SDSstring From MODIS::runGdal help: "Default is extract all SDS
#' (layers). See: ?getSds."
#' @param exclList List pairing SDS name (e.g., "Lai_1km", "FparLai_QC") with
#' values to exclude from interpolation (i.e., convert to "nodata"). For
#' example, for the Lai_1km product, valid range is 0-100, values 249-255 are
#' fill values, but only one value (255) is reported as _FillValue in the HDF
#' metadata. The exclude list should specify value ranges to exclude, and
#' should be in the form: \deqn{list("SDS name"="exclude range", ...)} where
#' "exclude range" can be a single number (e.g., 254), "gt <somenumber>" to
#' exclude all values over some threshold (e.g., "gt 100"), or "lt
#' <somenumber>" to exclude all values less than some threshold (e.g., "lt
#' 0").
#' @param resampList List pairing SDS name (e.g., "Lai_1km", "FparLai_QC") with
#' reasmpling method to use. GDAL possible options (depends on version):
#' 'near', 'bilinear', 'cubic', 'cubicspline', 'lanczos', 'mode', 'average'
#' (DEFAULT='near', unless otherwise specified in
#' \code{MODISoptions(resamplingType="<somemethod>")}). The resampling list
#' should be in the form: \deqn{list("SDS name"="method", ...)} where "method"
#' is one of the GDAL methods listed above.
#' @param quiet From MODIS::runGdal help: Logical, Default FALSE. Some progress
#' informations.
#' @param scriptPath OPTIONAL path to where gdal_calc.py script lives in case
#' it is not in the same default path as gdal executables
#' @return Empty
#'
#' @examples
#' ## First, specify the target directories for the MODIS package.
#' \dontrun{
#' MODISoptions(localArcPath="/d1/WRF_Hydro/RS/MODIS_ARC",
#' outDirPath="/d1/WRF_Hydro/RS/MODIS_ARC/PROCESSED",
#' stubbornness="low")
#'
#' ## Then, use GetMODIS to download the MODIS MOD15A2 (FPAR/LAI) product for
#' ## all tiles that overlap our Fourmile Creek domain (as specified by the
#' ## geo_em.d01.nc file), mosaic if multiple tiles, clip to the domain extent,
#' ## and resample to the domain grid cells. The raw HDF files will be stored in
#' ## /d1/WRF_Hydro/RS/MODIS_ARC/ and the final processed TIF files will be
#' ## stored in /d1/WRF_Hydro/RS/MODIS_ARC/PROCESSED/Fourmile_LAI/.
#'
#' GetMODIS(extent="/d1/WRF_Hydro/Fourmile_fire/DOMAIN/geo_em.d01.nc",
#' prodName="MOD15A2", outDir="Fourmile_LAI",
#' begin="2011.01.01", end="2011.01.31",
#' 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"))
#' }
#' @keywords IO
#' @concept MODIS dataGet
#' @family MODIS
#' @export
GetMODIS <- function(extent, prodName, outDir, begin=NULL, end=NULL,
collection=NULL, buffer=0.04,
SDSstring=NULL,
exclList=NULL, resampList=NULL, quiet=FALSE, scriptPath=NULL) {
# Check packages
if (!(require("rgdal") & require("raster") & require("ncdf4") & require("MODIS"))) {
stop("Required packages not found. Must have R packages: rgdal (requires GDAL system install), raster, ncdf4, and MODIS")
}
# Get paths
locPath <- paste0(options("MODIS_localArcPath"))
if (!file.exists(locPath)) {
dir.create(locPath, showWarnings = FALSE)
}
# Get geogrid and projection info
if (inherits(try(suppressWarnings(ncdf4::nc_open(extent))), "try-error")) {
if (file.exists(extent)) {
if ( grepl(".shp", extent, ignore.case = TRUE) ) {
message("Assuming file is a shapefile.")
hgt.r <- extent
} else if ( grepl(".tif", extent, ignore.case = TRUE) ) {
message("Assuming file is a georeferenced TIF.")
hgt.r <- raster::raster(extent)
} else {
stop("Must provide a .shp or .tif file.")
}
} else {
message("Assuming file is a country name. Requires 'mapdata' package.")
hgt.r <- extent
}
} else {
message("Assuming file is a geogrid in netcdf format.")
geogrd.nc <- ncdf4::nc_open(extent)
map_proj <- ncdf4::ncatt_get(geogrd.nc, varid=0, attname="MAP_PROJ")$value
cen_lat <- ncdf4::ncatt_get(geogrd.nc, varid=0, attname="CEN_LAT")$value
cen_lon <- ncdf4::ncatt_get(geogrd.nc, varid=0, attname="STAND_LON")$value
truelat1 <- ncdf4::ncatt_get(geogrd.nc, varid=0, attname="TRUELAT1")$value
truelat2 <- ncdf4::ncatt_get(geogrd.nc, varid=0, attname="TRUELAT2")$value
if (map_proj==1) {
geogrd.crs <- paste0("+proj=lcc +lat_1=", truelat1, " +lat_2=", truelat2,
" +lat_0=", cen_lat, " +lon_0=", cen_lon,
" +x_0=0 +y_0=0 +a=6370000 +b=6370000 +units=m +no_defs")
} else {
stop('Error: Projection type not supported (currently this tool only works for Lambert Conformal Conic projections).')
}
dx <- ncdf4::ncatt_get(geogrd.nc, varid=0, attname="DX")$value
dy <- ncdf4::ncatt_get(geogrd.nc, varid=0, attname="DY")$value
if ( dx != dy ) {
stop(paste0('Error: Asymmetric grid cells not supported. DX=', dx, ', DY=', dy))
}
# Create a readable TIF from geogrid
ExportGeogrid(extent, "HGT_M", paste0(locPath, "/geogrid_tmp.tif"))
hgt.r <- raster::raster(paste0(locPath, "/geogrid_tmp.tif"))
system(paste0("rm ", paste0(locPath, "/geogrid_tmp.tif")))
}
# Run the download & processing
mod.list <- MODIS::runGdal(product=prodName, collection=collection,
begin=begin, end=end,
extent=hgt.r, buffer=buffer,
SDSstring=SDSstring, job=outDir,
exclList=exclList, resampList=resampList,
quiet=quiet, scriptPath=scriptPath)
}
#' Convert a set of MODIS images to a raster stack and, optionally, a NetCDF
#' file.
#'
#' \code{ConvertRS2Stack} takes a set of pre-processed MODIS TIF files and
#' creates a raster stack and, optionally, a NetCDF file.
#'
#' \code{ConvertRS2Stack} scans the specified directory and imports
#' pre-processed MODIS TIF files matching the specified expression and combines
#' them into an R raster stack object and, optionally, an output NetCDF file.
#' Files in the input directory should be already processed through the
#' \code{\link{GetMODIS}} tool or follow the same file naming convention used by
#' the MODIS \code{\link[MODIS]{runGdal}} tool. See the R
#' \href{http://r-forge.r-project.org/projects/modis/}{MODIS} package for more
#' on those specifications.
#'
#' @param inPath The path to the directory that holds the already processed
#' MODIS TIF files.
#' @param matchStr The regular expression for filename matching (e.g.,
#' "*Lai_1km.tif").
#' @param begin Date string for the start date to include. The date string
#' should follow \href{http://r-forge.r-project.org/projects/modis/}{MODIS}
#' package convention (e.g., "2011.06.01"). (DEFAULT=NULL, all files are
#' processed).
#' @param end Date string for the end date to include. The date string should
#' follow \href{http://r-forge.r-project.org/projects/modis/}{MODIS} package
#' convention (e.g., "2011.06.01"). (DEFAULT=NULL, all files are processed).
#' @param noData Value for specifying "no data" per the MODIS product. Should be
#' combined with a qualifier. For example, for the MOD15A2 (LAI) product,
#' valid value range is 0-100. So the noData value should be set to 100 and
#' the noDataQual to "max". (DEFAULT=NULL, no additional "no data" screening
#' is applied)
#' @param noDataQual Qualifier for specifying "no data" per the MODIS product.
#' Should be combined with a nodata value. For example, for the MOD15A2 LAI
#' product, valid value range is 0-100. So the noData value should be set to
#' 100 and the noDataQual to "max". Options are: "exact" (value == noData are
#' converted to NA), "min" (value < noData are converted to NA), "max" (value
#' > noData are converted to NA). (DEFAULT="exact", only exact matches are
#' converted to NA)
#' @param valScale The scale factor to apply to the image values per the MODIS
#' product. The adjustment is applied as VALUE * valScale + valAdd. For
#' example, for the MOD15A2 LAI product, the scale factor is 0.1. (DEFAULT =
#' 1)
#' @param valAdd The addition value to apply to the image values per the MODIS
#' product. The adjustment is applied as VALUE * valScale + valAdd. (DEFAULT =
#' 0)
#' @param outFile OPTIONAL name for an output NetCDF file. A NetCDF file will
#' only be created if this file name is provided. Images will be exported
#' after the "no data" and scale adjustments are made. If you want to do
#' smoothing or other time series processing, do not export a NetCDF file here
#' but do processing on the raster stack and then use ConvertStack2NC to
#' export the processed brick to a NetCDF file.
#' @param varName Name for the NetCDF export variable. Only required if outFile
#' is provided.
#' @param varUnit Units for the NetCDF export variable. Only required if outFile
#' is provided.
#' @param varLong Long name for the NetCDF export variable. Only required if
#' outFile is provided.
#' @param varNA Value to set for "NA" or "no data". Default is -1.e+36. Only
#' required if outFile is provided.
#' @param pos1 From MODIS orgTime: Start position of date in the filename
#' (DEFAULT=10).
#' @param pos2 From MODIS orgTime: End position of date in the filename
#' (DEFAULT=16).
#' @param format From MODIS orgTime: How is the date formatted in the file.
#' Default is "\%Y\%j" ('YYYYDDD'). Read ?as.Date for for more information.
#' @return A raster stack.
#' @examples
#' ## Import the already processed LAI TIF images into a raster stack. Use
#' ## the full time series of images.
#'
#' \dontrun{
#' lai.b <- ConvertRS2Stack("/d6/adugger/WRF_Hydro/RS/MODIS_ARC/PROCESSED/BCNED_LAI",
#' "*Lai_1km.tif", noData=100, noDataQual="max",
#' valScale=0.1, valAdd=0)
#'
#' ## Export a subset of the already processed LAI TIF images into an output netcdf file
#'
#' lai.b <- ConvertRS2Stack("/d6/adugger/WRF_Hydro/RS/MODIS_ARC/PROCESSED/BCNED_LAI",
#' "*Lai_1km.tif", begin="2011.06.01", end="2011.06.30",
#' noData=100, noDataQual="max", valScale=0.1, valAdd=0,
#' outFile="BCNED_LAI.nc", varName="LAI",
#' varUnit="(m^2)/(m^2)", varLong="Leaf area index")
#' }
#' @keywords IO
#' @concept MODIS dataMgmt
#' @family MODIS
#' @export
ConvertRS2Stack <- function(inPath, matchStr, begin=NULL, end=NULL,
noData=NULL, noDataQual="exact", valScale=1, valAdd=0,
outFile=NULL, varName=NULL, varUnit=NULL, varLong=NULL,
varNA=(-1.e+36),
pos1=10, pos2=16, format="%Y%j") {
# Get file list
if (!is.null(begin) & !is.null(end)) {
beginDt <- format(as.POSIXct(begin, format="%Y.%m.%d", tz="UTC"),
"%Y%j", tz="UTC")
endDt <- format(as.POSIXct(end, format="%Y.%m.%d", tz="UTC"),
"%Y%j", tz="UTC")
timeInfo <- MODIS::orgTime(MODIS::preStack(path=inPath, pattern=matchStr),
begin=beginDt, end=endDt, pillow=0,
pos1=pos1, pos2=pos2, format=format)
} else {
timeInfo <- MODIS::orgTime(MODIS::preStack(path=inPath, pattern=matchStr),
pillow=0, pos1=pos1, pos2=pos2, format=format)
}
rsFilelist <- MODIS::preStack(path=inPath, pattern=matchStr, timeInfo=timeInfo)
#rsFilelist <- list.files(path=inPath, pattern=glob2rx(matchStr), full.names=TRUE)
# Loop through files
n <- 1
dtInts <- c()
dtNames <- c()
for (rsFile in rsFilelist) {
# Check dates
fileStr <- sub(inPath, rsFile, replacement="")
dtStr <- substr(unlist(strsplit(fileStr,"[.]"))[2],
nchar(unlist(strsplit(fileStr,"[.]"))[2])-(7-1),
nchar(unlist(strsplit(fileStr,"[.]"))[2]))
#dtPOSIXct <- as.POSIXct(dtStr, format="%Y%j", tz="UTC")
#begPOSIXct <- if(is.null(begin)) { dtPOSIXct } else { as.POSIXct(begin, format="%Y.%m.%d", tz="UTC") }
#endPOSIXct <- if(is.null(end)) { dtPOSIXct } else { as.POSIXct(end, format="%Y.%m.%d", tz="UTC") }
#if ( (dtPOSIXct >= begPOSIXct) & (dtPOSIXct <= endPOSIXct) ) {
rsRast <- raster::raster(rsFile)
# Remove MODIS nodata values
if (!is.null(noData)) {
if (noDataQual == "min") {
rsRast[rsRast[] < noData]<-NA
} else if (noDataQual == "max") {
rsRast[rsRast[] > noData]<-NA
} else {
rsRast[rsRast[] == noData]<-NA
}
}
# Apply MODIS scaling
if ((valScale != 1) | (valAdd != 0)) {
rsRast <- rsRast * valScale + valAdd
}
# Track dates
dtInts[n] <- as.integer(difftime(as.POSIXct(dtStr, format="%Y%j", tz="UTC"),
as.POSIXct("198001", format="%Y%j", tz="UTC",
units="days")))
dtNames[n] <- paste0("DT", format(as.POSIXct(dtStr, format="%Y%j", tz="UTC"),
"%Y.%m.%d"))
if (n==1) {
rsStack <- raster::stack(rsRast)
} else {
rsStack <- raster::addLayer(rsStack, rsRast)
}
n <- n+1
# } # end date check
} # end for loop
names(rsStack) <- dtNames
if (!is.null(outFile)) {
ConvertStack2NC(rsStack, outFile, varName, varUnit, varLong, varNA)
return(rsStack)
}
return(rsStack)
}
#' Convert a raster stack to a NetCDF file.
#'
#' \code{ConvertStack2NC} takes a raster stack of RS images and outputs a NetCDF
#' file with a time dimension and specified variable name & metadata.
#'
#' \code{ConvertStack2NC} converts a raster stack to an output NetCDF file. The
#' raster stack should be already processed through the
#' \code{\link{ConvertRS2Stack}} tool or follow the same layer (date) naming
#' convention.
#'
#' @param inStack The name of the raster stack to export.
#' @param outFile Name for an output NetCDF file.
#' @param varName Name for the NetCDF export variable.
#' @param varUnit Units for the NetCDF export variable.
#' @param varLong Long name for the NetCDF export variable.
#' @param varNA Value to set for "NA" or "no data". Default is -1.e+36.
#' @param flipY Reverse y indices, e.g., to match geogrid S->N
#' orientation (DEFAULT=TRUE).
#' @return NULL
#'
#' @examples
#' ## Export the raster stack of LAI images created through ConvertRS2Stack
#' ## to a NetCDF file. Use the full time series of images.
#'
#' \dontrun{
#' ConvertStack2NC(lai.b, outFile="BCNED_LAI.nc", varName="LAI",
#' varUnit="(m^2)/(m^2)", varLong="Leaf area index")
#' }
#' @keywords IO
#' @concept MODIS dataMgmt
#' @family MODIS
#' @export
ConvertStack2NC <- function(inStack, outFile=NULL, varName=NULL, varUnit=NULL,
varLong=NULL, varNA=-1.e+36, flipY=TRUE) {
# Get dates
dtInts <- c()
dtNames <- names(inStack)
for (i in 1:length(dtNames)) {
dtStr <- sub("DT", dtNames[i], replacement="")
dtInts[i] <- as.integer(difftime(as.POSIXct(dtStr, format="%Y.%m.%d", tz="UTC"),
as.POSIXct("198001", format="%Y%j", tz="UTC",
units="days")))
}
# Output NetCDF file
raster::writeRaster(inStack, outFile, "CDF", overwrite=TRUE,
varname=varName, varunit=varUnit, longname=varLong,
xname="west_east", yname="south_north", zname="Time",
zunit="days since 1980-01-01",
bylayer=FALSE, NAflag=varNA)
# Set the time variable
ncFile <- ncdf4::nc_open(outFile, write=TRUE)
ncdf4::ncvar_put(ncFile, "Time", dtInts)
ncdf4::nc_close(ncFile)
if (flipY) {
system(paste0('ncpdq -O -a -south_north ', outFile, ' ', outFile))
}
}
#' Run MODIS-R Whittaker smoothing over pre-processed raster stack.
#'
#' \code{SmoothStack} takes a raster stack of RS images and outputs a smoothed
#' raster brick over the same time period.
#'
#' \code{SmoothStack} converts a raster stack of RS images (as processed through
#' \code{\link{ConvertRS2Stack}}) and calls the MODIS-R
#' \code{\link[MODIS]{whittaker.raster}} smoothing function to generate a
#' smoothed raster brick over the same time series as the input. All function
#' parameters are per the \code{\link[MODIS]{whittaker.raster}} function except
#' we force the timeInfo to be derived from the input raster stack so the names
#' match and the smoothed brick can be used in other tools. The
#' \code{\link[MODIS]{whittaker.raster}} tool also exports a set of smoothed
#' TIFs, so also specify an output file directory.
#'
#' @param inStack The name of the raster stack to smooth.
#' @param w FROM \code{\link[MODIS]{whittaker.raster}}: In case of MODIS
#' composite the 'VI_Quality' raster-Brick, Stack or filenames. Use preStack
#' functionality to ensure the right input.
#' @param t FROM \code{\link[MODIS]{whittaker.raster}}: In case of MODIS
#' composite the 'composite_day_of_the_year' raster-Brick, Stack or filenames.
#' Use preStack functionality to ensure the right input.
#' @param lambda FROM \code{\link[MODIS]{whittaker.raster}}: _Yearly_ lambda
#' value passed to ?ptw:::wit2. If set as character (i.e. lambda="600"), it is
#' not adapted to the time serie length but used as a fixed value (see
#' details). High values = stiff/rigid spline
#' @param nIter FROM \code{\link[MODIS]{whittaker.raster}}: Number of iteration
#' for the upper envelope fitting.
#' @param outputAs FROM \code{\link[MODIS]{whittaker.raster}}: Character.
#' Organisation of output files: single each date one RasterLayer, yearly a
#' RasterBrick for each year, one one RasterBrick for the entire time-serie.
#' @param collapse FROM \code{\link[MODIS]{whittaker.raster}}: logical, if TRUE
#' the input data is treated as _1_single_Year_ collapsing the data using the
#' Julian date information without the year.
#' @param outDirPath FROM \code{\link[MODIS]{whittaker.raster}}: Output path
#' default is the current directory.
#' @param removeOutlier FROM \code{\link[MODIS]{whittaker.raster}}: Logical. See
#' details
#' @param outlierThreshold FROM \code{\link[MODIS]{whittaker.raster}}: Numerical
#' in the same unit as vi, used for outliers removal. See details
#' @param mergeDoyFun FROM \code{\link[MODIS]{whittaker.raster}}: Especially
#' when using argument collapse=TRUE, multiple measurements for one day can be
#' present, here you can choose how those values are merged to one single
#' value: "max" use the highest value, "mean" or "weighted.mean" use the mean
#' if no weighting "w" is available and weighted.mean if it is.
#' @param ... nFROM \code{\link[MODIS]{whittaker.raster}}: Arguments passed to
#' ?writeRaster (except filename is automatic), NAflag, datatype,
#' overwrite,...
#' @return raster brick of smoothed images
#'
#' @examples
#' ## Take the raster stack of LAI images created through ConvertRS2Stack and
#' ## apply a smoothing filter that also removes outliers, which we specify to
#' ## be more than 0.5 LAI from the smoothed value.
#' \dontrun{
#' lai.b.sm <- SmoothStack(lai.b,
#' outDirPath="/Volumes/d1/adugger/RS/MODIS_ARC/PROCESSED/FRNTRNG_LAI_SMOOTHED",
#' outputAs="one", removeOutlier=TRUE, outlierThreshold=0.5, lambda=1000,
#' overwrite=TRUE)
#' }
#' @keywords smooth
#' @concept MODIS dataAnalysis
#' @family MODIS
#' @export
SmoothStack <- function(inStack, w=NULL, t=NULL, lambda = 5000, nIter= 3,
outputAs="one", collapse=FALSE, outDirPath = "./",
removeOutlier=FALSE, outlierThreshold=NULL,
mergeDoyFun="max", ...) {
if (!file.exists(outDirPath)) {
dir.create(outDirPath, showWarnings = FALSE)
}
timeInfo <- MODIS::orgTime(inStack, pos1 = 3, pos2 = 13, format = "%Y.%m.%d",
pillow=0)
resultList <- MODIS::whittaker.raster(vi=inStack, w=w, t=t,
timeInfo=timeInfo, lambda=lambda,
nIter=nIter, outputAs=outputAs,
collapse=collapse,
outDirPath=outDirPath,
removeOutlier=removeOutlier,
outlierThreshold=outlierThreshold,
mergeDoyFun=mergeDoyFun, ...)
resultBrick <- resultList[[1]]
names(resultBrick) <- names(inStack)
resultBrick
}
#' Inserts pre-processed images into appropriate forcing NetCDF files by date.
#'
#' \code{InsertRS} takes a raster stack or brick of RS images and exports
#' individual images to matching (by date) forcing NetCDF files.
#'
#' \code{InsertRS} takes a raster stack or brick (as created by
#' \code{\link{ConvertRS2Stack}} or \code{\link{SmoothStack}}) or a NetCDF file
#' (as created by \code{\link{ConvertStack2NC}}) of RS images and exports each
#' layer (time step) to the appropriate time step forcing file. Only looks for
#' the date (not time) and inserts at the 00:00 hour on that date. The input
#' stack, brick, or file should be already processed through the
#' \code{\link{ConvertRS2Stack}}, \code{\link{SmoothStack}}, or
#' \code{\link{ConvertStack2NC}} tools or follow the same layer (date) naming
#' convention.
#'
#' @param inFile The name of the raster stack/brick or NetCDF file (full
#' pathname) to export.
#' @param forcPath Path to the forcing data you want to modify. Forcing data
#' files MUST match the size/resolution of the images in the inFile.
#' @param forcName The suffix for the forcing data files to modify
#' (DEFAULT="LDASIN_DOMAIN1")
#' @param varName Name for the NetCDF variable to export. The varibale will be
#' copied as-is, so make sure it matches the variable name needed in the
#' forcing data.
#' @param varUnit Units for the NetCDF export variable. Only required if the
#' inFile is a raster stack/brick. If the inFile is a NetCDF file, the units
#' will carry over.
#' @param varLong Long name for the NetCDF export variable. Only required if the
#' inFile is a raster stack/brick. If the inFile is a NetCDF file, the
#' longname will carry over.
#' @param varNA Value to set for "NA" or "no data". Default is -1.e+36.
#' @param overwrite Boolean to allow the tool to overwrite existing variables if
#' found in the forcing data. (DEFAULT=TRUE)
#' @return NULL
#'
#' @examples
#' ## Export the raster stack of LAI images created through ConvertRS2Stack to
#' ## the forcing data.
#'
#' \dontrun{
#' InsertRS(lai.b, forcPath="FORCING", forcName="LDASIN_DOMAIN3",
#' varName="LAI", varUnit="(m^2)/(m^2)", varLong="Leaf area index")
#'
#' ## Export the NetCDF of LAI images created through ConvertStack2NC to the
#' ## forcing data.
#'
#' InsertRS("BCNED_LAI.nc", forcPath="FORCING", forcName="LDASIN_DOMAIN3",
#' varName="LAI")
#' }
#' @keywords IO
#' @concept MODIS dataMgmt
#' @family MODIS
#' @export
InsertRS <- function(inFile, forcPath, forcName="LDASIN_DOMAIN1",
varName=NULL, varUnit=NULL, varLong=NULL, varNA=-1.e+36,
overwrite=TRUE) {
# Raster stack/brick case
if (inherits(inFile, "RasterStack") | inherits(inFile, "RasterBrick")) {
dtNames <- names(inFile)
for (i in dtNames) {
dtStr <- sub("DT", i, replacement="")
dtStrForc <- paste0(unlist(strsplit(dtStr,"[.]"))[1],
unlist(strsplit(dtStr,"[.]"))[2],
unlist(strsplit(dtStr,"[.]"))[3], "00")
ncFile <- ncdf4::nc_open(paste0(forcPath,"/",dtStrForc,".",forcName),
write=TRUE)
dimT <- ncdf4::ncdim_def( "Time", "", 1, unlim=TRUE, create_dimvar=T)
dimY <- ncdf4::ncdim_def( "south_north", "", 1:dim(inFile)[1], create_dimvar=T)
dimX <- ncdf4::ncdim_def( "west_east", "", 1:dim(inFile)[2], create_dimvar=T)
# NOTE: ncdf4 reads dimensions in reverse order from ncdump!
varNew <- ncdf4::ncvar_def(name=varName, units=varUnit,
dim=list(dimX, dimY, dimT),
missval=varNA, longname=varLong)
if ( (varName %in% names(ncFile$var)) ) {
if (overwrite) {
ncdf4::nc_close(ncFile)
system(paste0('ncks -O -x -v ', varName, ' ',
paste0(forcPath,"/",dtStrForc,".",forcName), ' ',
paste0(forcPath,"/",dtStrForc,".",forcName)))
ncFile <- ncdf4::nc_open(paste0(forcPath,"/",dtStrForc,".",forcName),
write=TRUE)
} else {
stop(paste0('Error: Variable ', varName, ' exists but overwite is set to FALSE. Exiting.'))
}
}
ncdf4::ncvar_add(ncFile, varNew)
# Close and re-open, otherwise ncdf4 throws error
ncdf4::nc_close(ncFile)
ncFile <- ncdf4::nc_open(paste0(forcPath,"/",dtStrForc,".",forcName),
write=TRUE)
ncdf4::ncvar_put(ncFile, varNew, RotateCw(raster::as.matrix(inFile[[i]])))
ncdf4::nc_close(ncFile)
}
# NetCDF file case
} else if ( inherits(inFile, "character") & file.exists(inFile) ) {
inNC <- ncdf4::nc_open(inFile)
nTime <- inNC$var[[varName]]$dim[[3]]$len
for (i in 1:(nTime)) {
dtNum <- inNC$var[[varName]]$dim[[3]]$val[[i]]
dtStr <- as.POSIXlt("198001", format = "%Y%j", tz = "UTC")
dtStr$mday <- dtStr$mday + dtNum
dtStrForc <- paste0(format(dtStr, "%Y%m%d"), "00")
if (overwrite) {
system(paste0("ncks -A -v ", varName, " -d Time,", i-1, ",", i-1, " ",
inFile , " ", paste0(forcPath,"/",dtStrForc,".",forcName)))
} else {
stop('Error: The NetCDF insert option uses ncks append which does NOT check for duplicate variables. Existing variables will automatically be overwritten. To OK this option, please set overwrite=TRUE.')
}
}
ncdf4::nc_open(inFile)
}
}
#' Calculates summary statistics from a remote sensing time series
#'
#' \code{CalcStatsRS} takes a raster stack/brick of RS images and generates a
#' time series dataframe of summary statistics across the domain for each time
#' step.
#'
#' \code{CalcStatsRS} takes a raster stack or brick of remote sensing images
#' over a time period (as created by \code{\link{ConvertRS2Stack}} or
#' \code{\link{SmoothStack}}) and generates a dataframe object that summarizes
#' all cells in the RS image at each time step. The statistics calculated are
#' mean, min, max, and standard deviation. This tool is useful for evaluating
#' how a smoothing function is impacting the time series of images.
#'
#' @param inStack The name of the raster stack or brick to calculate statistics
#' on.
#' @return A dataframe of statistics by time period (date).
#'
#' @examples
#' ## Calculate domain statistics and plot the mean.
#'
#' \dontrun{
#' stats.lai.b <- CalcStatsRS(lai.b)
#' with(stats.lai.b, plot(POSIXct, mean, typ='l'))
#' }
#' @keywords univar
#' @concept MODIS dataAnalysis
#' @family MODIS
#' @export
CalcStatsRS <- function(inStack) {
minDf <- as.data.frame(raster::cellStats(inStack, stat="min", na.rm=TRUE))
maxDf <- as.data.frame(raster::cellStats(inStack, stat="max", na.rm=TRUE))
meanDf <- as.data.frame(raster::cellStats(inStack, stat="mean", na.rm=TRUE))
sdDf <- as.data.frame(raster::cellStats(inStack, stat="sd", na.rm=TRUE))
statDf <- cbind(meanDf, minDf, maxDf, sdDf)
statDf$POSIXct <- as.POSIXct("1980-01-01", format="%Y-%m-%d", tz="UTC")
for (i in 1:nrow(statDf)) {
statDf$POSIXct[i] <- as.POSIXct(sub("DT",row.names(statDf)[i], replacement=""),
format="%Y.%m.%d", tz="UTC")
}
colnames(statDf) <- c("mean","min","max","sd","POSIXct")
rownames(statDf) <- NULL
statDf
}
#' Gap fill a series of MODIS images using interpolation.
#'
#' \code{GapFillRS} takes a set of pre-processed MODIS TIF files and
#' fills no-data gaps using interpolation.
#'
#' \code{GapFillRS} scans the specified directory and runs the GDAL python
#' script gdal_fillnodata.py. If a mask file is provided, GapFillRS will
#' convert specified mask cells to specified values AFTER the gap-filling.
#'
#' @param inPath The path to the directory that holds the already processed
#' MODIS TIF files.
#' @param matchStr The regular expression for filename matching (e.g.,
#' "*Lai_1km.tif").
#' @param outDir Pathname to the directory to store the output files.
#' @param begin Date string for the start date to include. The date string
#' should follow \href{http://r-forge.r-project.org/projects/modis/}{MODIS}
#' package convention (e.g., "2011.06.01"). (DEFAULT=NULL, all files are
#' processed).
#' @param end Date string for the end date to include. The date string should
#' follow \href{http://r-forge.r-project.org/projects/modis/}{MODIS} package
#' convention (e.g., "2011.06.01"). (DEFAULT=NULL, all files are processed).
#' @param maskFile Mask grid file. File must be in a GDAL-friendly format and
#' should match the dimensions and resolution of the inputs.
#' @param maskVals List of value sets to pull from mask file and convert.
#' The list should be in the form: list("maskValue"=outputValue),
#' e.g., list("1"=0, "16"=0) to convert all cells in the mask grid with mask
#' values of 1 and 16 to 0 in the output files.
#' @param nodataVal Value to use to represent no data. Make sure this value
#' is NOT in your mask raster. (DEFAULT=-9999)
#' @param scriptPath OPTIONAL path to where gdal_calc.py script lives in case
#' it is not in the same default path as gdal executables
#' @param pos1 From MODIS orgTime: Start position of date in the filename (DEFAULT=10).
#' @param pos2 From MODIS orgTime: End position of date in the filename (DEFAULT=16).
#' @param format From MODIS orgTime: How is the date formatted in the file, default
#' expects: YYYYDDD (\%Y\%j). Read ?as.Date for for more information.
#' @return null
#'
#' @examples
#' ## Gap fill a set of MODIS LAI images. We will convert all non-vegetation
#' ## land use types into 0.
#'
#' \dontrun{
#' maskVals <- list("1"=0, "16"=0, "19"=0, "23"=0,
#' "24"=0, "25"=0, "26"=0, "27"=0)
#' GapFillRS("~/wrfHydroTestCases/MODIS_ARC/PROCESSED/FOURMILE_LAI",
#' outDir="~/wrfHydroTestCases/MODIS_ARC/PROCESSED/FOURMILE_LAI_GAPFILL",
#' matchStr="*Lai_1km.tif",
#' begin="2013.05.01", end="2013.07.31",
#' maskFile="~/wrfHydroTestCases/Fourmile_Creek/DOMAIN/geo_LUINDEX.tif",
#' maskVals=maskVals)
#' }
#' @keywords IO
#' @concept MODIS dataMgmt
#' @family MODIS
#' @export
GapFillRS <- function(inPath, matchStr, outDir,
begin=NULL, end=NULL,
maskFile=NULL, maskVals=NULL,
nodataVal=-9999, scriptPath=NULL,
pos1=10, pos2=16, format="%Y%j") {
# Setup
if (!is.null(scriptPath)) scriptPath <- paste0(scriptPath, "/")
inPath <- path.expand(inPath)
outDir <- path.expand(outDir)
if (!is.null(maskFile)) maskFile <- path.expand(maskFile)
if (!file.exists(outDir)) {
dir.create(outDir, showWarnings = FALSE)
}
# Get file list
if (!is.null(begin) & !is.null(end)) {
beginDt <- format(as.POSIXct(begin, format="%Y.%m.%d", tz="UTC"), "%Y%j", tz="UTC")
endDt <- format(as.POSIXct(end, format="%Y.%m.%d", tz="UTC"), "%Y%j", tz="UTC")
timeInfo <- MODIS::orgTime(MODIS::preStack(path=inPath, pattern=matchStr),
begin=beginDt, end=endDt, pillow=0,
pos1=pos1, pos2=pos2, format=format)
} else {
timeInfo <- MODIS::orgTime(MODIS::preStack(path=inPath, pattern=matchStr), pillow=0,
pos1=pos1, pos2=pos2, format=format)
}
rsFilelist <- MODIS::preStack(path=inPath, pattern=matchStr, timeInfo=timeInfo)
#rsFilelist <- list.files(path=inPath, pattern=glob2rx(matchStr), full.names=TRUE)
# Mask case
if (!is.null(maskFile)) {
# Setup mask calc string
calcStr1 <- ''
calcStr2 <- ''
for (i in 1:length(maskVals)) {
maskIn <- names(maskVals)[i]
maskOut <- maskVals[[maskIn]]
if (calcStr1=='') {
calcStr1 <- paste0(maskOut, "*(A==", maskIn, ")")
calcStr2 <- paste0("+(", nodataVal, ")*((A!=", maskIn)
} else {
calcStr1 <- paste0(calcStr1, "+", maskOut, "*(A==", maskIn, ")")
calcStr2 <- paste0(calcStr2, ")&(A!=", maskIn)
}
}
calcStr <- paste0(calcStr1, calcStr2, "))")
# Create mask
tmpMask <- tempfile(fileext=".tif")
cmd_calc <- paste0(scriptPath,"gdal_calc.py ",
"-A ", "'", maskFile, "'",
" --outfile=", "'", tmpMask, "'",
" --calc='", calcStr,
"' --NoDataValue=", as.character(nodataVal),
" --overwrite")
system(cmd_calc)
}
# Loop through files
for (rsFile in rsFilelist) {
cmd_fill <- paste0(scriptPath,"gdal_fillnodata.py ",
rsFile, " ", outDir, "/", basename(rsFile))
system(cmd_fill)
# Mask case
if (!is.null(maskFile)) {
cmd_merge <- paste0(scriptPath,"gdal_merge.py ",
" -o ", "'", paste0(outDir, "/", basename(rsFile)), "'",
" -n ", as.character(nodataVal),
" '", paste0(outDir, "/", basename(rsFile)),
"' '", tmpMask, "'")
system(cmd_merge)
}
}
file.remove(tmpMask)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.