LOAD PACKAGES AND FUNCTIONS

Packages

pacman::p_load(tidyverse,maptools,spatialEco,RColorBrewer,extrafont)

Functions

I put all of functions in other files for easy management. So I need to load them all before working.

source("Point_translation_function.R")

Set theme for plot

my_color <-"#159858"
my_theme <- function(base_size =5, base_family = "CMU Serif" ){
    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 = 14, color = "white")
        )
}
theme_set(my_theme())

LOAD SHAPFILES OF PARACOU

Load

crswgsUTM <-CRS("+proj=utm +zone=22 +datum=WGS84 +units=m +no_defs")

path <- "D:/GIS/FTH/paracou_shpfiles/Rotate/"
individu <- readShapePoints(paste0(path,"Cecropia_individu_paracou5.shp"),proj4string=crswgsUTM,verbose=TRUE,repair = T)
## Shapefile type: Point, (1), # of Shapes: 2362
plot <- readShapePoly(paste0(path,"Plots.shp"),proj4string=crswgsUTM,verbose=TRUE,repair = T)
## Shapefile type: Polygon, (5), # of Shapes: 16
disturbed <- readShapePoly(paste0(path,"DisturbedAreas.shp"),proj4string=crswgsUTM,verbose=TRUE)
## Shapefile type: Polygon, (5), # of Shapes: 405
topo <- readShapePoly(paste0(path,"TopographicLevels_bached.shp"),proj4string=crswgsUTM,verbose=TRUE)
## Shapefile type: Polygon, (5), # of Shapes: 48
individu <- individu[individu$sp!="Indet.",]
disturbed@data$Type<-"Disturbed"

Show map

my.color<- brewer.pal(3,"Blues")
plot(plot, col="gray", xlab="Long", ylab="Lat");axis(1); axis(2); box()
topolist <-as.character(unique(topo$Type))
for(i in 1:length(topolist)){
    plot(topo[topo$Type==topolist[i],], col=my.color[i], add=T)
}
plot(disturbed, col="yellow", add=T, border=NA)
plot(individu, add=T, cex=0.5, pch=19)

IDEA

Move point A by a random vector V. But don’t move point A out of plot. The excess will be translated to the opposite side. B is the new position of the point.

TEST FOR PLOT 5 IN PARACOU

Subset data for plot 5

#Get value xmin, xmax, ymin, ymax of each plot
minmax<-rlist::list.load("minmax.rds")
p=5
subindi <- individu[individu$Plot==p,]
df <- subindi@data %>% dplyr::select(code,sp,Plot) %>% cbind(subindi@coords)
colnames(df) <- c("code","sp","Plot","x","y")
df$Plot <- as.character(df$Plot)

Test

We have plot 5

b.plot.map(5)

We make a random vector

set.seed(123)
b.plot.map(5)
vector <- b.random.vector(286800,582900,xlim=286400:286800,ylim = 582600:583200)

We move all of points in object df by this vector. Red points are the new position.

set.seed(123)
b.plot.map(5)
vector <- b.random.vector(286800,582900,xlim=286400:286800,ylim = 582600:583200)
for(i in 1:nrow(df)){
    b.move.one.point(df$x[i],df$y[i],vector,minmax[[df$Plot[i]]],plot = T)
}

We repet n times with n vector. Each time we will calculate the number of stem by habitat. We wil have n values which will be used for independent testing.

ANALYSIS FOR ALL PLOTS AT PARACOU

I do the same thing 2500 times for 16 plot at Paracou. Because this simulation will take long times. So I ran this before and saved the result in R. Now we just load and use it to continue the analysis. 

#load resultat from Toric translation Points
load("res_toric_translation2500.Rdata")

Observation data

obs.habit <- individu@data %>% group_by(sp,Topo,Disturbed) %>% tally(wt=NULL)
obs.topo <- individu@data %>% group_by(sp,Topo) %>% tally(wt=NULL)
obs.disturbed <- individu@data %>% group_by(sp,Disturbed) %>% tally(wt=NULL)

Test

For habitat topo and zone disturbed

C. obtusa

b.plot.density(obs.habit$n[1:6],obtusa_habit)
## # A tibble: 6 x 5
##              habitat   lower   upper   obs test_2_tail
##                <chr>   <dbl>   <dbl> <int>       <chr>
## 1 Bas fond-Disturbed 111.000 222.525   346           +
## 2       Bas fond-Non 204.000 320.525   206           n
## 3    Pente-Disturbed 120.000 192.000   216           +
## 4          Pente-Non 370.000 517.000   287           -
## 5  Plateau-Disturbed  66.000 125.525   152           +
## 6        Plateau-Non 193.475 329.000   152           -

C. sciadophylla

b.plot.density(obs.habit$n[7:12],sciado_habit)
## # A tibble: 6 x 5
##              habitat lower   upper   obs test_2_tail
##                <chr> <dbl>   <dbl> <int>       <chr>
## 1 Bas fond-Disturbed    83 146.000   151           +
## 2       Bas fond-Non   139 229.000   142           n
## 3    Pente-Disturbed    85 144.000   120           n
## 4          Pente-Non   272 402.000   239           -
## 5  Plateau-Disturbed    41 109.525   125           +
## 6        Plateau-Non   124 225.000   208           n

For habitat topography

C. obtusa

b.plot.density(obs.topo$n[1:3],obtusa_topo,ncol=1)
## # A tibble: 3 x 5
##    habitat lower upper   obs test_2_tail
##      <chr> <dbl> <dbl> <int>       <chr>
## 1 Bas fond   334   528   552           +
## 2    Pente   523   680   503           -
## 3  Plateau   282   432   304           n

C. sciadophylla

b.plot.density(obs.topo$n[4:6],sciado_topo,ncol=1)
## # A tibble: 3 x 5
##    habitat lower upper   obs test_2_tail
##      <chr> <dbl> <dbl> <int>       <chr>
## 1 Bas fond   241   357   293           n
## 2    Pente   377   519   359           -
## 3  Plateau   180   321   333           +

For zone disturbed

C. obtusa

b.plot.density(obs.disturbed$n[1:2],obtusa_disturbed)
## # A tibble: 2 x 5
##     habitat lower upper   obs test_2_tail
##       <chr> <dbl> <dbl> <int>       <chr>
## 1 Disturbed   349   481   714           +
## 2       Non   878  1010   645           -

C. sciadophylla

b.plot.density(obs.disturbed$n[3:4],sciado_disturbed)
## # A tibble: 2 x 5
##     habitat   lower   upper   obs test_2_tail
##       <chr>   <dbl>   <dbl> <int>       <chr>
## 1 Disturbed 245.475 351.000   396           +
## 2       Non 634.000 739.525   589           -