PPC - Analysis of CompuTree outputs
Roberto Scotti
NuoroForestrySchool19 novembre 2018
Abstract
TO BE COMPLETEDlibrary(tidyverse)## -- Attaching packages --------------------------------------- tidyverse 1.2.1 --## v ggplot2 3.0.0 v purrr 0.2.5
## v tibble 1.4.2 v dplyr 0.7.6
## v tidyr 0.8.1 v stringr 1.3.1
## v readr 1.1.1 v forcats 0.3.0## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()library(googlesheets)
suppressMessages(library(dplyr))
# gs_ls()
catalog <-
"https://docs.google.com/spreadsheets/d/1CBrzWp85K1Gj5zzljlVwzMUWvAgGd5_WP5vJY1Iq2A8/edit#gid=0" %>%
gs_url()## Sheet-identifying info appears to be a browser URL.
## googlesheets will attempt to extract sheet key from the URL.## Putative key: 1CBrzWp85K1Gj5zzljlVwzMUWvAgGd5_WP5vJY1Iq2A8## Auto-refreshing stale OAuth token.## Sheet successfully identified: "RilieviPiccolo SpRiViFo"catalog %>%
gs_ws_ls()## [1] "stazione" "FC" "Ore" "scaling mean"catalog_FC <- catalog %>% gs_read("FC")## Accessing worksheet titled 'FC'.## Parsed with column specification:
## cols(
## Id_foto = col_integer(),
## elab_AB = col_character(),
## elab_A = col_character(),
## Id_staz = col_integer(),
## Id_fusti = col_integer(),
## Id_dbhAB = col_integer(),
## Id_primo_cilindroAB = col_integer(),
## Id_dbhA = col_integer(),
## Id_primo_cilindroA = col_integer(),
## rotelle = col_character(),
## grad_area = col_character(),
## progressivo = col_character(),
## verso_prog = col_character(),
## DBH = col_integer(),
## specie = col_character(),
## Note = col_character()
## )# write.dcf(catalog_FC, "catalog_FC_offline.dcf")
# catalog_FC <- read.dcf("catalog_FC_offline.dcf")
# just to test data fetching!
catalog_FC[1:2,] %>%
select(Id_foto, Id_dbhAB)## # A tibble: 2 x 2
## Id_foto Id_dbhAB
## <int> <int>
## 1 1 1236559
## 2 2 1234171data_path <- "G:\\Il mio Drive\\DocumentiGDuniss\\02-RICERCA\\2017_Arzana-pineta\\RilieviPiccolo\\"
ds <- c("A", "AB") # disegno sperimentale
ps <- c("dbh", "cyl") # processing strategy
catalog_csv <- NULL
for(d in ds) { # A o AB
p <- paste0(data_path, "ds", d)
for(ft in ps) { # dbh or cyl
print(paste("design:", d, " - strategy:", ft))
fl <- dir(path = p, pattern = paste0(".+",ft,".+csv"))
# list of the CSV files for the design 'd' and the strategy 'ft'
print(fl)
catalog_csv <-
bind_rows(catalog_csv,
tibble( design = d, c_strat = ft, path = p, fileName = fl))
}
}## [1] "design: A - strategy: dbh"
## [1] "fc01A_dbh.csv" "fc02A_dbh.csv" "fc03A_dbh.csv" "fc08A_dbh.csv"
## [5] "fc20A_dbh.csv" "fc21A_dbh.csv" "fc22A_dbh.csv" "fc23A_dbh.csv"
## [9] "fc25A_dbh.csv" "fc26A_dbh.csv" "fc37A_dbh.csv" "fc42A_dbh.csv"
## [13] "fc44A_dbh.csv" "fc48A_dbh.csv" "fc49A_dbh.csv"
## [1] "design: A - strategy: cyl"
## [1] "fc01A_cylinder.csv" "fc02A_cylinder.csv" "fc03A_cylinder.csv"
## [4] "fc08A_cylinder.csv" "fc20A__cylinder.csv" "fc21A_cylinder.csv"
## [7] "fc22A_cylinder.csv" "fc23A_cylinder.csv" "fc25A_cylinder.csv"
## [10] "fc26A_cylinder.csv" "fc37A_cylinder.csv" "fc42A_cylinder.csv"
## [13] "fc44A_cylinder.csv" "fc48A_cylinder.csv" "fc49A_cylinder.csv"
## [1] "design: AB - strategy: dbh"
## [1] "fc01AB_dbh.csv" "fc02AB_dbh.csv" "fc03AB_dbh.csv" "fc08AB_dbh.csv"
## [5] "fc20AB_dbh.csv" "fc21AB_dbh.csv" "fc22AB_dbh.csv" "fc23AB_dbh.csv"
## [9] "fc25AB_dbh.csv" "fc26AB_dbh.csv" "fc37AB_dbh.csv" "fc42AB_dbh.csv"
## [13] "fc44AB_dbh.csv" "fc48AB_dbh.csv" "fc49AB_dbh.csv"
## [1] "design: AB - strategy: cyl"
## [1] "fc01AB_cylinder.csv" "fc02AB_cylinder.csv"
## [3] "fc03AB_cylinder.csv" "fc08AB_cylinder.csv"
## [5] "fc20AB_cylinder.csv" "fc21AB_cylinder.csv"
## [7] "fc22AB_cylinder (1).csv" "fc23AB_cylinder.csv"
## [9] "fc25AB_cylinder.csv" "fc26AB_cylinder.csv"
## [11] "fc37AB_cylinder.csv" "fc42AB_cylinder.csv"
## [13] "fc44AB_cylinder.csv" "fc48AB_cylinder.csv"
## [15] "fc49AB_cylinder.csv"catalog_csv <- catalog_csv %>%
add_column(.before = 1, IdFC = as.numeric(str_sub(catalog_csv$fileName, 3,4)))
cat_all <- catalog_FC %>%
select(starts_with("Id_")) %>%
select(-contains('fusti'), -contains('staz')) %>%
gather(key = "Id_key", value = "Id_value", -Id_foto) %>%
filter(!is.na(Id_value)) %>%
mutate(design = str_sub(Id_key, ifelse(str_sub(Id_key, -1)=="A", -1, -2)),
c_strat = ifelse(str_sub(Id_key, 4, 6)=="dbh", !!ps[1], !!ps[2])) %>%
rename(IdFC = Id_foto) %>%
full_join(catalog_csv, .)## Joining, by = c("IdFC", "design", "c_strat")# %>% filter(IdFC == 42) manca Id_dbhABlibrary(data.table)## data.table 1.11.8 Latest news: r-datatable.com##
## Attaching package: 'data.table'## The following objects are masked from 'package:dplyr':
##
## between, first, last## The following object is masked from 'package:purrr':
##
## transposeretrieve_masures <- function (path, fileName, Id_value) {
print(paste(fileName, "Id_value = ", Id_value))
c_strat <- if_else(str_detect(fileName, "dbh"), "dbh",
if_else(str_detect(fileName, "cyl"), "cyl",
NA_character_))
if (is.na(c_strat)) stop("Missing CompuTree strategy!")
switch (c_strat,
dbh = {delim <- ";" ; skip <- 0},
cyl = {delim <- "\t"; skip <- 1})
rd <- suppressMessages(
suppressWarnings(
read_delim(paste0(path,"\\", fileName), delim = delim, skip = skip)))
# print(rd[1,1:10])
switch (c_strat,
dbh = {
# print(paste("searching for:", Id_value, "in:", fileName))
out <- rd %>%
filter(ID == Id_value) %>%
select(Rayonducercle, CenterZ, MinZ, MaxZ, SizeZ, Rayonducercle, Erreurdajustementducercle)
# out <- if_else(nrow(out) > 0, unlist(out), NaN)
},
cyl = {
delim <- "\t"
fn <- paste0(path,"\\", fileName)
rd <- suppressMessages(
suppressWarnings(
read_delim(fn , delim = delim, skip = 0, n_max = 2, col_names = F)))
rd <- t(rd)
cols2read <- c(1, 2, which(!is.na(str_locate(rd, "Up to Branch Order 0")[,1])))
rd <- fread(file = fn, sep = delim, header = T, skip = 1, select = cols2read)
tree_id <- rd[ID == Id_value,]$ParentID
out <- rd[ParentID == tree_id & !is.na(ID),]
# print(out)
}
)
# print(out)
return(out)
}
# test
# with(cat_all[16,], retrieve_masures(path, fileName, Id_value))
meas0 <- cat_all %>%
select(path, fileName, Id_value) %>%
mutate(meas = pmap(., retrieve_masures)) %>%
mutate(nr = map_int(meas, nrow)) %>%
select(-path) %>%
full_join(select(cat_all, fileName, IdFC, design, c_strat), .) %>%
full_join(select(catalog_FC, Id_foto, specie, DBH), ., by = c("Id_foto" = "IdFC")) %>%
rename(IdFC = Id_foto)## [1] "fc01A_dbh.csv Id_value = 4334336"
## [1] "fc02A_dbh.csv Id_value = 5421038"
## [1] "fc03A_dbh.csv Id_value = 6609090"
## [1] "fc08A_dbh.csv Id_value = 7684519"
## [1] "fc20A_dbh.csv Id_value = 962857"
## [1] "fc21A_dbh.csv Id_value = 707659"
## [1] "fc22A_dbh.csv Id_value = 1991893"
## [1] "fc23A_dbh.csv Id_value = 838138"
## [1] "fc25A_dbh.csv Id_value = 1581415"
## [1] "fc26A_dbh.csv Id_value = 2191494"
## [1] "fc37A_dbh.csv Id_value = 3099991"
## [1] "fc42A_dbh.csv Id_value = 660024"
## [1] "fc44A_dbh.csv Id_value = 205670"
## [1] "fc48A_dbh.csv Id_value = 1436936"
## [1] "fc49A_dbh.csv Id_value = 2197239"
## [1] "fc01A_cylinder.csv Id_value = 3503082"
## [1] "fc02A_cylinder.csv Id_value = 4348386"
## [1] "fc03A_cylinder.csv Id_value = 5445123"
## [1] "fc08A_cylinder.csv Id_value = 6626845"
## [1] "fc20A__cylinder.csv Id_value = 13176"
## [1] "fc21A_cylinder.csv Id_value = 13207"
## [1] "fc22A_cylinder.csv Id_value = 722191"
## [1] "fc23A_cylinder.csv Id_value = 20144"
## [1] "fc25A_cylinder.csv Id_value = 855963"
## [1] "fc26A_cylinder.csv Id_value = 1589562"
## [1] "fc37A_cylinder.csv Id_value = 2206920"
## [1] "fc42A_cylinder.csv Id_value = 10475"
## [1] "fc44A_cylinder.csv Id_value = 5582"
## [1] "fc48A_cylinder.csv Id_value = 505427"
## [1] "fc49A_cylinder.csv Id_value = 1450620"
## [1] "fc01AB_dbh.csv Id_value = 1236559"
## [1] "fc02AB_dbh.csv Id_value = 1234171"
## [1] "fc03AB_dbh.csv Id_value = 1721759"
## [1] "fc08AB_dbh.csv Id_value = 1257462"
## [1] "fc20AB_dbh.csv Id_value = 1081551"
## [1] "fc21AB_dbh.csv Id_value = 833611"
## [1] "fc22AB_dbh.csv Id_value = 1307437"
## [1] "fc23AB_dbh.csv Id_value = 889794"
## [1] "fc25AB_dbh.csv Id_value = 1684460"
## [1] "fc26AB_dbh.csv Id_value = 938556"
## [1] "fc37AB_dbh.csv Id_value = 1728859"
## [1] "fc42AB_dbh.csv Id_value = NA"
## [1] "fc44AB_dbh.csv Id_value = 1587040"
## [1] "fc48AB_dbh.csv Id_value = 2540970"
## [1] "fc49AB_dbh.csv Id_value = 3488621"
## [1] "fc01AB_cylinder.csv Id_value = 21506"
## [1] "fc02AB_cylinder.csv Id_value = 21526"
## [1] "fc03AB_cylinder.csv Id_value = 46045"
## [1] "fc08AB_cylinder.csv Id_value = 14517"
## [1] "fc20AB_cylinder.csv Id_value = 14667"
## [1] "fc21AB_cylinder.csv Id_value = 12691"
## [1] "fc22AB_cylinder (1).csv Id_value = 22929"
## [1] "fc23AB_cylinder.csv Id_value = 15416"
## [1] "fc25AB_cylinder.csv Id_value = 909137"
## [1] "fc26AB_cylinder.csv Id_value = 11086"
## [1] "fc37AB_cylinder.csv Id_value = 955996"
## [1] "fc42AB_cylinder.csv Id_value = 1740910"
## [1] "fc44AB_cylinder.csv Id_value = 1316495"
## [1] "fc48AB_cylinder.csv Id_value = 1605681"
## [1] "fc49AB_cylinder.csv Id_value = 2561522"## Joining, by = "fileName"dbh_meas <- meas0 %>%
filter(c_strat == 'dbh' & nr > 0) %>%
unnest %>%
bind_rows(filter(select(meas0, -meas), c_strat == 'dbh' & nr == 0)) %>%
arrange(specie, IdFC, design, c_strat) %>%
mutate(M_dbh = DBH, TPPC_dbh = Rayonducercle*200) %>%
select(-DBH, -Rayonducercle)
library(ggpmisc)## For news about 'ggpmisc', please, see https://www.r4photobiology.info/## For on-line documentation see https://docs.r4photobiology.info/ggpmisc/formula <- y ~ x
dbh_meas %>%
filter(!is.na(TPPC_dbh)) %>%
# filter(TPPC_dbh < 60) %>%
# c'è una osservazione fuori range!!
ggplot(aes(TPPC_dbh, M_dbh, shape = design, linetype = design, colour= design)) +
geom_abline(slope = 1, intercept = 0, colour = 'brown', size = 1.5) +
geom_point() +
scale_shape_manual(values=c(1,2)) +
geom_smooth(method=lm, formula = formula) +
stat_poly_eq(aes(label = paste(..eq.label.., ..adj.rr.label.., sep = "~~~~")), formula = formula, parse = TRUE, size = 8) +
theme(legend.justification=c(1,0), legend.position=c(1,0),
legend.title = element_text(size=20, face = "bold"),
legend.text = element_text(size = 20, face = "bold"))
lmod <- function(df){
lm(DBH ~ I(TPPC_dbh), data = df)
}
# Test H0 : M_dbh = 0 + 1 * TPPC_dbh
lin_mods <- dbh_meas %>%
filter(!is.na(TPPC_dbh)) %>%
# filter(TPPC_dbh < 60) %>%
mutate(d = TPPC_dbh - M_dbh) %>%
group_by(design) %>%
nest() %>%
mutate(lmod = map(data, ~lm(M_dbh ~ TPPC_dbh, data = .)),
lmod1 = map(data, ~lm(M_dbh ~ 1+TPPC_dbh+offset(TPPC_dbh), data = .))) %>%
mutate(ellip90 = map(lmod, ~ellipse::ellipse(., level = .90)),
ellip99 = map(lmod, ~ellipse::ellipse(., level = .99)))
# tramite H0 : M_dbh = p + (q + 1) * TPPC_dbh, con p = q = 0
# map(lin_mods$lmod1, summary)
# tramite ellisse di confidenza congiunta: il punto (0, 1) è interno all'ellisse?
lin_mods %>%
select(-data, -lmod, -lmod1) %>%
mutate(ellip90 = map(ellip90, ~as.tibble(.)),
ellip99 = map(ellip99, ~as.tibble(.))) %>%
unnest %>%
{bind_rows(transmute(., design, conf.level =.9, int = `(Intercept)`, spl = TPPC_dbh),
transmute(., design, conf.level =.99, int = `(Intercept)1`, spl = TPPC_dbh1))} %>%
mutate(conf.level = as.factor(conf.level), design = as.factor(design)) %>%
ggplot(aes(int, spl)) +
geom_point(aes(x =0, y =1), shape = 3, size = 10) +
geom_text(aes(x =0, y =1), label ="M_dbh = TPPC_dbh * 1 + 0", position = position_nudge(x = 3, y = 0), angle = -45, size = 8) +
geom_path(aes(linetype = conf.level, color = design), size = 1) +
scale_linetype_manual(values=c("dotted", "dashed")) +
labs(title = "H0: TCPP_dbh is equal to M_dbh",
subtitle = "Confidence ellipse for M_dbh = TCPP_dbh * p + q") +
xlab("intercept 'q'") + ylab("slope 'p'")
# può essere utile la 'stima per differenza'?
test <- function(x) {
t.test(x$TPPC_dbh - x$M_dbh, mu = 0, conf.level = .99) %>% broom::tidy()
}
see <- function(df) {
x <- df$TPPC_dbh - df$M_dbh
if(any(is.na(x))) stop("NA non ammessi!")
sqrt((var(x)/(length(x) -1)))
}
dbh_meas %>%
filter(!is.na(TPPC_dbh)) %>%
# filter(TPPC_dbh < 60) %>%
group_by(design) %>%
nest() %>%
mutate(testIFzero = map(data, test),
sd.est = map_dbl(data, see)) %>%
select(-data) %>%
unnest %>%
mutate('SEE%' = sprintf("%.1f%%", 100 * sd.est / estimate)) %>%
select(design, conf.low, conf.high, 'SEE%') %>%
kableExtra::kable(caption = "99% confidence test for 'T-PPC_dbh - M_dbh'") %>%
kableExtra::kable_styling()| design | conf.low | conf.high | SEE% |
|---|---|---|---|
| A | -0.3963355 | 1.683777 | 56.2% |
| AB | -0.2592402 | 2.344321 | 43.0% |
cyl_meas <- meas0 %>%
filter(c_strat == 'cyl') %>%
select(IdFC, design, meas) %>%
unnest
mins <-
cyl_meas %>%
group_by(IdFC, design) %>%
summarise(minMinZ = min(`Min Z`), minCenterZ = min(`Center Z`), minMaxZ = min(`Max Z`)) %>%
ungroup %>%
inner_join(rename(cyl_meas, minMinZ = `Min Z`)) %>%
# filter(minCenterZ != `Center Z` | minMaxZ != `Max Z`)
# filter(`Max Z` - minMinZ != `Size Z`)
# per verifica, tutti i minimi corrispondono! e Size = Max - Min
inner_join(select(dbh_meas, IdFC, design, CenterZ)) %>%
transmute(IdFC = IdFC, design = design, SoilZ = CenterZ - 1.3, minMinZ = minMinZ, delta = minMinZ - SoilZ )## Joining, by = c("IdFC", "design", "minMinZ")## Joining, by = c("IdFC", "design")formula <- y ~ 1 + offset(x)
mins %>% ggplot(aes(minMinZ, SoilZ, colour = design)) + geom_point() +
geom_smooth(method=lm, formula = formula) +
stat_poly_eq(aes(label = paste(..eq.label.., sep = "~~~~")),
formula = formula, parse = TRUE)## Warning: Removed 1 rows containing non-finite values (stat_smooth).## Warning: Removed 1 rows containing non-finite values (stat_poly_eq).## Warning: Removed 1 rows containing missing values (geom_point).
# min('Min Z') risulta posizionato 14-15 cm sotto SoilZ = (CenterZ)_dbh - 1.3 c
cyl_meas %>%
inner_join(mins) %>%
mutate(h_sez = `Center Z` - minMinZ, d_sez = Rayon * 200) %>%
ggplot(aes(d_sez, h_sez, colour = design)) +
geom_path() +
facet_wrap( ~ IdFC, scales="free") +
theme(legend.justification=c(1,0), legend.position=c(1,0),
legend.title = element_text(size=20, face = "bold"),
legend.text = element_text(size = 20, face = "bold"))## Joining, by = c("IdFC", "design")
rmarkdown::render("ProcessingComputreeOutput.Rmd")
markdown::rpubsUpload("T-PPC test Ver.1", "ProcessingComputreeOutput.html",
id = "https://api.rpubs.com/api/v1/document/436483/e2a62534558d4c39ae0f92b5067ced8d")