require(ncdf4)
require(raster)
require(data.table)
require(reshape2)
require(ggplot2)
transform.dt.to.brick = function(my.dt, variable.name) {
array.from.dt = acast(my.dt, lat ~ lon ~ time, value.var = variable.name)
my.brick = brick(array.from.dt, ymn = min(as.numeric(rownames(array.from.dt))),
ymx = max(as.numeric(rownames(array.from.dt))), xmn = min(as.numeric(colnames(array.from.dt))),
xmx = max(as.numeric(colnames(array.from.dt))))
my.brick = flip(my.brick, direction = "2")
return(my.brick)
}
####OWDA
test <- readRDS("../OWDA/Data/OWDA_1000.rds")
colnames(test) = c("time", "lat", "lon", "scPDSI")
owda_nc = transform.dt.to.brick(test, variable.name = "scPDSI")
owda_nc_1001 <- owda_nc$X1001
owda_nc_1005 <- owda_nc$X1005
owda_nc_rnorm <- owda_nc$X1005
owda_nc_rnorm[,] <- rnorm(length(owda_nc_rnorm[,]), mean = 1, sd = 1.5)
ncells <- length(owda_nc_1001)
owda_nc_1001_agg = list(owda_nc_1001)
owda_nc_1005_agg = list(owda_nc_1005)
owda_nc_rnorm_agg = list(owda_nc_rnorm)
i <- 2
while(ncells > 50){
owda_nc_1001_agg[[i]] <- aggregate(owda_nc_1001, fact = i)
owda_nc_1005_agg[[i]] <- aggregate(owda_nc_1005, fact = i)
owda_nc_rnorm_agg[[i]] <- aggregate(owda_nc_rnorm, fact = i)
ncells <- length(owda_nc_1001_agg[[i]])
ncells <- length(owda_nc_1005_agg[[i]])
ncells <- length(owda_nc_rnorm_agg[[i]])
i <- i + 1
print(i)
}
aa <- sapply(sapply(owda_nc_1001_agg, getValues), sd, na.rm = T)
bb <- sapply(sapply(owda_nc_1005_agg, getValues), sd, na.rm = T)
######NASA MERRA PRCP
test <- brick("merra_tp_-10-50E_30-70N.nc")
test_x7 <- test$X7
test_x7[,] <- scale(test_x7[,], center = T, scale = T)
test_x10 <- test$X10
test_x10[,] <- scale(test_x10[,], center = T, scale = T) #For monthly (mm) * 3600 * 24 * 30
ncells <- length(test_x7)
test_7_agg = list(test_x7)
test_10_agg = list(test_x7)
i <- 2
while(ncells > 50){
test_7_agg[[i]] <- aggregate(test_x7, fact = i)
test_10_agg[[i]] <- aggregate(test_x10, fact = i)
ncells <- length(test_7_agg[[i]])
ncells <- length(test_10_agg[[i]])
i <- i + 1
print(i)
}
aa <- sapply(sapply(test_7_agg, getValues), sd, na.rm = T)
bb <- sapply(sapply(test_10_agg, getValues), sd, na.rm = T)
aa_mat <- cbind(1:14, aa)
bb_mat <- cbind(1:14, bb)
scalegram_plot(aa_mat)
aabb <- rbind(cbind(aa_mat, 7), cbind(bb_mat, 10))
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.