#Análisis Geoestadístico del Plomo en el Valle del Meuse
###El siguiente flujo de trabajo implementa un análisis exploratorio y geoestadístico de concentraciones de plomo (Pb) en el valle del río Meuse. Se emplean objetos espaciales tipo sf, análisis variográfico y técnicas básicas de estadística espacial descriptiva.
##1. Importación y Preparación de Datos Espaciales
setwd("C:/Users/vicen/Downloads/RMarkDown")
library(gstat)
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 4.0.0 ✔ tibble 3.2.1
## ✔ lubridate 1.9.4 ✔ tidyr 1.3.1
## ✔ purrr 1.0.4
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(sf)
## Linking to GEOS 3.13.0, GDAL 3.10.1, PROJ 9.5.1; sf_use_s2() is TRUE
library(terra)
## terra 1.8.42
##
## Adjuntando el paquete: 'terra'
##
## The following object is masked from 'package:tidyr':
##
## extract
meuse <- st_read("meuse.shp")
## Reading layer `meuse' from data source
## `C:\Users\vicen\Downloads\RMarkDown\meuse.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 1 field
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 178440 ymin: 329600 xmax: 181560 ymax: 333760
## Projected CRS: Double_Stereographic
st_crs(meuse)
## Coordinate Reference System:
## User input: Double_Stereographic
## wkt:
## PROJCRS["Double_Stereographic",
## BASEGEOGCRS["GCS_Bessel 1841",
## DATUM["unknown",
## ELLIPSOID["bessel",6377397.155,299.1528128,
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["Degree",0.0174532925199433]]],
## CONVERSION["unnamed",
## METHOD["Oblique Stereographic",
## ID["EPSG",9809]],
## PARAMETER["Latitude of natural origin",52.1561605555556,
## ANGLEUNIT["Degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",5.38763888888889,
## ANGLEUNIT["Degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["Scale factor at natural origin",0.9999079,
## SCALEUNIT["unity",1],
## ID["EPSG",8805]],
## PARAMETER["False easting",155000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",463000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]]],
## CS[Cartesian,2],
## AXIS["(E)",east,
## ORDER[1],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]],
## AXIS["(N)",north,
## ORDER[2],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]]]
plomo <- read.csv("datos_plomo.csv") %>%
st_as_sf(coords = c("x","y"), crs = 9001)
st_crs(plomo) <- st_crs(plomo)
###Se carga la capa vectorial meuse.shp, que representa la geometría del valle.
###Se importan los puntos de muestreo de plomo desde un CSV y se convierten a un objeto sf utilizando sus coordenadas.
###Se establece el sistema de referencia espacial (CRS 9001, generalmente asociado a proyecciones planas métricas).
##2. Visualización Espacial del Dominio de Estudio
plot(meuse$geometry, main="plomo - valle del meuse")
plot(plomo$geometry, col="red", add = TRUE)
###Se genera un mapa base con los límites del valle y se superponen los puntos de muestreo de plomo, permitiendo visualizar la distribución espacial de la red de observaciones.
##3. Exploración Estadística Univariada
copia_plomo <- plomo %>% st_drop_geometry() %>% dplyr::select(-c(X))
hist(copia_plomo$plomo, xlab="Pl")
summary(copia_plomo)
## plomo
## Min. : 37.0
## 1st Qu.: 72.5
## Median :123.0
## Mean :153.4
## 3rd Qu.:207.0
## Max. :654.0
boxplot(copia_plomo$plomo, main="Boxplot - Plomo en el valle meuse")
##Las geometrías se separan para obtener una tabla puramente atributiva.
###Se realizan:
### Histograma → Distribución de frecuencias del plomo.
###Medidas descriptivas → Media, mediana, cuartiles y valores extremos.
###Diagrama de caja → Identificación de valores atípicos y dispersión general.
##4. Análisis de Continuidad Espacial mediante Cloud de Semivarianzas
hscat(formula = plomo ~ 1, data = plomo,
breaks = c(20,30,90,120,200,400,700,1000,1500,2500,5000))
###El gráfico h-scatterplot evalúa la relación espacial entre pares de puntos a intervalos crecientes de distancia, proporcionando evidencia preliminar de autocorrelación espacial.
##5. Tendencias Espaciales en los Ejes Cardinales
par(mfrow=c(2,1))
df_plomo <- read.csv("datos_plomo.csv", header = TRUE, sep=",")
plot(df_plomo$x, df_plomo$plomo, xlab="oeste-este", ylab="plomo")
abline(lm(df_plomo$plomo ~ df_plomo$x), col="red")
plot(df_plomo$y, df_plomo$plomo, xlab="sur-norte", ylab="plomo conc")
abline(lm(df_plomo$plomo ~ df_plomo$y), col="red")
##Se examina la existencia de tendencia direccional (gradient spatial) de la concentración de plomo:
### Eje Oeste–Este (X)
###Eje Sur–Norte (Y)
###La recta de regresión permite identificar si hay pendiente significativa que sugiera un componente estructural antes de ajustar el variograma.
##6. Construcción del Mapa Variográfico
gm_map <- variogram(plomo ~ 1, data = plomo,
cutoff = 1000, width = 100, map = TRUE)
plot(gm_map, col.regions = viridis::viridis(100),
main = "Mapa variografico de plomo en valle Meuse",
xlab = "Variacion eje X",
ylab = "Variacion eje Y")
###Se genera un variograma direccional en formato de mapa, representando la estructura de dependencias espaciales:
### Cutoff = 1000 m: distancia máxima considerada.
###Width = 100 m: ancho de los bins espaciales.
###El resultado es un mapa que muestra anisotropías o patrones direccionales en la semivarianza.
###Resumen Técnico
###Este flujo corresponde a un análisis clásico en geoestadística:
### Preparación y limpieza de datos espaciales
###Exploración estadística y visualización temática
###Evaluación preliminar de autocorrelación espacial
###Detección de tendencias regionales
###Cálculo de un variograma empírico en formato direccional
###Este pipeline es la base para la modelación variográfica y posterior interpolación mediante kriging.