xuanlongma

## This expression will list all files within the dir.hdf 
## directory which contains ‘h29v11’ AND ended with the ‘.hdf’ extension.
files.hdf.h29v11 <- list.files(dir.hdf, pattern = '*.h29v11.*.hdf$', 
                               recursive = TRUE, full.names = TRUE)

xuanlongma

## invoke system command to convert shapefiles to raster

## Please install gdal before using this snippet
## http://www.gdal.org

## 12000: the width of the output raster
## 6400: the height of the output raster
## layer: the attribute field
system(command = 'gdal_rasterize -a input -ts 12000 6400 -l 
       layer /Users/me/Desktop/input.shp /Users/me/Desktop/output.tif')

xuanlongma

## Convert shapefile to raster

## Packages
require(sp)
require(raster)

## Import shapefile
shp.x = readOGR('shapefile.shp', layer = 'shapefile')

## Create an empty raster

xuanlongma

## RMSE: Root Mean Squared Error

## obs: observations
## pred: predictions
fun.rmse <- function(obs, pred) {
  sqrt(mean((obs-pred)^2))
}

xuanlongma

## m is the nonlinear model, x is the variable

r2 <- 1 - (deviance(m) / sum((x - mean(x)) ^ 2))

xuanlongma

## Thanks RobertH for sharing this very useful function
## http://stackoverflow.com/a/9937083

fun.quickstack <- function(f) {
  r <- raster(f[1])
  ln <- extension(basename(f), '')
  s <- stack(r)
  s@layers <- sapply(1:length(f), function(x){ r@file@name = f[x]; 
  r@layernames=ln[x]; r@data@haveminmax=FALSE ; r })
  s@layernames <- ln

xuanlongma

## Original arthur: vikrant @ R help mailing list
## return the maximum CCF value along with corresponding lag value
## http://r.789695.n4.nabble.com/ccf-function-td2288257.html

fun.ccfmax <- function(a,b) { 
 d <- ccf(a, b, plot = FALSE) 
 cor = d$acf[,,1] 
 lag = d$lag[,,1] 
 res = data.frame(cor,lag) 
 res_max = res[which.max(res$cor),] 

xuanlongma

## Compute the cross-correlation Pearson product-moment correlation coefficient
## and check the p.value 
## This function is an alternative of the native ccf() function which does not 
## return p.value by default

fun.ccfmax <- function(x) {
 a <- x[1:(length(x) / 2)]
 b <- x[-(1:(length(x) / 2))]
 if (all(is.na(a)) | all(is.na(b))) {
   return(c(NA, NA, NA))

xuanlongma

## x: daily rainfall time series
## lindex: date stamps of the MOD13 VI product time series
## rindex = lindex + 15
## TODO: increase the efficiency by replace the loop use sapply

fun.aggre16d <- function(x, lindex, rindex) {
  v.16d.rain <- vector('numeric', length(lindex))
  for (i in 1:length(lindex)) {
    v.16d.rain[i] <- sum(x[lindex[i]:rindex[i]], na.rm = T)
  }

xuanlongma

## Arthur: Roman Luštrik
## http://stackoverflow.com/questions/9017070/set-the-color-in-plot-xts

plot.xts2 <- function (x, y = NULL, type = "l", auto.grid = TRUE, major.ticks = "auto", 
    minor.ticks = TRUE, major.format = TRUE, bar.col = "grey", 
    candle.col = "white", ann = TRUE, axes = TRUE, col = "black", ...) 
{
    series.title <- deparse(substitute(x))
    ep <- axTicksByTime(x, major.ticks, format = major.format)
    otype <- type