WSL/SLF GitLab Repository

Commit a1eeee4b authored by Adrien Michel's avatar Adrien Michel
Browse files

Implemented spatial correlation in R

parent 6c792c82
......@@ -29,6 +29,7 @@ export(plt.snow.days.maps)
export(plt.snow.days.maps.calib)
export(plt.snow.season.length)
export(plt.snow.season.length.calib)
export(plt.spatial.correlation.calib)
export(read.asc)
export(read.mask.asc)
export(read.mask.tif)
......
......@@ -63,8 +63,8 @@ snowSeasonLengthCpp <- function(data, indices, ncores, swe_threshold = 0.005, ma
#' everything else as FALSE.
#' @param x An integer vector
#' @export
spatialCorrelationCpp <- function(x_coords, y_coords, ncores) {
.Call(`_SnowQM_spatialCorrelationCpp`, x_coords, y_coords, ncores)
spatialCorrelationCpp <- function(x_coords, y_coords, ncores, dist_threshold) {
.Call(`_SnowQM_spatialCorrelationCpp`, x_coords, y_coords, ncores, dist_threshold)
}
#' xxx
......
......@@ -443,8 +443,8 @@ plt.seasonnal.summary.calib <- function(seasonnal.summary) {
for(season in seasons) {
for(variant in c("training","corrected")) {
text(1.005+i,-lim,
paste("Mean\n", round(out[["mean"]][[phase]][[variant]][[season]], 3), "\n",
"MAE\n", round(out[["mae"]][[phase]][[variant]][[season]], 3)),
paste("Mean\n", round(out[[phase]][[variant]][[season]][["mean"]], 3), "\n",
"MAE\n", round(out[[phase]][[variant]][[season]][["mae"]], 3)),
cex=0.5)
i=i+1
}
......
......@@ -2,6 +2,15 @@
############################# INNER FUNCTIONS ##################################
################################################################################
get.correlation <- function(data,indices, ncores) {
data.tmp <- rowMeanCpp(data$data, ncores)
diff=rep(0,times=length(data$x.indices))
for (pixel in 1:length(data$x.indices)) {
diff[pixel] = data.tmp[pixel]/mean(data.tmp[indices[[pixel]]])
}
return(matrix(diff))
}
################################################################################
############################# SINGLE MODEL RUN #################################
################################################################################
......@@ -28,26 +37,16 @@
#' @param swe.threshold The threshold SWE value to define snow days
#'
#' @author Adrien Michel, 2021 (WSL/SLF)
plt.spatial.correlation <- function(data.target.calib, data.training.calib,
data.corrected.calib, data.target.valid,
data.training.valid, data.corrected.valid,
out.plot.dir, ncores = 1, dist.threshold = 5) {
#' @export
plt.spatial.correlation.calib <- function(data.target.calib, data.training.calib,
data.corrected.calib, data.target.valid,
data.training.valid, data.corrected.valid,
out.plot.dir, ncores = 1, dist.threshold = 5) {
x.coords.vect <- data.target.calib$x.coords[data.target.calib$x.indices]/1000
y.coords.vect <- data.target.calib$y.coords[data.target.calib$y.indices]/1000
time <- Sys.time()
indices = spatialCorrelationCpp(x.coords.vect, y.coords.vect, ncores, dist.threshold)
print(Sys.time()-time)
get.correlation <- function(data,indices, ncores) {
data.tmp <- rowMeanCpp(data$data, ncores)
diff=rep(0,times=length(data$x.indices))
for (pixel in 1:length(data$x.indices)) {
diff[pixel] = data.tmp[pixel]/mean(data.tmp[indices[[pixel]]])
}
return(matrix(diff))
}
data.target.calib$correlation=get.correlation(data.target.calib, indices, ncores)
data.target.valid$correlation=get.correlation(data.target.valid, indices, ncores)
......@@ -73,18 +72,30 @@ plt.spatial.correlation <- function(data.target.calib, data.training.calib,
out = list()
par(mfrow=c(2,2))
out$training$calib$me <- round(mean(data.training.calib$correlation.rel),4)
out$training$calib$mae <- round(mean(abs(data.training.calib$correlation.rel)),3)
plot(data.training.calib,var="correlation.rel", main=paste("Traning data calib", me, mae), zlim=c(-1,1))
plot(data.training.calib,var="correlation.rel",
main=paste("Traning data calib", out$training$calib$me, out$training$calib$mae),
zlim=c(-1,1))
out$corrected$calib$me <- round(mean(data.corrected.calib$correlation.rel),4)
out$corrected$calib$mae <- round(mean(abs(data.corrected.calib$correlation.rel)),3)
plot(data.corrected.calib,var="correlation.rel", main=paste("Corrected data calib", me, mae),zlim=c(-1,1))
plot(data.corrected.calib,var="correlation.rel",
main=paste("Corrected data calib", out$corrected$calib$me, out$corrected$calib$mae),
zlim=c(-1,1))
out$training$valid$me <- round(mean(data.training.valid$correlation.rel),4)
out$training$valid$mae <- round(mean(abs(data.training.valid$correlation.rel)),3)
plot(data.training.valid,var="correlation.rel", main=paste("Target data valid", me, mae),zlim=c(-1,1))
plot(data.training.valid,var="correlation.rel",
main=paste("Target data valid", out$training$valid$me, out$training$valid$mae),
zlim=c(-1,1))
out$corrected$valid$me <- round(mean(data.corrected.valid$correlation.rel),4)
out$corrected$valid$mae <- round(mean(abs(data.corrected.valid$correlation.rel)),3)
plot(data.corrected.valid,var="correlation.rel", main=paste("Corrected data valid", me, mae),zlim=c(-1,1))
plot(data.corrected.valid,var="correlation.rel",
main=paste("Corrected data valid", out$corrected$valid$me, out$corrected$valid$mae),
zlim=c(-1,1))
dev.off()
return(out)
}
\ No newline at end of file
......@@ -255,18 +255,19 @@ assessment.plot.calib.internal <- function (data.target.calib,
as.numeric(round(difftime(Sys.time(),start.time.2, units = "mins"),2)),
"minutes"))
# print("[I] Plotting spatial correlation")
# start.time.2=Sys.time()
# plt.spatial.correlation(data.target.calib = data.target.calib,
# data.training.calib = data.training.calib,
# data.corrected.calib = data.corrected.calib,
# data.target.valid = data.target.valid,
# data.training.valid = data.training.valid,
# data.corrected.valid = data.corrected.valid,
# out.plot.dir = out.plot.dir, ncores = ncores, dist.threshold = dist.threshold)
# print(paste("[I] Snow spatial correlation plot produced in",
# as.numeric(round(difftime(Sys.time(),start.time.2, units = "mins"),2)),
# "minutes"))
print("[I] Plotting spatial correlation")
start.time.2=Sys.time()
plt.spatial.correlation.calib(data.target.calib = data.target.calib,
data.training.calib = data.training.calib,
data.corrected.calib = data.corrected.calib,
data.target.valid = data.target.valid,
data.training.valid = data.training.valid,
data.corrected.valid = data.corrected.valid,
out.plot.dir = out.plot.dir, ncores = ncores,
dist.threshold = dist.threshold)
print(paste("[I] Snow spatial correlation plot produced in",
as.numeric(round(difftime(Sys.time(),start.time.2, units = "mins"),2)),
"minutes"))
#Clean maps data from output
snow.season$data.target.calib <- NULL
......
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/assessment_spatial_correlation.R
\name{plt.spatial.correlation}
\alias{plt.spatial.correlation}
\name{plt.spatial.correlation.calib}
\alias{plt.spatial.correlation.calib}
\title{XXXXX}
\usage{
plt.spatial.correlation(
plt.spatial.correlation.calib(
data.target.calib,
data.training.calib,
data.corrected.calib,
......
......@@ -4,7 +4,7 @@
\alias{spatialCorrelationCpp}
\title{xxx}
\usage{
spatialCorrelationCpp(x_coords, y_coords, ncores)
spatialCorrelationCpp(x_coords, y_coords, ncores, dist_threshold)
}
\arguments{
\item{x}{An integer vector}
......
......@@ -89,15 +89,16 @@ BEGIN_RCPP
END_RCPP
}
// spatialCorrelationCpp
Rcpp::List spatialCorrelationCpp(const Rcpp::NumericVector x_coords, const Rcpp::NumericVector y_coords, const size_t ncores);
RcppExport SEXP _SnowQM_spatialCorrelationCpp(SEXP x_coordsSEXP, SEXP y_coordsSEXP, SEXP ncoresSEXP) {
Rcpp::List spatialCorrelationCpp(const Rcpp::NumericVector x_coords, const Rcpp::NumericVector y_coords, const size_t ncores, const double dist_threshold);
RcppExport SEXP _SnowQM_spatialCorrelationCpp(SEXP x_coordsSEXP, SEXP y_coordsSEXP, SEXP ncoresSEXP, SEXP dist_thresholdSEXP) {
BEGIN_RCPP
Rcpp::RObject rcpp_result_gen;
Rcpp::RNGScope rcpp_rngScope_gen;
Rcpp::traits::input_parameter< const Rcpp::NumericVector >::type x_coords(x_coordsSEXP);
Rcpp::traits::input_parameter< const Rcpp::NumericVector >::type y_coords(y_coordsSEXP);
Rcpp::traits::input_parameter< const size_t >::type ncores(ncoresSEXP);
rcpp_result_gen = Rcpp::wrap(spatialCorrelationCpp(x_coords, y_coords, ncores));
Rcpp::traits::input_parameter< const double >::type dist_threshold(dist_thresholdSEXP);
rcpp_result_gen = Rcpp::wrap(spatialCorrelationCpp(x_coords, y_coords, ncores, dist_threshold));
return rcpp_result_gen;
END_RCPP
}
......@@ -120,7 +121,7 @@ static const R_CallMethodDef CallEntries[] = {
{"_SnowQM_computeSeasonnalSummaryCpp", (DL_FUNC) &_SnowQM_computeSeasonnalSummaryCpp, 4},
{"_SnowQM_snowDaysComparisonCpp", (DL_FUNC) &_SnowQM_snowDaysComparisonCpp, 5},
{"_SnowQM_snowSeasonLengthCpp", (DL_FUNC) &_SnowQM_snowSeasonLengthCpp, 5},
{"_SnowQM_spatialCorrelationCpp", (DL_FUNC) &_SnowQM_spatialCorrelationCpp, 3},
{"_SnowQM_spatialCorrelationCpp", (DL_FUNC) &_SnowQM_spatialCorrelationCpp, 4},
{"_SnowQM_rowMeanCpp", (DL_FUNC) &_SnowQM_rowMeanCpp, 2},
{NULL, NULL, 0}
};
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment