LOAD PACKAGES AND FUNCTION
Packages
pacman::p_load(tidyverse,Thermimage,RColorBrewer,extrafont) Functions
I put all of functions in other files for easy management. So I need to load them all before working.
source("torus.paracou.function.R")Set theme for plot
color.code <- brewer.pal(n = 6, name = "Set1")
my_color <-"#159858"
my_theme <- function(base_size =5, base_family = "sans"){
theme_bw(base_size = base_size, base_family = base_family) +
theme(
panel.grid.major = element_line(color = NA),
strip.background = element_rect(fill = my_color, color = my_color, size =0.5),
strip.text = element_text(face = "bold", size = 15, color = "white")
)
}
theme_set(my_theme())LOAD DATA OF PARACOU PLOTS
The plots 1-15 have the same size of 250m x 250m. Particularly, plot 16 is 500m x 500m. So I split the plot 16 into 4 subplots with 250m x 250m. So we have a total of 19 plots 250m.
At each plot, I divided into grids 10m x 10m quadrat.We have 625 quadrats of 10m for each plot. A total of 16 plots are 625*19 = 11875 quadrats.
At each quadrat 10m, I calculated the slope and altitude value based on the DEM 2015 raster. At the same time, I determined the topography habitat as well as the stems of the two species Cecropia obtusa and C. sciadophylla.
From habitat vector GIS we define 6 habitats: 1-Valley+Disturbed; 2-Valley; 3-Slope+Disturbed; 4-Slope; 5-Hill+Disturbed; 6-Hill
In the following table, each line is the data of a quadrat 10m.
Habitat_all <- read_csv("Quadrat10-OK3.csv")
Habitat_all## # A tibble: 11,875 x 13
## ID Plot ID_quadrat DEM2015 SLOPE_DEM2 Topo Topo_num Disturbed
## <int> <chr> <chr> <dbl> <dbl> <chr> <int> <int>
## 1 1 P01 P01-1 8.98 14.88633 Pente 2 2
## 2 2 P01 P01-2 9.67 15.57686 Pente 2 2
## 3 3 P01 P01-3 9.04 23.21437 Pente 2 2
## 4 4 P01 P01-4 8.30 24.04573 Bas fond 1 2
## 5 5 P01 P01-5 6.96 25.25134 Bas fond 1 2
## 6 6 P01 P01-6 5.90 25.26593 Bas fond 1 2
## 7 7 P01 P01-7 4.49 12.12114 Bas fond 1 2
## 8 8 P01 P01-8 4.06 9.16037 Bas fond 1 2
## 9 9 P01 P01-9 3.93 6.70303 Bas fond 1 2
## 10 10 P01 P01-10 4.07 5.74594 Bas fond 1 2
## # ... with 11,865 more rows, and 5 more variables: Habitat <chr>,
## # Habitat_num <int>, obtusa <int>, sciadophylla <int>, total <int>ANALYSIS
Calculate value observe for 3 case
1-Topo and Disturbed; 2-Topo; 3-Disturbed
#calculate value observe
(obs.habit <- Habitat_all %>% group_by(Habitat) %>% summarise(obtusa=sum(obtusa),
sciadophylla=sum(sciadophylla)))## # A tibble: 6 x 3
## Habitat obtusa sciadophylla
## <chr> <int> <int>
## 1 Bas fond-Disturbed 327 148
## 2 Bas fond-No 222 149
## 3 Pente-Disturbed 218 125
## 4 Pente-No 287 224
## 5 Plateau-Disturbed 150 138
## 6 Plateau-No 155 201(obs.topo <- Habitat_all %>% group_by(Topo) %>% summarise(obtusa=sum(obtusa),
sciadophylla=sum(sciadophylla)))## # A tibble: 3 x 3
## Topo obtusa sciadophylla
## <chr> <int> <int>
## 1 Bas fond 549 297
## 2 Pente 505 349
## 3 Plateau 305 339(obs.disturbed <- Habitat_all %>% group_by(Disturbed) %>% summarise(obtusa=sum(obtusa),
sciadophylla=sum(sciadophylla)))## # A tibble: 2 x 3
## Disturbed obtusa sciadophylla
## <int> <int> <int>
## 1 1 695 411
## 2 2 664 574So we have true data on the stem of trees per habitat of each species and their proportions.
TORUS Tranlastions idea
I take the example in plot P03 to present the idea of this method “torus translation test”
#Get data from plot P03
hab_p3 <- Habitat_all %>% filter(Plot =="P03")
hab_p3_mat <- matrix(hab_p3$Habitat_num, nrow=25)Show 625 quadrats of plot 03
b.image(hab_p3_mat,color.code)
with
Get number of stems of C. obtusa and show it in plot. This is observational data.
#get number individu of obtusa
obtusa <- matrix(hab_p3$obtusa, nrow=25)
b.image(hab_p3_mat,color.code)
b.show.numer.indi(obtusa)
We have a number of stem per habitat of C. obtusa at Plot 03:
obs.obtusa<-c()
for(i in 1:6){obs.obtusa[i]<-sum(obtusa[hab_p3_mat==i])}
names(obs.obtusa) <- paste0("Hab-",1:6)
obs.obtusa## Hab-1 Hab-2 Hab-3 Hab-4 Hab-5 Hab-6
## 20 18 22 30 27 9Next, we move the habitat map to the right with distance x = 5. While still fixing the position of the individual species of C. obtusa.
# Now we can make new.habitat by x=5 and y=0
newhab <- b.make.new.habitat.mat(hab_p3_mat, x=5, y=0, 25,25)
par(mfrow=c(1,2))
b.image(hab_p3_mat,color.code)
b.show.numer.indi(obtusa)
b.image(newhab,color.code)
b.show.numer.indi(obtusa)
We have Torus number of stem per habitat of C. obtusa at Plot 03. It has changed
torus.obtusa<-c()
for(i in 1:6){obs.obtusa[i]<-sum(obtusa[newhab==i])}
names(obs.obtusa) <- paste0("Hab-",1:6)
obs.obtusa## Hab-1 Hab-2 Hab-3 Hab-4 Hab-5 Hab-6
## 19 32 7 45 6 17Next, from the new habitat we can make three other habitat maps by: rotating 180 degrees, mirror and flip.
### make 3 other map from newhab rotate 180, mirror and flip
newhab.180 <- rotate180.matrix(newhab)
newhab.mirror <- mirror.matrix(newhab)
newhab.flip <- flip.matrix(newhab)
par(mfrow=c(2,2))
b.image(newhab,color.code,"New habitat")
b.image(newhab.180,color.code,"New habitat 180 rotate")
b.image(newhab.mirror,color.code,"New habitat mirror")
b.image(newhab.flip,color.code,"New habitat flip")
So with 1 times to moving of the habitat map, we have 4 new habitat maps. With each plot have 25x25 quadrat, we can move 625 times and create 2,500 habitat maps.
We do the same thing for the 18 plots remaining (250mx250m). In summary, we will have 2500 the number simulations of stems per habitat of all sample plots of each species.
