Dear all,
I have a working code to compute trends on a NetCDF time series. However, due to the size of the dataset (global scale), I need assistance in running this code in parallel on an HPC system.
library(ncdf4)
library(trend)
library(Kendall)
library(doParallel)
library(foreach)
# Open the NetCDF file
prec_month_max <- nc_open("Prec_1981-2010.nc")
lon <- ncvar_get(prec_month_max, varid = "lon")
lat <- ncvar_get(prec_month_max, varid = "lat")
time <- ncvar_get(prec_month_max, varid = "time")
nx <- length(lon)
ny <- length(lat)
nt <- length(time)
obs <- ncvar_get(prec_month_max, "pr")
print("calc rl")
rl = array(0, c(nx, ny))
mk = array(0, c(nx, ny))
pv = array(0, c(nx, ny))
for (x in 1:nx) {
print(x)
for (y in 1:ny) {
tryCatch({
Rr = ts(obs[x, y, ], start = c(1981), frequency = 12)
Rr=Rr[-(1:6)]
slp = sens.slope(Rr)
rl[x, y] = slp$estimates
mk_res = MannKendall(Rr)
tau = mk_res$tau[1]
pvalue = mk_res$sl[1]
mk[x, y]=tau
pv[x, y]=pvalue
}, error=function(e){cat("ERROR :",conditionMessage(e), "\n")})
}
}
time <- ncdim_def("Time", "months", 1, unlim=TRUE)
lons <- ncdim_def("longitude", "degrees_east", lon)
lats <- ncdim_def("latitude", "degrees_north", lat)
var_temp_rl <- ncvar_def("sens_slope", "years", list(lons, lats, time), longname="Sen's slope")
var_temp_mk <- ncvar_def("mk", "years", list(lons, lats, time), longname="Modified MK trend")
var_temp_pv <- ncvar_def("pv", "years", list(lons, lats, time), longname="Pvalue")
ncnew <- nc_create("Global_MK_trend.nc", list(var_temp_rl, var_temp_mk, var_temp_pv))
ncvar_put(ncnew, var_temp_rl, rl, start = c(1, 1, 1), count = c(nx, ny, 1))
ncvar_put(ncnew, var_temp_mk, mk, start = c(1, 1, 1), count = c(nx, ny, 1))
ncvar_put(ncnew, var_temp_pv, pv, start = c(1, 1, 1), count = c(nx, ny, 1))
nc_close(ncnew)
Thank you for your kind support.
Best wishes,
Solomon