Introduction

2014 - Arborea - Pinus pinea plantation 4 (50 m * 50 m) plots Each plot is divided in (10 m side) squared subplots identified using the following scheme:

“A5” “B5” “C5” “D5” “E5”
“A4” “B4” “C4” “D4” “E4”
“A3” “B3” “C3” “D3” “E3”
“A2” “B2” “C2” “D2” “E2”
“A1” “B1” “C1” “D1” “E1”

In each subplot, point positions have been registered measuring the distances from the extremes of one side. Let ‘a’ be the side used as reference for the positions measurements. We will define a local cartesian sistem for each subplot using ‘a’ as ‘x axis’ and the left extreme of ‘a’, say ‘A’, as origin (0,0) The coordinates of the point on opposite extreme of ‘a’, say ‘B’, are hence (10, 0) For eac registered point, say ‘C’, the distances (AC), say ‘c’, and (BC), say ‘b’, have been measured and registered. To compute local cartesian coordinates of point ‘C’, say (lx, ly), one can exploit Erone’s formula, actually reversing it. (see: https://it.wikipedia.org/wiki/Formula_di_Erone ‘Dimostrazione alternativa’) \[ \begin{eqnarray} x & = & \frac {(a^{2}-b^{2}+c^{2})} {2a} \\ y & = & \frac {\sqrt{(4a^{2}c^{2}-(a^{2}-b^{2}+c^{2})^{2})}} {2a} \end{eqnarray} \]

Access to data

library(tidyverse)
## Warning: package 'tidyverse' was built under R version 3.3.3
## Loading tidyverse: ggplot2
## Loading tidyverse: tibble
## Loading tidyverse: tidyr
## Loading tidyverse: readr
## Loading tidyverse: purrr
## Loading tidyverse: dplyr
## Warning: package 'ggplot2' was built under R version 3.3.3
## Warning: package 'tibble' was built under R version 3.3.3
## Warning: package 'tidyr' was built under R version 3.3.3
## Warning: package 'readr' was built under R version 3.3.3
## Warning: package 'purrr' was built under R version 3.3.3
## Warning: package 'dplyr' was built under R version 3.3.3
## Conflicts with tidy packages ----------------------------------------------
## filter(): dplyr, stats
## lag():    dplyr, stats
library(magrittr)
## 
## Attaching package: 'magrittr'
## The following object is masked from 'package:purrr':
## 
##     set_names
## The following object is masked from 'package:tidyr':
## 
##     extract
library(googlesheets)
## Warning: package 'googlesheets' was built under R version 3.3.3
library(spatstat)
## Warning: package 'spatstat' was built under R version 3.3.3
## Loading required package: nlme
## 
## Attaching package: 'nlme'
## The following object is masked from 'package:dplyr':
## 
##     collapse
## Loading required package: rpart
## 
## spatstat 1.52-1       (nickname: 'Apophenia') 
## For an introduction to spatstat, type 'beginner'
## 
## Note: R version 3.3.2 (2016-10-31) is more than 9 months old; we strongly recommend upgrading to the latest version
suppressMessages(library(dplyr))
gs_ls()
## Warning: package 'bindrcpp' was built under R version 3.3.3
## # A tibble: 245 x 10
##                 sheet_title        author  perm version
##                       <chr>         <chr> <chr>   <chr>
##  1               DatiConcas dsch93tortoli    rw     old
##  2                   PLOT 4        scotti    rw     old
##  3                   PLOT 1        scotti    rw     old
##  4                   PLOT 3        scotti    rw     old
##  5                   PLOT 2        scotti    rw     old
##  6 BaseDati_AnimatedMovies…        scotti    rw     old
##  7                 Syllabus        scotti    rw     old
##  8 FATTURA N3_ACCONTO SANNA  irenepiredda    rw     old
##  9              GO_baseDati        scotti    rw     old
## 10              RotelleVer3        scotti    rw     old
## # ... with 235 more rows, and 6 more variables: updated <dttm>,
## #   sheet_key <chr>, ws_feed <chr>, alternate <chr>, self <chr>,
## #   alt_key <chr>
tpp <- gs_title("DatiConcas")
## Sheet successfully identified: "DatiConcas"
gs_ws_ls(tpp)
## [1] "RaccoltaFogliDati" "Correzioni"
ppwm <- tpp %>% gs_read(ws="RaccoltaFogliDati")
## Accessing worksheet titled 'RaccoltaFogliDati'.
## Parsed with column specification:
## cols(
##   plot = col_integer(),
##   subplot = col_character(),
##   id_sogetto = col_character(),
##   species = col_character(),
##   c = col_integer(),
##   b = col_integer(),
##   dbhNS = col_double(),
##   dbhEO = col_integer(),
##   h_tot = col_integer(),
##   h_crown_ins = col_integer(),
##   crown_rN = col_integer(),
##   crown_rE = col_integer(),
##   crown_rS = col_integer(),
##   crown_rO = col_integer(),
##   age = col_integer()
## )
# Nei tabelloni originari (1 per plot!!) il calcolo delle coordinate cartesiane è sviluppato nei fogli 'cooridinate' che, in generale, pescano i valori da 'Dati'. I casi per i quali 'b'+'c'<10m sono però corretti direttamente in 'cooridinate'!! La tabella 'corr' raccolglie le differenze riscontrate tra 'Dati' e 'cooridinate'.
corr <- tpp %>% gs_read(ws="Correzioni",  col_types=cols(id_sogetto = col_character()))
## Accessing worksheet titled 'Correzioni'.
# apply corrections
ppwm01 <- ppwm %>%
  left_join(corr) %>%
  mutate(b = ifelse(is.na(b_corr), b, b_corr),
         c = ifelse(is.na(c_corr), c, c_corr)) %>%
  select(-ends_with("_corr"))
## Joining, by = c("plot", "subplot", "id_sogetto")

Compute plot-wise coordinates (x, y)

# compute subplots origins
side <- 1:5
subplot_origin <- expand.grid(c=side, r=side) %>%
  mutate(subplot = paste(LETTERS[c], r, sep=""),
         x0 = 1000 * (c - 1),
         y0 = 1000 * (r - 1)) %>%
  select(-r, -c)

# distances are registered as cm (integers)
a = 10*100

ppwm02 <- ppwm01 %>%
# compute local cartesian coordiantes (lx, ly)
  mutate(lx = (a^2 - b^2 + c^2)/(2*a),
         ly = sqrt((4*a^2*c^2 - (a^2 - b^2 + c^2)^2))/(2*a)) %>%
#(alternativa di Laura) ly = sqrt((2*(a^2*b^2 + a^2*c^2 + b^2*c^2)- (a^4 + b^4 + c^4))/(2*a)))
# compute plot-wise coordinates (x, y)
  full_join(subplot_origin) %>%
  mutate(x = (x0 + lx), 
         y = (y0 + ly))
## Joining, by = "subplot"
# CREATE one PPP object for each 'plot'
side_l = a * length(side)
cf <- 1/100 # conversion of "cm" to "m" 
#               marks = .[,7:16],
Arborea_permanent_plots <- ppwm02 %>%
  group_by(plot) %>%
#  filter(plot==1 & subplot=="A1") %$%
  do(ppp = ppp(.$x * cf, .$y * cf, 
               window = owin(xrange=c(0, side_l * cf),
                             yrange=c(0, side_l * cf), 
                               unitname="m"))
  ) %>%
  ungroup()
## Warning: 1 point was rejected as lying outside the specified window
## Warning: 7 points were rejected as lying outside the specified window
## Warning: data contain duplicated points

## Warning: data contain duplicated points
## Warning: 1 point was rejected as lying outside the specified window
## Warning: data contain duplicated points

Verify distribution of points

plot(as.solist(Arborea_permanent_plots$ppp), main="2014Concas-Monitoring Pinus pinea regeneration, Arborea", main.panel=paste("Plot",Arborea_permanent_plots$plot))
## Warning: Interpretation of arguments maxsize and markscale has changed (in
## spatstat version 1.37-0 and later). Size of a circle is now measured by its
## diameter.
## Warning in plot.ppp(x = structure(list(window = structure(list(type =
## "rectangle", : 1 illegal points also plotted
## Warning in plot.ppp(x = structure(list(window = structure(list(type =
## "rectangle", : 7 illegal points also plotted
## Warning in plot.ppp(x = structure(list(window = structure(list(type =
## "rectangle", : 1 illegal points also plotted

plot.with.Axis <- function(i) {
  ungroup(Arborea_permanent_plots) %$% 
    plot(.$ppp[[i]], main=paste("Plot", .$plot[[i]]))
  at <- seq(0,50,10)
  Axis(side=1, at=at)
  Axis(side=4, at=at)
}
op <- par(no.readonly=T)
par(pin= c(par()$pin[1],.7 * par()$fin[2]))
for(i in 1:4) plot.with.Axis(i)
## Warning in plot.ppp(.$ppp[[i]], main = paste("Plot", .$plot[[i]])): 1
## illegal points also plotted

## Warning in plot.ppp(.$ppp[[i]], main = paste("Plot", .$plot[[i]])): 7
## illegal points also plotted

## Warning in plot.ppp(.$ppp[[i]], main = paste("Plot", .$plot[[i]])): 1
## illegal points also plotted

par(op)
  
#          dbh = ifelse(is.na(dbhEO), dbhNS, sqrt(mean(dbhNS^2, dbhEO^2)))