Análise da distribuição das árvores considerando área de copa: adaptação do índice de morisita

Pacotes necessários

library(sf)

Warning: package ‘sf’ was built under R version 4.3.3Linking to GEOS 3.11.2, GDAL 3.8.2, PROJ 9.3.1; sf_use_s2() is TRUE

library(dplyr)

Warning: package ‘dplyr’ was built under R version 4.3.3
Attaching package: ‘dplyr’

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

    filter, lag

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

    intersect, setdiff, setequal, union

library(readr)

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

library(spatstat.geom)

Warning: package ‘spatstat.geom’ 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
Warning: package ‘spatstat.univar’ was built under R version 4.3.3spatstat.univar 3.1-2
spatstat.geom 3.3-6

library(spatstat.explore)

Warning: package ‘spatstat.explore’ was built under R version 4.3.3Carregando pacotes exigidos: spatstat.random
Warning: package ‘spatstat.random’ was built under R version 4.3.3spatstat.random 3.3-2
Carregando pacotes exigidos: nlme

Attaching package: ‘nlme’

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

    collapse

spatstat.explore 3.4-2

library(spatstat.random)
library(ggplot2)

Use suppressPackageStartupMessages() to eliminate package startup messages

1- Importar dados

grid = st_read("C:/Users/Samsung/OneDrive/Documentos/1 - PósDoc/2 - FAPESP - Silvipastoril/3 - Execução/Dados/Fase 1/Distribuição/grid_20x20_singlepart.geojson")

Reading layer `grid_20x20_singlepart' from data source 
  `C:\Users\Samsung\OneDrive\Documentos\1 - PósDoc\2 - FAPESP - Silvipastoril\3 - Execução\Dados\Fase 1\Distribuição\grid_20x20_singlepart.geojson' 
  using driver `GeoJSON'
Simple feature collection with 3550 features and 2 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 676072.3 ymin: 8868882 xmax: 678343.2 ymax: 8870838
Projected CRS: WGS 84 / UTM zone 21S

grid

Simple feature collection with 3550 features and 2 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 676072.3 ymin: 8868882 xmax: 678343.2 ymax: 8870838
Projected CRS: WGS 84 / UTM zone 21S
First 10 features:
   area ID.Grid                       geometry
1   400       1 POLYGON ((677185.1 8870838,...
2   400       2 POLYGON ((677173.6 8870822,...
3   400       3 POLYGON ((677162.1 8870805,...
4   400       4 POLYGON ((677150.6 8870789,...
5   400       5 POLYGON ((677139.2 8870773,...
6   400       6 POLYGON ((677127.7 8870756,...
7   400       7 POLYGON ((677116.2 8870740,...
8   400       8 POLYGON ((677104.8 8870723,...
9   400       9 POLYGON ((677093.3 8870707,...
10  400      10 POLYGON ((677081.8 8870691,...

pontos = read_csv("C:/Users/Samsung/OneDrive/Documentos/1 - PósDoc/2 - FAPESP - Silvipastoril/3 - Execução/Dados/Fase 1/Especies_QGIS.csv")

Rows: 2457 Columns: 7── Column specification ───────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (1): Espécie
dbl (6): N, Piquete, LAT, LONG, LONG1, LAT1
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

pontos

copas = st_read("C:/Users/Samsung/OneDrive/Documentos/1 - PósDoc/2 - FAPESP - Silvipastoril/3 - Execução/Dados/Fase 1/GRIDs/Grids Sombreados.shp")

Reading layer `Grids Sombreados' from data source 
  `C:\Users\Samsung\OneDrive\Documentos\1 - PósDoc\2 - FAPESP - Silvipastoril\3 - Execução\Dados\Fase 1\GRIDs\Grids Sombreados.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 3550 features and 6 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 676072.3 ymin: 8868882 xmax: 678343.2 ymax: 8870838
Projected CRS: WGS 84 / UTM zone 21S

copas                

Simple feature collection with 3550 features and 6 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 676072.3 ymin: 8868882 xmax: 678343.2 ymax: 8870838
Projected CRS: WGS 84 / UTM zone 21S
First 10 features:
   area ID.Grid Área.sem Área.somb X..N..Sombr X.Sombreada                       geometry
1   400       1  400.000     0.000     100.000       0.000 MULTIPOLYGON (((677185.1 88...
2   400       2  400.000     0.000     100.000       0.000 MULTIPOLYGON (((677173.6 88...
3   400       3  400.000     0.000     100.000       0.000 MULTIPOLYGON (((677162.1 88...
4   400       4  400.000     0.000     100.000       0.000 MULTIPOLYGON (((677150.6 88...
5   400       5  400.000     0.000     100.000       0.000 MULTIPOLYGON (((677139.2 88...
6   400       6  400.000     0.000     100.000       0.000 MULTIPOLYGON (((677127.7 88...
7   400       7  382.960    17.040      95.740       4.260 MULTIPOLYGON (((677132.6 88...
8   400       8  358.073    41.927      89.518      10.482 MULTIPOLYGON (((677109.2 88...
9   400       9  396.272     3.728      99.068       0.932 MULTIPOLYGON (((677097.9 88...
10  400      10  383.605    16.395      95.901       4.099 MULTIPOLYGON (((677098.2 88...

2- Converter pontos para objeto sf e plotar

pontos_sf <- st_as_sf(pontos, coords = c("LONG1", "LAT1"), crs = 4326)
pontos_sf <- st_transform(pontos_sf, crs = st_crs(grid))
copas <- st_transform(copas, st_crs(grid))
area_sombreada_real = st_difference(grid, st_union(copas))

Warning: attribute variables are assumed to be spatially constant throughout all geometries

plot1=ggplot() +
  geom_sf(data = grid, fill = NA, color = "black") +
  geom_sf(data = pontos_sf, color = "red", size = 1) +
  theme_minimal() +
  labs(title = "Pontos de árvores sobre o grid 20x20 m")
plot1

plot2 = ggplot() +
  geom_sf(data = grid, fill = NA, color = "black", size=0.2) +
  geom_sf(data = area_sombreada_real, fill = "darkgreen", color = NA, alpha = 0.7) +
  theme_minimal() +
  labs(title = "Áreas reais ocupadas por copa dentro do grid 20x20 m")
plot2

Aparentemente deu tudo certo com os grids e com as coordenadas

3 - Metodos de de determinação de distribuição

3.1 - Clark-Evans

O índice de Clark-Evans é uma medida estatística descritiva usada para determinar se um padrão espacial de pontos é: Agregado (clusterizado), Aleatório ou Disperso (uniformemente distribuído). Formula: R= robs/rexp. Robs: distância média entre cada ponto e seu vizinho mais próximo (na amostra); Rexp: distância esperada entre vizinhos mais próximos se os pontos fossem distribuídos aleatoriamente em um plano com densidade constante. Interpretação dos valores de R: R < 1 Agregado; R ≈ 1 Aleatório (Poisson); R > 1 Disperso (uniforme). Limitações: Não funciona bem com poucos pontos.

calcular_padroes <- function(celula, pontos_celula) {
  if (nrow(pontos_celula) < 2) {
    return(data.frame(padrao = "Insuficiente", R = NA))
  }
  
  janela <- as.owin(st_geometry(celula))
  coords <- st_coordinates(pontos_celula)
  ppp_obj <- ppp(x = coords[,1], y = coords[,2], window = janela)

  ce <- clarkevans(ppp_obj)
  
  # Verifica se ce é uma lista com R
  R_val <- if ("R" %in% names(ce)) ce$R else ce
  
  padrao <- ifelse(R_val < 1, "Agregado",
                   ifelse(R_val > 1, "Disperso", "Aleatório"))
  
  return(data.frame(padrao = padrao, R = R_val))
}

# 1. Loop sobre células do grid
resultados <- lapply(1:nrow(grid), function(i) {
  celula <- grid[i, ]
  pontos_dentro <- pontos_sf[st_within(pontos_sf, celula, sparse = FALSE), ]
  dados <- calcular_padroes(celula, pontos_dentro)
  # Garante que sempre retorna uma única linha com o ID
  dados$ID <- i
  return(dados[1, ])
})

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

# 2. Juntar resultados ao grid
resultados_df <- do.call(rbind, resultados)
grid_resultado <- cbind(grid, resultados_df)
View(grid_resultado)

# 3. Exportar resultado para o QGIS
st_write(grid_resultado, "grid_resultado_clarkevans.shp", delete_dsn = TRUE)

Warning: GDAL Error 1: grid_resultado_clarkevans.shp does not appear to be a file or directory.

Deleting source `grid_resultado_clarkevans.shp' failed
Writing layer `grid_resultado_clarkevans' to data source 
  `grid_resultado_clarkevans.shp' using driver `ESRI Shapefile'
Writing 3550 features with 5 fields and geometry type Polygon.

3.2 - Índice de Morisita

O Morisita é ideal para padrões de contagem em subáreas, mesmo com poucas árvores por célula. Cada celula do grid é dividida em uma subdivisão virtual em grade (ex: 4x4 quadrados = 16 subáreas) e é contado quantas árvores que caem em cada uma. A partir disso, é calculado o índice de Morisita por célula.

names(grid_resultado1)

[1] "area_m2"         "id_grid"         "padrao_morisita" "id_morisita"     "id_celula"      
[6] "geometry"       

4 - Comparação entre Clark-Evans vs Morisita

3.3 - Índice de Morisita adaptado

Nesta abordagem fiz uma mudança na equação para considerar a area ocupada pela copa, em vez de apenas a coordenada.

# 4. Exportar para uso no QGIS
st_write(grid_resultado_area, "grid_resultado_morisita_area.gpkg", layer = "morisita_area", delete_layer = TRUE)

Writing layer `morisita_area' to data source `grid_resultado_morisita_area.gpkg' using driver `GPKG'
Writing 3550 features with 5 fields and geometry type Polygon.

5 - Plot Morisita adaptado