Background

We are using WRF-Hydro to predict streamflow for Fourmile Creek at the Orodell USGS gage for the 2013 snowmelt period. We want to obtain MODIS LAI images for our region of interest and insert the processed images into our forcing time series so we can use dynamic LAI in future WRF-Hydro model runs.

Load the rwrfhydro package.

library("rwrfhydro")
## To check rwrfhydro updates run: CheckForUpdates()

Set a data path to the Fourmile Creek test case.

dataPath <- '~/wrfHydroTestCases/Fourmile_Creek'

Setup the MODIS-R tool.

NOTE: These tools rely heavily on the R-Forge MODIS-R package, so please see their development page and R help pages for more on the package tools and options. Currently, we are making some edits to the MODIS tool to add some capacity specific to WRF-Hydro applications. So please install the rwrfhydro version of MODIS from github (devtools::install_github("aubreyd/modis4rwrfhydro")). A local system GDAL install is also required.

To start, we need to load a few libraries and then tell MODIS-R where our working directories are. If this is the first time using MODIS-R, we will also need to tell it where the local GDAL libraries live. Your GDAL pathname may vary; type “which gdalinfo” in a terminal window to find the appropriate path for your system. If gdalinfo is not found, make sure you have GDAL installed on your local system.

library(rgeos)
## rgeos version: 0.3-8, (SVN revision (unknown))
##  GEOS runtime version: 3.4.2-CAPI-1.8.2 r3921 
##  Polygon checking: TRUE
library(rgdal)
## Loading required package: sp
## rgdal: version: 0.9-1, (SVN revision 518)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 1.11.0, released 2014/04/16
## Path to GDAL shared files: /Library/Frameworks/GDAL.framework/Versions/1.11/Resources/gdal
## Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480]
## Path to PROJ.4 shared files: (autodetected)
library(raster)
library(MODIS)
## MODIS_manual: https://ivfl-rio.boku.ac.at/owncloud/public.php?service=files&t=660dc830afb091237cc40b3dea2fdf6b
library(ncdf4)

## if gdalPath returns 1=error, then try these two lines, substiting your shell profile script as appropriate.
# bashPath <- system('. ~/.bash_profile; echo $PATH', intern = TRUE)
# Sys.setenv(PATH=bashPath)
gdalPath <- sub("gdalinfo", replacement="", system("which gdalinfo", intern=TRUE))
MODISoptions(localArcPath="~/wrfHydroTestCases/MODIS_ARC", 
             outDirPath="~/wrfHydroTestCases/MODIS_ARC/PROCESSED", gdalPath=gdalPath)
## To install all required and suggested packages run:
##  setRepositories() # activate CRAN, R-forge, and Omegahat and then: 
##  install.packages(c(' mapdata '),dependencies=TRUE)
## Setting 'localArcPath' to '/Volumes/d1/adugger/WRF_Hydro/test_cases/MODIS_ARC'
## If you already have downloaded some HDF-files to '/Volumes/d1/adugger/WRF_Hydro/test_cases/MODIS_ARC' you can use '?orgStruc()' to re-arrange your HDF-data!
## 'outDirPath' has been changed from '/Volumes/d1/adugger/WRF_Hydro/test_cases/MODIS_ARC/PROCESSED' to '/Volumes/d1/adugger/WRF_Hydro/test_cases/MODIS_ARC/PROCESSED'
##   'MRT_HOME' not set/found! MRT is NOT enabled! See: 'https://lpdaac.usgs.gov/tools/modis_reprojection_tool'
## 
## STORAGE:
## _______________
## localArcPath : /Volumes/d1/adugger/WRF_Hydro/test_cases/MODIS_ARC 
## outDirPath   : /Volumes/d1/adugger/WRF_Hydro/test_cases/MODIS_ARC/PROCESSED 
## 
## 
## DOWNLOAD:
## _______________
## MODISserverOrder : LPDAAC, LAADS 
## dlmethod         : auto 
## stubbornness     : high 
## 
## 
## PROCESSING:
## _______________
## GDAL           : GDAL 1.11.0, released 2014/04/16 
## MRT            : Not available. Use 'MODIS:::checkTools('MRT')' for more information! 
## pixelSize      : asIn 
## outProj        : asIn 
## resamplingType : NN 
## dataFormat     : GTiff 
## 
## 
## DEPENDENCIES:
## _______________
## 
## NULL

Download and process MODIS tiles

Here we download the relevant MODIS tiles and resample to the Fourmile Creek geogrid. The GetMODIS tool calls the MODIS-R runGdal tool to identify the required MODIS tiles to cover the domain, download the appropriate product tiles, mosaic them, and clip/resample (nearest neighbor) to the geogrid extent and resolution. The result will be a directory of processed TIF images in the outDirPath specified in MODISoptions. We specify the geogrid, a start and end date, and the MOD15A2 (FPAR/LAI) MODIS product name.

GetMODIS(geogrdPath=paste0(dataPath, '/DOMAIN/geo_em.d01_1km.nc'), prodName="MOD15A2", 
         outDir="FOURMILE_LAI", begin="2013.03.01", end="2013.07.31")
## ########################
## outProj          =  +proj=lcc +lat_1=30 +lat_2=50 +lat_0=39.9400100708008 +lon_0=-105 +x_0=0 +y_0=0 +a=6370000 +b=6370000 +units=m +no_defs  (Specified by raster*/spatial* object)
## pixelSize        =  1000 1000  (Specified by raster* object)
## resamplingType   =  near 
## Output Directory =  /Volumes/d1/adugger/WRF_Hydro/test_cases/MODIS_ARC/PROCESSED/FOURMILE_LAI 
## ########################

Create a raster stack of the downloaded and processed MODIS LAI images. This tool does some additional processing of the images to get them into a more usable format. Specifically, it removes “no data” values (in this case, any value over 100) and applies product-specific scale factors (in this case, 0.1). See the MODIS data page for specifics on valid value ranges and scale factors.

lai.b <- ConvertRS2Stack(paste0(options("MODIS_outDirPath"), '/FOURMILE_LAI'), "*Lai_1km.tif$", 
                         begin="2013.03.01", end="2013.07.31", noData=100, 
                         noDataQual="max", valScale=0.1, valAdd=0)
## Found 19 files!
## Found 19 files!

See the help on ConvertRS2Stack for more details on the processing options.

?ConvertRS2Stack
ConvertRS2Stack R Documentation

Convert a set of MODIS images to a raster stack and, optionally, a NetCDF file.

Description

ConvertRS2Stack takes a set of pre-processed MODIS TIF files and creates a raster stack and, optionally, a NetCDF file.

Usage

ConvertRS2Stack(inPath, matchStr, begin = NULL, end = NULL, noData = NULL,
  noDataQual = "exact", valScale = 1, valAdd = 0, outFile = NULL,
  varName = NULL, varUnit = NULL, varLong = NULL, varNA = -1e+36)

Arguments

inPath

The path to the directory that holds the already processed MODIS TIF files.

matchStr

The regular expression for filename matching (e.g., “*Lai_1km.tif“).

begin

Date string for the start date to include. The date string should follow MODIS package convention (e.g., “2011.06.01”). (DEFAULT=NULL, all files are processed).

end

Date string for the end date to include. The date string should follow MODIS package convention (e.g., “2011.06.01”). (DEFAULT=NULL, all files are processed).

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)

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)

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)

valAdd

The addition value to apply to the image values per the MODIS product. The adjustment is applied as VALUE * valScale + valAdd. (DEFAULT = 0)

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.

varName

Name for the NetCDF export variable. Only required if outFile is provided.

varUnit

Units for the NetCDF export variable. Only required if outFile is provided.

varLong

Long name for the NetCDF export variable. Only required if outFile is provided.

varNA

Value to set for “NA” or “no data”. Default is -1.e+36. Only required if outFile is provided.

Details

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 GetMODIS tool or follow the same file naming convention used by the MODIS runGdal tool. See the R MODIS package for more on those specifications.

Value

A raster stack.

See Also

Other MODIS: CalcStatsRS; ConvertStack2NC; GetMODIS; InsertRS; SmoothStack

Examples

## Import the already processed LAI TIF images into a raster stack. Use
## the full time series of images.

## Not run: 
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=c("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")

## End(Not run)


We can export the processed raster stack to a NetCDF file for use outside of R.

ConvertStack2NC(lai.b, outFile=paste0(dataPath, '/ANALYSIS/FOURMILE_LAI_SUM13.nc'), 
                varName="LAI", varUnit="(m^2)/(m^2)", varLong="Leaf area index", varNA=-1.e+36)

Let’s take a look at the raster stack we created. List the available layers (labelled by date).

names(lai.b)
##  [1] "DT2013.03.06" "DT2013.03.14" "DT2013.03.22" "DT2013.03.30"
##  [5] "DT2013.04.07" "DT2013.04.15" "DT2013.04.23" "DT2013.05.01"
##  [9] "DT2013.05.09" "DT2013.05.17" "DT2013.05.25" "DT2013.06.02"
## [13] "DT2013.06.10" "DT2013.06.18" "DT2013.06.26" "DT2013.07.04"
## [17] "DT2013.07.12" "DT2013.07.20" "DT2013.07.28"

Plot a few of the images in the LAI raster stack. We can choose which layer to plot by index or name.

plot(lai.b, 1)

plot(lai.b, "DT2013.07.12")

Apply smoothing filters

Apply a whittaker filter to smooth the LAI time series. This is generally not necessary for a derived product like LAI, but may be useful for products like NDVI, snow cover fraction, vegetated fraction, etc.

Try a highly smoothed filter.

lai.b.sm1 <- SmoothStack(lai.b, 
                         outDirPath=paste0(options("MODIS_outDirPath"), '/FOURMILE_LAI_SMOOTH1'), 
                         outputAs="one", removeOutlier=TRUE, outlierThreshold=0.5, 
                         lambda=1000, overwrite=TRUE)
## Loading required package: ptw
## Yearly 'lambda' is: 1000 
## Now changed with lambda*('length of input data period in days'/365) to: 391.7808 
## Data is in, start processing!

Try a less aggressive filter.

lai.b.sm2 <- SmoothStack(lai.b, 
                         outDirPath=paste0(options("MODIS_outDirPath"), '/FOURMILE_LAI_SMOOTH2'), 
                         outputAs="one", lambda=20, overwrite=TRUE)
## Yearly 'lambda' is: 20 
## Now changed with lambda*('length of input data period in days'/365) to: 7.835616 
## Data is in, start processing!

Calculate statistics and plot the domain means over time.

stats.lai.b <- CalcStatsRS(lai.b)
stats.lai.b.sm1 <- CalcStatsRS(lai.b.sm1)
stats.lai.b.sm2 <- CalcStatsRS(lai.b.sm2)
head(stats.lai.b)
## Warning: package 'pander' was built under R version 3.1.2
mean min max sd POSIXct
0.5407 0.1 1.3 0.3119 2013-03-06
0.5386 0.1 1.2 0.2389 2013-03-14
0.465 0.1 0.9 0.2081 2013-03-22
0.6664 0 1.7 0.2591 2013-03-30
0.6886 0 1.3 0.2916 2013-04-07
0.2236 0 0.6 0.1084 2013-04-15
with(stats.lai.b, plot(POSIXct, mean, typ='l'))
with(stats.lai.b.sm1, lines(POSIXct, mean, col='red'))
with(stats.lai.b.sm2, lines(POSIXct, mean, col='blue'))
legend("topleft", c("MODIS Raw", "Smooth Filter 1", "Smooth Filter 2"), 
       col=c("black", "red", "blue"), lty=c(1,1,1))

Export to forcing

Export the second smoothed time series to the forcing data for use in future model runs.

InsertRS(lai.b.sm2, forcPath=paste0(dataPath, '/FORCING'), forcName="LDASIN_DOMAIN1", 
         varName="LAI", varUnit="(m^2)/(m^2)", varLong="Leaf area index", overwrite=TRUE)