Análisis de Información Geográfica y Espacial

Maestría en Ciencia de Datos - Javeriana Cali

Caso de Estudio

Se realiza un análisis de patrones puntales tomando como referencia la región centro occidente del país, que cubre departamentos como Huila, Tolima, Valle del Cauca, Eje Cafetero, Cundinamarca, Bogotá y parte de Meta y Antioquia respectivamente.

Se segmenta esta región del país para analizar el comportamiento de la actividad sísmica, tomando como referencia de tiempo entre enero de 2016 y septiembre de 2023.

Desarrollo del Caso

Importación de Librerías y cargue de la base de datos

Se realizar la importación de las librerías requeridas para el análisis de patrones puntuales de los sismos en Colombia.

# importación de Librerías

require(spatstat)
require(rgdal)
require(maptools)
require(sp)
require(raster)
require(stpp)
require(plot3D)
require(KernSmooth)
require(leaflet)
library(sf) # Development version in the ppp branch
library(tmap)
library(here)
library(readr)

Posteriormente, se realiza el cargue de la base de datos del National Earthquake Hazards Reduction Program (NEHRP).

# Importación de la base de datos

sismos <- read_csv("C:/Users/User/Downloads/sismo.csv")

Describir los sismos

A partir de la librería Leaflet se construye un mapa para ubicar cada uno de los sucesos sísmicos presentados, para empezar a identificar patrones en la ocurrencia de estos.

leaflet() %>% addTiles() %>% 
  addCircleMarkers(
    lng = sismos$longitude,
    lat = sismos$latitude, 
    label = paste0("Ubicación: ",as.character(sismos$place),
                   " | Sismo de magnitud: ",as.character(sismos$mag),
                   " | Latitud: ",as.numeric(sismos$latitude),
                   " | Longitud: ", as.numeric(sismos$longitude))
  )

Convertir en patrón puntual

Se construye el patrón puntual de los sismos a partir de la ventana seleccionada para la obtención de los datos, logrando identificar una concentración entre los departamentos de Meta y Cundinamarca, así como entre Chocó y el eje cafetero.

zona=owin(xrange = c(-77.355, -73.114),yrange = c(2.089, 5.736))
patron_sismos=ppp(x = sismos$longitude,y = sismos$latitude,window =zona )
plot(patron_sismos)

Identificar el tipo de patrón

Con la gráfica correspondiente a la aplicación del Kest se puede determinar que el patrón es muy agregado ya que este se encuentra por encima de la línea azul.

plot(Kest(patron_sismos))

El patrón es muy agregado porque está por encima de la linea azul.

Predicción de la probabilidad

#prueba chi-cuadrado
qc.loc <- quadratcount(patron_sismos, nx=5, ny=5)
plot(patron_sismos, pch=3, cex=0.6)
plot(qc.loc, add=T, textargs = list(col='red'))

Con la predicción anterior, se aprecia dos puntos en los cuales se generan unos patrones de ocurrencia de estos fenómenos sismicos.

bw.diggle(patron_sismos)
##      sigma 
## 0.09367295
densidad_sismos=density(patron_sismos)
plot(densidad_sismos)

##intensidad proceso aleatorio##
#intensidad espacial
intensidad_sismos = density.ppp(patron_sismos,diggle=TRUE,sigma=bw.diggle(patron_sismos))
intensidad_sismos1 = density.ppp(patron_sismos,diggle=TRUE,sigma=bw.ppl(patron_sismos))
intensidad_sismos2 = density.ppp(patron_sismos,diggle=TRUE,sigma=bw.ppl(patron_sismos), kernel="epanechnikov")
intensidad_sismos3 = density.ppp(patron_sismos,diggle=TRUE,sigma=bw.ppl(patron_sismos), kernel="quartic")
intensidad_sismos4 = density.ppp(patron_sismos,diggle=TRUE,sigma=bw.scott(patron_sismos))

intensidad_sismos5=densityAdaptiveKernel.ppp(patron_sismos)

A continuación se construye la función de intensisdad a partir de distintos métodos

# Método Diggle - Usando validación cruzada para la selección del ancho de banda

par(mfrow=c(1,1))
plot(intensidad_sismos,xlab="x",ylab="y",main=paste("Modelo Gaussiano sigma=",round(attr(intensidad_sismos,"sigma"),5)))
contour(density(patron_sismos),add=TRUE)

# Método Verosimilitud - Usando validación cruzada para la selección del ancho de banda


plot(intensidad_sismos1,xlab="x",ylab="y",main=paste("Modelo Gaussiano sigma=",round(attr(intensidad_sismos1,"sigma"),5)))
contour(density(patron_sismos),add=TRUE)

# Método Epanechnikov - Usando validación cruzada para la selección del ancho de banda


plot(intensidad_sismos2,xlab="x",ylab="y",main=paste("Modelo epanechnikov sigma=",round(attr(intensidad_sismos2,"sigma"),5)))
contour(density(patron_sismos),add=TRUE)

# Método Quartic - Usando validación cruzada para la selección del ancho de banda

plot(intensidad_sismos3,xlab="x",ylab="y",main=paste("Modelo método quartic sigma=",round(attr(intensidad_sismos3,"sigma")),5))
contour(density(patron_sismos),add=TRUE)

# Método Scott - Usando validación cruzada para la selección del ancho de banda

plot(intensidad_sismos4,xlab="x",ylab="y",main=paste("Modelo método scott sigmas=",round(bw.scott(patron_sismos),5)))
contour(density(patron_sismos),add=TRUE)

# Método Kernel Adaptativo

plot(intensidad_sismos5,xlab="x",ylab="y",main="Kernel Adaptativo")
contour(density(patron_sismos),add=TRUE)

# Método Scott Sigmas- Usando validación cruzada para la selección del ancho de banda

plot(intensidad_sismos4,xlab="x",ylab="y",main=paste("Modelo método scott sigmas=",round(bw.scott(patron_sismos),5)))
contour(density(patron_sismos),add=TRUE)
points(patron_sismos,col="red",pch=16)

##crear un raster con la densidad
tabla=data.frame(densidad_sismos)[,1:3]
require(raster)
raster_sismos=rasterFromXYZ(tabla)
# Sistema de Referencia

crs(raster_sismos) = "+init=EPSG:4326"

pos=which(raster_sismos[]<0.6)
raster_sismos[][pos]=NA
plot(raster_sismos)

# Mapa

leaflet() %>% addTiles() %>% addRasterImage(raster_sismos,opacity = 0.6)

Finalmente, se pudo determinar que las regiones incluidas dentro de la ventana de trabajo, los territorios que registran los mayores patrones de ocurrencia de sismos están hacia el oriente entre los departamentos de Meta, Huila y la zona rural de Bogotá, mientras que hacia el occidente se observa un patrón en los límites entre Chocó y los departamentos de Valle del Cauca, Risaralda y Quindío.