Patrones puntuales

Patrones puntuales

# PATRONES PUNTUALES SISMOS

require(spatstat)

## Loading required package: spatstat

## Loading required package: spatstat.data

## Loading required package: nlme

## Loading required package: rpart

## 
## spatstat 1.64-1       (nickname: 'Help you I can, yes!') 
## For an introduction to spatstat, type 'beginner'

## 
## Note: spatstat version 1.64-1 is out of date by more than 6 months; we recommend upgrading to the latest version.

require(readxl) #cargar la informacion de sismos (query.csv) desde import dataset - from text (readr) pq es .csv

## Loading required package: readxl

fresmo_1_Da <- read.csv("D:/ESPECIALIZACION/SEMESTRE_1/1. Tratamiento de datos/Clase12_patrones puntuales/datos.csv")
require(raster)

## Loading required package: raster

## Loading required package: sp

## 
## Attaching package: 'raster'

## The following objects are masked from 'package:spatstat':
## 
##     area, rotate, shift

## The following object is masked from 'package:nlme':
## 
##     getData

global=shapefile("D:/ESPECIALIZACION/SEMESTRE_1/1. Tratamiento de datos/Clase12_patrones puntuales/shape_global/g2008_0.shp")


##--- 1. Del shape global extrae solo Colombia

sismos= fresmo_1_Da
par(mfrow=c(1,1))
col=global[global$ADM0_NAME=="Colombia",0] #sacar del shape la zona de Colombia
plot(col)
points(sismos[,3:2],col="red")             #3:2 coord longu,latitud de los datos de sismos

##--- 2. Extrae solo los datos de sismos contenidos en Colombia 

require(sp)
require(rgeos)

## Loading required package: rgeos

## rgeos version: 0.5-5, (SVN revision 640)
##  GEOS runtime version: 3.8.0-CAPI-1.13.1 
##  Linking to sp version: 1.4-4 
##  Polygon checking: TRUE

coords=data.frame(sismos[,3:2])
sismos_sp=SpatialPoints(coords,proj4string = crs(col))
valids=extract(col,sismos_sp)

## Warning in proj4string(x): CRS object has comment, which is lost in output

## Warning in proj4string(y): CRS object has comment, which is lost in output

sismos_col=sismos_sp[which(valids$poly.ID==1),]
plot(col)
plot(sismos_col,add=T,col="red")

##--- 3. Crear un patron puntual para el área de estudio "Colombia"

coords_col=coordinates(sismos_col)
extent(col)

## class      : Extent 
## xmin       : -79.04723 
## xmax       : -66.86983 
## ymin       : -4.230484 
## ymax       : 12.46331

require(spatstat)
patron_sismos=ppp(x=coords_col[,1],y=coords_col[,2],window=owin(c(-79.04723 ,-66.86983),c(-4.230484,12.46331)))
plot(patron_sismos)

##--- 4. Encontrar el tipo de patron - Chi-Cuadrado, p-value 

quadrat.test(patron_sismos)

## Warning: Some expected counts are small; chi^2 approximation may be inaccurate

## 
##  Chi-squared test of CSR using quadrat counts
## 
## data:  patron_sismos
## X2 = 389.34, df = 24, p-value < 2.2e-16
## alternative hypothesis: two.sided
## 
## Quadrats: 5 by 5 grid of tiles

plot(Kest(patron_sismos), main = "")

#--- Dado que p-value es cercano a 0 lo que indica los sismos en el pais a nivel espacial no ocurren de forma aleotoria.
#--- Hay lugares que presentan mayor probabilidad o riesgo de presenar un sismo posiblemte por la gemorfologia de la zona. 
#--- el patron con los datos de sisimos en Colombia muestran un patron tipo Closter,
#--- donde la frecuencia observada esta muy por encima de la frecuencia aleatoria 


##--- 5. Estimar la intensidad

plot(density(patron_sismos,1),main="Intensidad estimada")
contour(density(patron_sismos,1),add=TRUE)

intensidad=density(patron_sismos,1)
intensidad=data.frame(intensidad)
intensidad_map=rasterFromXYZ(xyz =intensidad[,1:3])

plot(intensidad_map, col=colorRampPalette(c("white", "yellow","orange","red"))(100))
plot(col,add=T)
points(coords_col,col="red")

##--- 6. Generar una máscara para recortar la intensidad al limite de Colombia

intensidad_col=mask(intensidad_map,col)
plot(intensidad_col, col=colorRampPalette(c("white", "yellow","orange","red"))(100))
plot(col,add=T)
points(coords_col,col="black",pch=16,cex=0.3)