# 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)
<- st_as_sf(
meuse_sf
meuse,coords = c("x", "y"),
crs = 28992 # EPSG 28992 = RD New (Holanda)
)<- st_as_sf(meuse.grid, coords = c("x", "y"), crs = 28992)
meuse_grid_sf
# 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")
Visualização 3D com Rayshader
Curso de Análise Espacial e Interpolação de Dados no R
1. Preparação dos dados
2. Visualização 3D Básica com WebGL
library(plotly)
# Dados simplificados
<- data.frame(
plot_data 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
<- SpatialPointsDataFrame(meuse[, c("x", "y")],
meuse_sp data = meuse,
proj4string = CRS("+init=EPSG:28992"))
<- SpatialPixels(SpatialPoints(meuse.grid[, c("x", "y")],
meuse_grid_sp proj4string = CRS("+init=EPSG:28992")))
# 2. Análise de variograma
<- variogram(elev ~ 1, meuse_sp)
variogram_model <- fit.variogram(variogram_model, vgm(1, "Gau", 800, 0.1))
vgm_model
# 3. Execute a krigagem
<- krige(elev ~ 1, meuse_sp, newdata = meuse_grid_sp, model = vgm_model) krige_result
[using ordinary kriging]
plot(krige_result)
# 4. Converter para terra e selecionar apenas a predição
<- rast(krige_result)[[1]] # [[1]] seleciona a primeira camada
elev_krige names(elev_krige) <- "elevation"
4. Cálculo do Gradiente Hidráulico
<- elev_krige[[1]]
elev_layer # Calcular componentes do gradiente (direção do fluxo)
<- terrain(elev_layer, v = c("slope", "aspect"), unit = "radians")
gradient
# Componentes do vetor gradiente (PERPENDICULAR às isolinhas)
<- sin(gradient$aspect) * gradient$slope # Componente x
dx_rast <- cos(gradient$aspect) * gradient$slope # Componente y dy_rast
5. Geração das Setas de Fluxo
<- st_make_grid(meuse_sf, cellsize = 350, what = "centers") %>%
arrow_grid st_as_sf() %>%
vect()
# Extração e normalização dos vetores
<- data.frame(
arrow_data 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)
<- as.data.frame(elev_krige, xy = TRUE) # Converter para dataframe
dados_raster
## 2. Plot 2D (ggplot2) com elementos personalizados
<- ggplot() +
gg_plot
# 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)
::rglwidget() rgl