knitr::opts_chunk$set(echo=T, comment=NA, error=T, warning=F, message = F, fig.align = 'center', results="hold")
z
at selected Ppr
and Tpr
# get a z value using HY library(zFactor) z.Shell(pres.pr = 1.5, temp.pr = 2.0)
From the Standing-Katz chart we obtain a digitized point:
# get a z value from the SK chart at the same Ppr and Tpr library(zFactor) tpr_vec <- c(2.0) getStandingKatzMatrix(tpr_vector = tpr_vec, pprRange = "lp")[1, "1.5"]
z
at selected Ppr
and Tpr
library(zFactor) z.Shell(pres.pr = 1.5, temp.pr = 1.1)
From the Standing-Katz chart we obtain a digitized point:
library(zFactor) tpr_vec <- c(1.1) getStandingKatzMatrix(tpr_vector = tpr_vec, pprRange = "lp")[1, "1.5"]
We perceive a noticeable difference between the values of z from the HY calculation and the value read from the Standing-Katz chart.
z
for several Ppr
and Tpr
In this example we provide vectors instead of a single point.
library(zFactor) ppr <- c(0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5) tpr <- c(1.3, 1.5, 1.7, 2) z.Shell(ppr, tpr)
Which is equivalent to using the sapply function with the internal function .z.HallYarborough
, which we call adding the prefix zFactor:::
. That is, the package name and three dots.
# test HY with 1st-derivative using the values from paper ppr <- c(0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5) tpr <- c(1.3, 1.5, 1.7, 2) hy <- sapply(ppr, function(x) sapply(tpr, function(y) zFactor:::.z.Shell(pres.pr = x, temp.pr = y))) rownames(hy) <- tpr colnames(hy) <- ppr print(hy)
With the same ppr and tpr vector, we do the same for the Standing-Katz chart:
library(zFactor) sk <- getStandingKatzMatrix(ppr_vector = ppr, tpr_vector = tpr) sk
Subtract and find the difference:
err <- round((sk - hy) / sk * 100, 2) err
print(colSums(err))
print(rowSums(err))
Tpr
library(zFactor) tpr2 <- c(1.05, 1.1) ppr2 <- c(0.5, 1.5, 2.5, 3.5, 4.5, 5.5) sk2 <- getStandingKatzMatrix(ppr_vector = ppr2, tpr_vector = tpr2, pprRange = "lp") sk2
We do the same with the correlation:
# calculate z values at lower values of Tpr library(zFactor) tpr <- c(1.05, 1.1) ppr <- c(0.5, 1.5, 2.5, 3.5, 4.5, 5.5) corr2 <- z.Shell(pres.pr = ppr, temp.pr = tpr) print(corr2)
err2 <- round((sk2 - corr2) / sk2 * 100, 2) err2
We can see that using Hall-Yarborough correlation shows a very high error at values of Tpr
lower or equal than 1.1
being Tpr=1.05
the worst curve to calculate z values from.
t_err2 <- t(err2) t_err2
Applying the function summary
:
sum_t_err2 <- summary(t_err2) sum_t_err2
We can see that the errors in z
are considerable with a r sum_t_err2[1,1]
% and r sum_t_err2[6,1]
% for Tpr=1.05
, and a r sum_t_err2[1,2]
% and r sum_t_err2[6,2]
% for Tpr=1.10
SK
chart values vs HY
correlationlibrary(zFactor) library(tibble) tpr2 <- c(1.05, 1.1, 1.2, 1.3) ppr2 <- c(0.5, 1.0, 1.5, 2, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0, 6.5) sk_corr_2 <- createTidyFromMatrix(ppr2, tpr2, correlation = "SH") as.tibble(sk_corr_2)
Plotting the difference between the z values in the Standing-Katz and the values calculated by Hall-Yarborough:
library(ggplot2) p <- ggplot(sk_corr_2, aes(x=Ppr, y=z.calc, group=Tpr, color=Tpr)) + geom_line() + geom_point() + geom_errorbar(aes(ymin=z.calc-dif, ymax=z.calc+dif), width=.4, position=position_dodge(0.05)) print(p)
Tpr
curveslibrary(zFactor) library(ggplot2) library(tibble) # get all `lp` Tpr curves tpr_all <- getCurvesDigitized(pprRange = "lp") ppr <- c(0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5) sk_corr_all <- createTidyFromMatrix(ppr, tpr_all, correlation = "SH") as.tibble(sk_corr_all) p <- ggplot(sk_corr_all, aes(x=Ppr, y=z.calc, group=Tpr, color=Tpr)) + geom_line() + geom_point() + geom_errorbar(aes(ymin=z.calc-dif, ymax=z.calc+dif), width=.4, position=position_dodge(0.05)) print(p)
The greatest errors are localized in two of the Tpr curves: at 1.05 and 1.1
# MSE: Mean Squared Error # RMSE: Root Mean Sqyared Error # RSS: residual sum of square # ARE: Average Relative Error, % # AARE: Average Absolute Relative Error, % library(dplyr) grouped <- group_by(sk_corr_all, Tpr, Ppr) smry_tpr_ppr <- summarise(grouped, RMSE= sqrt(mean((z.chart-z.calc)^2)), MSE = sum((z.calc - z.chart)^2) / n(), RSS = sum((z.calc - z.chart)^2), ARE = sum((z.calc - z.chart) / z.chart) * 100 / n(), AARE = sum( abs((z.calc - z.chart) / z.chart)) * 100 / n() ) ggplot(smry_tpr_ppr, aes(Ppr, Tpr)) + geom_tile(data=smry_tpr_ppr, aes(fill=AARE), color="white") + scale_fill_gradient2(low="blue", high="red", mid="yellow", na.value = "pink", midpoint=12.5, limit=c(0, 25), name="AARE") + theme(axis.text.x = element_text(angle=45, vjust=1, size=11, hjust=1)) + coord_equal() + ggtitle("Shell", subtitle = "SH")
# get all `lp` Tpr curves tpr <- getCurvesDigitized(pprRange = "lp") ppr <- c(0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5) # calculate HY for the given Tpr all_corr <- sapply(ppr, function(x) sapply(tpr, function(y) z.Shell(pres.pr = x, temp.pr = y))) rownames(all_corr) <- tpr colnames(all_corr) <- ppr cat("Calculated correlation\n") print(all_corr) cat("\nStanding-Katz chart\n") all_sk <- getStandingKatzMatrix(ppr_vector = ppr, tpr_vector = tpr) all_sk # find the error cat("\n Errors in percentage \n") all_err <- round((all_sk - all_corr) / all_sk * 100, 2) # in percentage all_err cat("\n Errors in Ppr\n") summary(all_err) # for the transposed matrix cat("\n Errors for the transposed matrix: Tpr \n") summary(t(all_err))
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.