Índice de Morisita de Novo Mundo

1 - Pacotes necessários

library(sf)

Linking to GEOS 3.11.2, GDAL 3.7.2, PROJ 9.3.0; sf_use_s2() is TRUE

library(spatstat)

Warning: package ‘spatstat’ was built under R version 4.3.3Carregando pacotes exigidos: spatstat.data
Warning: package ‘spatstat.data’ was built under R version 4.3.3Carregando pacotes exigidos: spatstat.univar
spatstat.univar 3.1-1
Carregando pacotes exigidos: spatstat.geom
Warning: package ‘spatstat.geom’ was built under R version 4.3.3spatstat.geom 3.3-3
Carregando pacotes exigidos: spatstat.random
Warning: package ‘spatstat.random’ was built under R version 4.3.3spatstat.random 3.3-2
Carregando pacotes exigidos: spatstat.explore
Warning: package ‘spatstat.explore’ was built under R version 4.3.3Carregando pacotes exigidos: nlme

Attaching package: ‘nlme’

The following object is masked from ‘package:lme4’:

    lmList

spatstat.explore 3.3-3
Carregando pacotes exigidos: spatstat.model
Warning: package ‘spatstat.model’ was built under R version 4.3.3Carregando pacotes exigidos: rpart
spatstat.model 3.3-2
Carregando pacotes exigidos: spatstat.linnet
Warning: package ‘spatstat.linnet’ was built under R version 4.3.3spatstat.linnet 3.2-2

spatstat 3.2-1 
For an introduction to spatstat, type ‘beginner’ 

library(dplyr)

Attaching package: ‘dplyr’

The following object is masked from ‘package:nlme’:

    collapse

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union

library(ggplot2)
library(openxlsx)

Warning: package ‘openxlsx’ was built under R version 4.3.3

library(tibble)

2 - Dados

grids=readxl::read_excel("C:/Users/Samsung/Documents/3 - PósDoc/2 - FAPESP - Silvipastoril/3 - Execução/Dados/Novo Mundo/Limites dos grids_R.xlsx")
grids$Y1=as.numeric(grids$Y1)
grids$X1=as.numeric(grids$X1)
grids$Y2=as.numeric(grids$Y2)
grids$X2=as.numeric(grids$X2)
grids$Y3=as.numeric(grids$Y3)
grids$X3=as.numeric(grids$X3)
grids$Y4=as.numeric(grids$Y4)
grids$X4=as.numeric(grids$X4)
str(grids)

tibble [170 × 9] (S3: tbl_df/tbl/data.frame)
 $ Grid: num [1:170] 1 2 3 4 5 6 7 8 9 10 ...
 $ Y1  : num [1:170] 8937053 8937039 8937026 8937013 8936999 ...
 $ X1  : num [1:170] 709421 709436 709451 709466 709481 ...
 $ Y2  : num [1:170] 8937068 8937054 8937041 8937027 8937014 ...
 $ X2  : num [1:170] 709434 709449 709464 709479 709494 ...
 $ Y3  : num [1:170] 8937039 8937026 8937013 8936999 8936986 ...
 $ X3  : num [1:170] 709436 709451 709466 709481 709495 ...
 $ Y4  : num [1:170] 8937054 8937041 8937027 8937014 8937001 ...
 $ X4  : num [1:170] 709449 709464 709479 709494 709509 ...

print(grids)

arvores=read.csv("C:/Users/Samsung/Documents/3 - PósDoc/2 - FAPESP - Silvipastoril/3 - Execução/Dados/Novo Mundo/Arvores.csv")
arvores

3 - Plotando os grids e as arvores dentro dos grids

# Ajustar os pontos do database para corrigir os limites dos grids
grids <- grids %>%
  mutate(
    X_temp = X3, Y_temp = Y3, # Salvar temporariamente os valores de X3 e Y3
    X3 = X4, Y3 = Y4,         # Substituir X3 e Y3 pelos valores de X4 e Y4
    X4 = X_temp, Y4 = Y_temp  # Substituir X4 e Y4 pelos valores salvos de X3 e Y3
  ) %>%
  select(-X_temp, -Y_temp)    # Remover colunas temporárias

# Criar os polígonos dos grids
grids_sf <- do.call(rbind, lapply(1:nrow(grids), function(i) {
  coords <- matrix(
    c(
      grids$X1[i], grids$Y1[i],
      grids$X2[i], grids$Y2[i],
      grids$X3[i], grids$Y3[i],
      grids$X4[i], grids$Y4[i],
      grids$X1[i], grids$Y1[i] # Fechar o polígono
    ),
    ncol = 2, byrow = TRUE
  )
  st_polygon(list(coords)) %>% st_sfc(crs = 32721) %>% st_sf(Grid = grids$Grid[i])
}))

# Converter as coordenadas das árvores para o sistema de coordenadas UTM
arvores_sf <- st_as_sf(arvores, coords = c("LONG", "LAT"), crs = 4326)
arvores_sf <- st_transform(arvores_sf, crs = 32721)

# Plotar os grids e as árvores
plot <- ggplot() +
  geom_sf(data = grids_sf, fill = NA, color = "blue", size = 0.5) +
  geom_sf(data = arvores_sf, color = "darkgreen", size = 0.5) +
  theme_minimal() +
  labs(title = "Distribuição de árvores nos grids", 
       x = "Oeste", y = "Sul")
plot

4 - Função para calcular o índice de Morisita

calcular_morisita <- function(pontos, num_celulas) {
  grid <- quadratcount(pontos, nx = num_celulas, ny = num_celulas)
  N <- sum(grid)
  n <- num_celulas^2
  X_i <- as.vector(grid)
  I_M <- (n * sum(X_i * (X_i - 1))) / (N * (N - 1))
  return(I_M)
}

# Calcular o índice de Morisita para cada grid
resultados <- lapply(1:nrow(grids_sf), function(i) {
  grid_poligono <- st_geometry(grids_sf[i, ])
  pontos_no_grid <- arvores_sf[st_within(arvores_sf, grid_poligono, sparse = FALSE), ]
  if (nrow(pontos_no_grid) > 1) {
    pontos_ppp <- as.ppp(st_coordinates(pontos_no_grid), W = as.owin(st_bbox(grid_poligono)))
    calcular_morisita(pontos_ppp, num_celulas = 5) # Ajuste o número de células conforme necessário
  } else {
    NA
  }
})

Warning: data contain duplicated pointsWarning: data contain duplicated pointsWarning: data contain duplicated pointsWarning: data contain duplicated points

grids_sf$Morisita_Index <- unlist(resultados)

# Classificar os grids como "Agregado" ou "Disperso"
grids_sf <- grids_sf %>%
  mutate(Classification = ifelse(Morisita_Index > 1, "Agregado", "Disperso"))
grids_sf

Simple feature collection with 170 features and 3 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 709421.1 ymin: 8936598 xmax: 709993.3 ymax: 8937127
Projected CRS: WGS 84 / UTM zone 21S
First 10 features:
   Grid                              . Morisita_Index Classification
1     1 POLYGON ((709421.1 8937053,...             NA           <NA>
2     2 POLYGON ((709436 8937039, 7...             NA           <NA>
3     3 POLYGON ((709450.8 8937026,...             NA           <NA>
4     4 POLYGON ((709465.7 8937013,...             NA           <NA>
5     5 POLYGON ((709480.5 8936999,...             NA           <NA>
6     6 POLYGON ((709495.4 8936986,...             NA           <NA>
7     7 POLYGON ((709510.3 8936972,...             NA           <NA>
8     8 POLYGON ((709525.1 8936959,...             NA           <NA>
9     9 POLYGON ((709540 8936946, 7...             NA           <NA>
10   10 POLYGON ((709554.8 8936932,...             NA           <NA>

print(as_tibble(grids_sf), n = Inf)

# Exportar os resultados para um arquivo XLSX
write.xlsx(st_drop_geometry(grids_sf), "Resultado_Morisita.xlsx", rownames = FALSE)

3 - Plot do índice de Morisita

# Plotar os grids e as árvores
plot1 <- ggplot() +
  geom_sf(data = grids_sf, aes(fill = Classification), color = "blue", size = 0.5) +
  geom_sf(data = arvores_sf, color = "darkgreen", size = 0.5) +
  scale_fill_manual(values = c("Agregado" = "red", "Disperso" = "green"), name = "Classificação") +
  theme_minimal() +
  labs(title = "Distribuição de árvores e classificação dos grids", 
       x = "Oeste", y = "Sul")
plot1