7 Geometría 3D con R y PostGIS

A partir del MDT de nuestro municipio, vamos a obtener una muestra de puntos en 3D, a almacenarlos en PostGIS y a representarlos en R.

Mediante el fragmento de código siguiente, obtenemos una muestra de puntos en 3D.

file <- 'datos/Lanjaron/p03out/mdt-lanjaron-bbox-jsamos.tif'
mdt <- terra::rast(file)

# Extraer puntos y elevaciones del MDT
puntos <- as.data.frame(mdt, xy = TRUE, na.rm = TRUE)
colnames(puntos) <- c("x", "y", "z")

# Seleccionar una muestra aleatoria
muestra <- sample(1:nrow(puntos), size = 1000, replace = FALSE)
puntos <- puntos[muestra, ]

A continuación, podemos generar una capa a partir de las coordenadas de los puntos, mediante el código siguiente.

puntos$geom <- paste0("POINT Z (", puntos$x, " ", puntos$y, " ", puntos$z, ")")
puntos$id <- 1:nrow(puntos)
puntos <- puntos[, c('id', 'geom')]
sf_puntos <- sf::st_as_sf(puntos, wkt = "geom", crs = sf::st_crs(mdt))

Los podemos almacenar directamente en PostGIS, mediante el código siguiente.

con <- RPostgres::dbConnect(
  RPostgres::Postgres(),
  dbname = 'lanjaron_jsamos',
  host = 'localhost',
  port = '5432',
  user = 'postgres',
  password = 'postgres'
)

sf::st_write(obj = sf_puntos,
             dsn = con,
             layer = "p3d_lanjaron_bbox_jsamos")

RPostgres::dbDisconnect(con)

También podemos representar los puntos desde R, mediante el código siguiente.

xyz <- sf::st_coordinates(sf_puntos)
rgl::plot3d(xyz[, 'X'], xyz[, 'Y'], xyz[, 'Z'],
       col = "blue", size = 1.5, type = "s",
       xlab = "X", ylab = "Y", zlab = "Z")

Figura 7.1: Resultado de 3D de muestra de puntos desde R.

En la figura 7.1 se muestra el resultado.