# Limpiar la informacion almacenada en la memoria
rm(list=ls())

# Cargar los paquetes
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 7 months; we recommend upgrading to the latest version.
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
require(sp)

library(readr)
g2008_0 <- read_csv("C:/Users/linam/Desktop/g2008_0.shp")
## Registered S3 method overwritten by 'cli':
##   method     from    
##   print.boxx spatstat
## 
## -- Column specification --------------------------------------------------------
## cols(
##   col_character()
## )
## Warning: 106872 parsing failures.
## row col  expected        actual                                 file
##   1               embedded null 'C:/Users/linam/Desktop/g2008_0.shp'
##   1  -- 1 columns 3 columns     'C:/Users/linam/Desktop/g2008_0.shp'
##   2               embedded null 'C:/Users/linam/Desktop/g2008_0.shp'
##   3               embedded null 'C:/Users/linam/Desktop/g2008_0.shp'
##   4               embedded null 'C:/Users/linam/Desktop/g2008_0.shp'
## ... ... ......... ............. ....................................
## See problems(...) for more details.
View(g2008_0)


library(readr)
query <- read_csv("C:/Users/linam/Desktop/Doctorado-20201031T024248Z-001/Doctorado/Semestre V/Análisis Espacial de Datos/query.csv")
## 
## -- Column specification --------------------------------------------------------
## cols(
##   .default = col_double(),
##   time = col_datetime(format = ""),
##   magType = col_character(),
##   nst = col_logical(),
##   net = col_character(),
##   id = col_character(),
##   updated = col_datetime(format = ""),
##   place = col_character(),
##   type = col_character(),
##   status = col_character(),
##   locationSource = col_character(),
##   magSource = col_character()
## )
## i Use `spec()` for the full column specifications.
View(query)

#----------------------------------------------------------------------

par(mfrow=c(2,3)) #Graficos

# Datos de inhibicion que provienen de una muestra de c??lulas biol??gicas en una regi??n rectangular.
data(cells)
plot(cells, main = "Celulas",pch=16)

# Datos de aleatoriedad completa que provienen de las ubicaciones de pinos en Jap??n en una regi??n rectangular.
data(japanesepines) 
plot(japanesepines, main = "Pinos",pch=16)

# Datos de agregacion que provienen de ubicaciones de Secoyas (??rboles) en una regi??n rectangular.
data(redwood)
plot(redwood, main = "Secoyas",pch=16)

##Test de Cuadrantes
plot(quadratcount(cells,nx = 5,ny = 5),cex=2)
plot(quadratcount(japanesepines),cex=2)
plot(quadratcount(redwood),cex=2)

quadrat.test(cells)
## Warning: Some expected counts are small; chi^2 approximation may be inaccurate
## 
##  Chi-squared test of CSR using quadrat counts
## 
## data:  cells
## X2 = 9.1905, df = 24, p-value = 0.005679
## alternative hypothesis: two.sided
## 
## Quadrats: 5 by 5 grid of tiles
quadrat.test(japanesepines)
## Warning: Some expected counts are small; chi^2 approximation may be inaccurate
## 
##  Chi-squared test of CSR using quadrat counts
## 
## data:  japanesepines
## X2 = 16.923, df = 24, p-value = 0.2962
## alternative hypothesis: two.sided
## 
## Quadrats: 5 by 5 grid of tiles
quadrat.test(redwood)
## Warning: Some expected counts are small; chi^2 approximation may be inaccurate
## 
##  Chi-squared test of CSR using quadrat counts
## 
## data:  redwood
## X2 = 64.613, df = 24, p-value = 2.774e-05
## alternative hypothesis: two.sided
## 
## Quadrats: 5 by 5 grid of tiles
# K para los datos de las celulas
plot(Kest(cells), main = "K-estimada (Celulas)")
# K para los datos de los pinos
plot(Kest(japanesepines), main = "K-estimada (Pinos)")
# K para los datos de las secoyas
plot(Kest(redwood), main = "K-estimada (Secoyas)")

##Estimar la funcion de intensidad (Kernel)

par(mfrow=c(2,2))

plot(density(cells,sigma=5), main = "Densidad K=5 (Secoyas)")
plot(density(cells,sigma = 0.1), main = "Densidad K=0.1 (Secoyas)")
plot(density(cells,sigma = 0.05), main = "Densidad K=0.05 (Secoyas)")
plot(density(cells,sigma = 0.005), main = "Densidad K=0.01 (Secoyas)")

par(mfrow=c(1,3))
plot(redwood)
plot(density(redwood,sigma = 0.1), main = "Densidad K=0.1 (Secoyas)")
plot(density(redwood,sigma = 0.05), main = "Densidad K=0.05 (Secoyas)")

##Patrones Puntuales Sismos
require(raster)
global=shapefile("C:/Users/linam/Desktop/g2008_0.shp")
sismos=query
par(mfrow=c(1,1))
col=global[global$ADM0_NAME=="Colombia",0]
plot(col)
points(sismos[,3:2],col="red")

## Contorno
require(sp)
coords=data.frame(sismos[,3:2])
sismos_sp=SpatialPoints(coords,proj4string = crs(col))
valids=extract(col,sismos_sp)
## Loading required namespace: rgeos
## 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")

##Para cambiar sistema de referencia

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)

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

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")

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)

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.