Geosaber
  • Potenciometrico
  • Visualização 3D

Curso - Análise Espacial e Interpolação de Dados no R (9 a 13/06)

Visualização 3D com Rayshader

Conteúdo

  • 1. Preparação dos dados
  • 2. Visualização 3D Básica com WebGL
  • 3. Superfície Potenciométrica
  • 4. Cálculo do Gradiente Hidráulico
  • 5. Geração das Setas de Fluxo
  • 5. Visualização 3D com rayshader
Curso de Análise Espacial e Interpolação de Dados no R

Dias 09 a 13 de Junho, das 19h as 21h

Inscrições abertas, vagas limitadas.

https://geosaber.com.br

1. Preparação dos dados

# Carregar as bibliotecas
library(terra)
library(gstat)
library(rayshader)
library(rgl)
library(sf)
library(stars)
library(sp)
library(ggplot2)

# Carregar o dataset meuse
data(meuse)
data(meuse.grid)
meuse_sf <- st_as_sf(
  meuse,
  coords = c("x", "y"), 
  crs = 28992  # EPSG 28992 = RD New (Holanda)
)
meuse_grid_sf <- st_as_sf(meuse.grid, coords = c("x", "y"), crs = 28992)

# Visualizar os pontos de elevação
ggplot(meuse, aes(x, y, color = elev)) + 
  geom_point(size = 3) +
  scale_color_viridis_c() +
  coord_equal() +
  labs(title = "Pontos de amostra - Elevação")

2. Visualização 3D Básica com WebGL

library(plotly)

# Dados simplificados
plot_data <- data.frame(
  x = st_coordinates(meuse_sf)[,1],
  y = st_coordinates(meuse_sf)[,2],
  z = meuse_sf$elev
)

# Plot 3D interativo
plot_ly(plot_data, 
        x = ~x, y = ~y, z = ~z,
        type = "scatter3d",
        mode = "markers",
        color = ~z,
        colors = viridisLite::viridis(100)) |>
  layout(scene = list(
    xaxis = list(title = "Longitude"),
    yaxis = list(title = "Latitude"),
    zaxis = list(title = "Elevação"))
  )

3. Superfície Potenciométrica

# Carregar pacotes
library(terra)
library(gstat)
library(sp)

# 1. Carregar e preparar os dados
data(meuse)
data(meuse.grid)

# Converter para objetos espaciais
meuse_sp <- SpatialPointsDataFrame(meuse[, c("x", "y")], 
                                 data = meuse,
                                 proj4string = CRS("+init=EPSG:28992"))

meuse_grid_sp <- SpatialPixels(SpatialPoints(meuse.grid[, c("x", "y")],
                                          proj4string = CRS("+init=EPSG:28992")))

# 2. Análise de variograma
variogram_model <- variogram(elev ~ 1, meuse_sp)
vgm_model <- fit.variogram(variogram_model, vgm(1, "Gau", 800, 0.1))

# 3. Execute a krigagem
krige_result <- krige(elev ~ 1, meuse_sp, newdata = meuse_grid_sp, model = vgm_model)
[using ordinary kriging]
plot(krige_result)

# 4. Converter para terra e selecionar apenas a predição
elev_krige <- rast(krige_result)[[1]]  # [[1]] seleciona a primeira camada
names(elev_krige) <- "elevation"

4. Cálculo do Gradiente Hidráulico

elev_layer <- elev_krige[[1]]
# Calcular componentes do gradiente (direção do fluxo)
gradient <- terrain(elev_layer, v = c("slope", "aspect"), unit = "radians")

# Componentes do vetor gradiente (PERPENDICULAR às isolinhas)
dx_rast <- sin(gradient$aspect) * gradient$slope  # Componente x
dy_rast <- cos(gradient$aspect) * gradient$slope   # Componente y

5. Geração das Setas de Fluxo

arrow_grid <- st_make_grid(meuse_sf, cellsize = 350, what = "centers") %>% 
  st_as_sf() %>% 
  vect()

# Extração e normalização dos vetores
arrow_data <- data.frame(
  x = crds(arrow_grid)[,1],
  y = crds(arrow_grid)[,2],
  dx = as.numeric(terra::extract(dx_rast, arrow_grid)[,2]),
  dy = as.numeric(terra::extract(dy_rast, arrow_grid)[,2])
) %>% 
  filter(!is.na(dx) & !is.na(dy)) %>%
  mutate(
    # Normalização para tamanho fixo (mantendo direção)
    norm = sqrt(dx^2 + dy^2),
    dx_norm = dx/norm * 200,  # 50 = comprimento fixo
    dy_norm = dy/norm * 200,
    xend = x + dx_norm,
    yend = y + dy_norm
  )

5. Visualização 3D com rayshader

library(ggplot2)
library(rayshader)
library(terra)

## 1. Dados: raster de krigagem (SpatRaster)
dados_raster <- as.data.frame(elev_krige, xy = TRUE)  # Converter para dataframe

## 2. Plot 2D (ggplot2) com elementos personalizados
gg_plot <- ggplot() +
  
  # Camada raster (krigagem)
  geom_raster(data = dados_raster,
              aes(x = x, y = y, fill = elevation)) +
  # Gradiente natural (terreno)
  scale_fill_gradientn(colors = terrain.colors(256), name = "Elevação (m)") +
  
    # Vetores de fluxo (setas)
  geom_segment(
    data = arrow_data,
    aes(x = x, y = y, xend = xend, yend = yend),
    color = "red",
    linewidth = 0.6,
    arrow = arrow(
      angle = 20,
      length = unit(0.1, "inches"),
      type = "open"
    )
  ) +
  
  # Pontos de amostragem (ex: dados do 'meuse')
  geom_point(data = meuse, aes(x = x, y = y), 
             color = "red", size = 1.5, alpha = 0.7) +
  
  # Contornos (opcional)
  geom_contour(data = dados_raster, 
               aes(x = x, y = y, z = elevation),
               color = "black", alpha = 0.3) +
  
  # Layout
  labs(title = "Elevação com Setas de Fluxo",
       x = "Coordenada X (UTM)", 
       y = "Coordenada Y (UTM)") +
  coord_equal() +
  theme_bw()

## 3. Conversão para 3D (rayshader)
plot_gg(
  ggobj = gg_plot,
  width = 10, 
  height = 8,
  scale = 100,           # Aumente para exagero vertical mais dramático
  shadow_intensity = 0.8,
  sunangle = 225,        # 225° = luz do noroeste
  theta = 35,            # Ângulo de rotação
  phi = 25,              # Inclinação
  windowsize = c(1200, 900)
)

## 4. Visualização interativa (HTML)
rgl::rglwidget()

Copyright

Geosaber₢2025
 
  • https://geosaber.com.br