Análise da densidade das árvores considerando área de copa

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)
library(purrr)

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

View(grid)

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

Reading layer `Poligono_piquetes' from data source 
  `C:\Users\Samsung\OneDrive\Documentos\1 - PósDoc\2 - FAPESP - Silvipastoril\3 - Execução\Dados\Fase 1\Densidade\Poligono_piquetes.geojson' 
  using driver `GeoJSON'
Simple feature collection with 122 features and 4 fields
Geometry type: MULTIPOLYGON
Dimension:     XYZ
Bounding box:  xmin: 676096.9 ymin: 8868829 xmax: 678333.3 ymax: 8870821
z_range:       zmin: 0 zmax: 0
Projected CRS: WGS 84 / UTM zone 21S

View(Piquetes)

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...

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

Rows: 40 Columns: 5── Column specification ──────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (2): Densidade, Densidade real
dbl (3): Piquete, Grid, Rep
ℹ 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.

View(grids_selecionados)

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 - Densidade arborea

# 1. Gerar as novas geometrias invertidas (sombra real)
geometrias_corrigidas <- map2(st_geometry(grid), st_geometry(copas), st_difference)

# 2. Substituir a geometria no objeto copas1
copas1 <- copas
st_geometry(copas1) <- st_sfc(geometrias_corrigidas, crs = st_crs(copas))

# 3. Verificar quais ficaram vazias e substituir pela geometria original
geometrias_vazias <- st_is_empty(copas1)
st_geometry(copas1)[geometrias_vazias] <- st_geometry(grid)[geometrias_vazias]
View(copas1)

copas_classificado <- copas1 %>%
  mutate(classe_sombra = case_when(
    X.Sombreada <= 8 ~ "0 a 8%",
    X.Sombreada > 8.01 & X.Sombreada <= 16 ~ "8 a 16%",
    X.Sombreada > 16.01 & X.Sombreada <= 24 ~ "16 a 24%",
    X.Sombreada > 24 ~ "> 24%"
  ))
View(copas_classificado)

st_write(copas_classificado, "copas_classificado_sombra.gpkg", layer = "copas_sombra", delete_layer = TRUE)

Deleting layer `copas_sombra' using driver `GPKG'
Writing layer `copas_sombra' to data source `copas_classificado_sombra.gpkg' using driver `GPKG'
Writing 3550 features with 7 fields and geometry type Unknown (any).

#Plot
p2=ggplot(copas_classificado) +
  geom_sf(aes(fill = classe_sombra), color = "gray30", size = 0.1) +
  scale_fill_manual(values = c("0 a 8%" = "lightgray","8 a 16%" = "gold","16 a 24%" = "orange","> 24%" = "darkgreen"),name = "Porcentagem de Sombra") + theme_minimal() +
  labs(title = "Classificação dos Grids por Área Sombreada (%)",x = NULL, y = NULL) +
  theme(legend.position = "right",plot.title = element_text(size = 14, face = "bold"))
p2

8 - Filtrar grids e piquetes

st_write(pontos_clean, "pontos_filtrados.shp")

Warning: Field names abbreviated for ESRI Shapefile driver

Layer pontos_filtrados in dataset pontos_filtrados.shp already exists:
use either append=TRUE to append to layer or append=FALSE to overwrite layer

Erro: Dataset already exists.