Determinando a distribuição das Árvores com o índice de Morisita

---

Packages

---

library(spatstat)
library(sf)

---

Gerar dados de exemplo para árvores em dois grupos (agregado e disperso)

---

# Grupo 1: Árvores agregadas (clustered)
set.seed(1)
coords_agregado = data.frame(
  x = rnorm(50, mean = 10, sd = 1),  # coordenadas agregadas
  y = rnorm(50, mean = 10, sd = 1)
)

# Grupo 2: Árvores dispersas
set.seed(2)
coords_disperso = data.frame(
  x = runif(50, min = 5, max = 20),  # coordenadas mais espalhadas
  y = runif(50, min = 5, max = 20)
)

---

Definir os polígonos (áreas de pastagem) para cada grupo

---

# Polígono 1 para o grupo agregado
poligono1 = owin(xrange = c(8, 12), yrange = c(8, 12))

# Polígono 2 para o grupo disperso
poligono2 = owin(xrange = c(5, 20), yrange = c(5, 20))

---

Converter as coordenadas para objetos espaciais do spatstat

---

# Pontos agregados
pontos_agregado = ppp(coords_agregado$x, coords_agregado$y, window = poligono1)

Warning: 3 points were rejected as lying outside the specified window

str(pontos_agregado)

List of 5
 $ window    :List of 4
  ..$ type  : chr "rectangle"
  ..$ xrange: num [1:2] 8 12
  ..$ yrange: num [1:2] 8 12
  ..$ units :List of 3
  .. ..$ singular  : chr "unit"
  .. ..$ plural    : chr "units"
  .. ..$ multiplier: num 1
  .. ..- attr(*, "class")= chr "unitname"
  ..- attr(*, "class")= chr "owin"
 $ n         : int 47
 $ x         : num [1:47] 9.37 10.18 9.16 11.6 10.33 ...
 $ y         : num [1:47] 10.4 9.39 10.34 8.87 11.43 ...
 $ markformat: chr "none"
 - attr(*, "class")= chr "ppp"
 - attr(*, "rejects")=List of 5
  ..$ window    :List of 5
  .. ..$ type  : chr "polygonal"
  .. ..$ xrange: num [1:2] 7.64 12
  .. ..$ yrange: num [1:2] 8 12.5
  .. ..$ bdry  :List of 1
  .. .. ..$ :List of 2
  .. .. .. ..$ x: num [1:6] 11.62 9.13 7.64 7.88 9.99 ...
  .. .. .. ..$ y: num [1:6] 12.54 12.09 10.01 8.99 8.06 ...
  .. ..$ units :List of 3
  .. .. ..$ singular  : chr "unit"
  .. .. ..$ plural    : chr "units"
  .. .. ..$ multiplier: num 1
  .. .. ..- attr(*, "class")= chr "unitname"
  .. ..- attr(*, "class")= chr "owin"
  ..$ n         : int 3
  ..$ x         : num [1:3] 11.51 7.79 10.59
  ..$ y         : num [1:3] 12.4 10 12.2
  ..$ markformat: chr "none"
  ..- attr(*, "class")= chr "ppp"

# Pontos dispersos
pontos_disperso = ppp(coords_disperso$x, coords_disperso$y, window = poligono2)
str(pontos_disperso)

List of 5
 $ window    :List of 4
  ..$ type  : chr "rectangle"
  ..$ xrange: num [1:2] 5 20
  ..$ yrange: num [1:2] 5 20
  ..$ units :List of 3
  .. ..$ singular  : chr "unit"
  .. ..$ plural    : chr "units"
  .. ..$ multiplier: num 1
  .. ..- attr(*, "class")= chr "unitname"
  ..- attr(*, "class")= chr "owin"
 $ n         : int 50
 $ x         : num [1:50] 7.77 15.54 13.6 7.52 19.16 ...
 $ y         : num [1:50] 5.11 5.22 15.25 18.95 9.13 ...
 $ markformat: chr "none"
 - attr(*, "class")= chr "ppp"

---

Função para calcular o Índice de Morisita

---

calcular_morisita = function(pontos, num_celulas) {
  # Dividir o polígono em uma grade com `num_celulas` por `num_celulas`
  grid = quadratcount(pontos, nx = num_celulas, ny = num_celulas)
  # Número total de árvores
  N = sum(grid)
  # Número de células
  n = num_celulas^2
  # Contagem de árvores em cada célula
  X_i = as.vector(grid)
  # Calcular o Índice de Morisita
  I_M = (n * sum(X_i * (X_i - 1))) / (N * (N - 1))
  return(I_M)
}
calcular_morisita

function(pontos, num_celulas) {
  # Dividir o polígono em uma grade com `num_celulas` por `num_celulas`
  grid = quadratcount(pontos, nx = num_celulas, ny = num_celulas)
  # Número total de árvores
  N = sum(grid)
  # Número de células
  n = num_celulas^2
  # Contagem de árvores em cada célula
  X_i = as.vector(grid)
  # Calcular o Índice de Morisita
  I_M = (n * sum(X_i * (X_i - 1))) / (N * (N - 1))
  return(I_M)
}

---

Calcular o Índice de Morisita para cada grupo

---

indice_morisita_agregado = calcular_morisita(pontos_agregado, num_celulas = 5)
indice_morisita_agregado

[1] 1.480111

indice_morisita_disperso = calcular_morisita(pontos_disperso, num_celulas = 5)
indice_morisita_disperso

[1] 1.142857

# Exibir os resultados
cat("Índice de Morisita (Agregado):", indice_morisita_agregado, "\n")

Índice de Morisita (Agregado): 1.480111 

cat("Índice de Morisita (Disperso):", indice_morisita_disperso, "\n")

Índice de Morisita (Disperso): 1.142857 

Valores de Morisita maiores que 1 indicam agregação, aproximadamente 1 indicam aleatoriedade, e menores que 1 indicam dispersão.

---

Visualizar a distribuição espacial com a grade sobreposta

---

# Grupo Agregado
plot(pontos_agregado, main = "Grupo Agregado", pch = 20)

grid_agregado = quadratcount(pontos_agregado, nx = 5, ny = 5)  # Cria a grade
# Plotar a grade em um gráfico separado
plot(grid_agregado, main = "Contagem de Quadrantes (Agregado)", col = "red")

# Grupo Disperso
plot(pontos_disperso, main = "Grupo Disperso", pch = 20)

grid_disperso = quadratcount(pontos_disperso, nx = 5, ny = 5)  # Cria a grade
# Plotar a grade em um gráfico separado
plot(grid_disperso, main = "Contagem de Quadrantes (Disperso)", col = "red")

---

Tentando adaptar o indice de morisita usando o raio

---

Data

set.seed(1)
coords_agregado = data.frame(
  x = rnorm(50, mean = 10, sd = 1),
  y = rnorm(50, mean = 10, sd = 1),
  raio = runif(50, 0.1, 0.5)
)
coords_agregado$superficie_coberta = pi * (coords_agregado$raio^2)
coords_agregado

set.seed(2)
coords_disperso = data.frame(
  x = runif(50, min = 5, max = 20),
  y = runif(50, min = 5, max = 20),
  raio = runif(50, 0.1, 0.5)
)
coords_disperso$superficie_coberta = pi * (coords_disperso$raio^2)
coords_disperso

---

Definir limites dos polígonos

---

limites_agregado = owin(xrange = c(8, 12), yrange = c(8, 12))
limites_disperso = owin(xrange = c(5, 20), yrange = c(5, 20))

---

Função para calcular o Índice de Morisita considerando superfície coberta

---

calcular_morisita_superficie = function(coords, num_celulas, poligono) {
  largura_celula = (poligono$xrange[2] - poligono$xrange[1]) / num_celulas
  altura_celula = (poligono$yrange[2] - poligono$yrange[1]) / num_celulas
  # Inicializar a grade
  grid_superficie = matrix(0, nrow = num_celulas, ncol = num_celulas)
  # Preencher a grade com a superfície coberta
  for (i in 1:nrow(coords)) {
    x_pos = min(num_celulas, max(1, floor((coords$x[i] - poligono$xrange[1]) / largura_celula) + 1))
    y_pos = min(num_celulas, max(1, floor((coords$y[i] - poligono$yrange[1]) / altura_celula) + 1))
    grid_superficie[x_pos, y_pos] <- grid_superficie[x_pos, y_pos] + coords$superficie_coberta[i]
  }
  # Calcular o Índice de Morisita modificado
  N = sum(grid_superficie) # Total de área coberta
  n = num_celulas^2        # Número de células
  X_i = as.vector(grid_superficie)
  
  I_M = if (N > 1) {
    (n * sum(X_i * (X_i - 1))) / (N * (N - 1))
  } else {
    0  # Caso especial onde N é muito pequeno
  }
  return(I_M)
}

---

Calcular o Índice de Morisita para ambos os grupos

---

indice_morisita_agregado = calcular_morisita_superficie(coords_agregado, num_celulas = 5, limites_agregado)
indice_morisita_disperso = calcular_morisita_superficie(coords_disperso, num_celulas = 5, limites_disperso)
cat("Índice de Morisita (Agregado com Superfície Coberta):", indice_morisita_agregado, "\n")

Índice de Morisita (Agregado com Superfície Coberta): 0.2192138 

cat("Índice de Morisita (Disperso com Superfície Coberta):", indice_morisita_disperso, "\n")

Índice de Morisita (Disperso com Superfície Coberta): 0.2984291 

---

Plot

---

plot_com_grade = function(coords, poligono, num_celulas, titulo) {
  # Configuração da janela gráfica
  plot(NA, xlim = poligono$xrange, ylim = poligono$yrange, xlab = "X", ylab = "Y", main = titulo)
  largura_celula <- (poligono$xrange[2] - poligono$xrange[1]) / num_celulas
  altura_celula <- (poligono$yrange[2] - poligono$yrange[1]) / num_celulas
  # Adicionar a grade
  for (i in 0:num_celulas) {
    abline(v = poligono$xrange[1] + i * largura_celula, col = "gray", lty = 2)
    abline(h = poligono$yrange[1] + i * altura_celula, col = "gray", lty = 2)
  }
  # Plotar cada ponto e o raio (raiz da superfície coberta)
  points(coords$x, coords$y, pch = 16, col = "blue")
  symbols(coords$x, coords$y, circles = sqrt(coords$superficie_coberta / pi), add = TRUE, inches = FALSE, fg = "red")
}

# Plotar grupo agregado e disperso com grade e raios
plot_com_grade(coords_agregado, limites_agregado, 5, "Grupo Agregado com Superfície Coberta")

plot_com_grade(coords_disperso, limites_disperso, 5, "Grupo Disperso com Superfície Coberta")