Patron Puntual Cluster

En el siguiente código de R vamos a determinar las zonas donde se presentaron hurtos a motocicletas y establecimientos en el municipio de Santiago de Cali, con el fin de establecer posibles políticas públicas que sirvan para mitigar este problema con mayor enfoque en las áreas más afectadas, tomando como datos representativos los obtenidos para el año 2017.

Primero definimos el espacio de trabajo en el directorio personal donde se encuentran las carpetas con los shapefiles de Comunas y de Hurtos reportados por la Policía Nacional en el año 2017, verificando que quede bien direccionado.

setwd(“C:/Users/Andres/Desktop/ESPECIALIZACION UNIVALLE/1. TRATAMIENTO DE DATOS ESPACIALES - LUNES 7-10 AM/CLASE 13 - 19-05-2021 Patrones Puntuales - Cluster”)

getwd()

##Llamamos las librerías que son necesarias para el procesamiento y visualización de la información.

require(raster)

## Loading required package: raster

## Loading required package: sp

require(rasterVis)

## Loading required package: rasterVis

## Warning: package 'rasterVis' was built under R version 4.0.5

## Loading required package: terra

## Warning: package 'terra' was built under R version 4.0.5

## terra version 1.1.4

## Loading required package: lattice

## Loading required package: latticeExtra

## Warning: package 'latticeExtra' was built under R version 4.0.5

##Cargamos los shps de hurtos / comunas y generamos la tabla de atributos del shp hurtos2017 nombrándola "tabla".

hurtos=shapefile("Hurtos2017/2017_0.shp")

## Warning in OGRSpatialRef(dsn, layer, morphFromESRI = morphFromESRI, dumpSRS
## = dumpSRS, : Discarded datum Marco Geocentrico Nacional de Referencia in
## Proj4 definition: +proj=col_urban +lat_0=3.44188333333333 +lon_0=-76.5205625
## +x_0=1061900.18 +y_0=872364.63 +h_0=1000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0
## +units=m +no_defs

comunas=shapefile("Comunas/bcs_lim_comunas.shp")

## Warning in OGRSpatialRef(dsn, layer, morphFromESRI = morphFromESRI, dumpSRS
## = dumpSRS, : Discarded datum Marco Geocentrico Nacional de Referencia in
## Proj4 definition: +proj=col_urban +lat_0=3.44188333333333 +lon_0=-76.5205625
## +x_0=1061900.18 +y_0=872364.63 +h_0=1000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0
## +units=m +no_defs

tabla=hurtos@data

Se observa que la tabla de atributos del shp de hurtos contiene información como, coordenadas en el sistema de referencia geográfico (WGS84) en las columnas 6 (Longitud) y 7 (Latitud), además del tipo de hurto en la columna 11 con nombre “DELITO”, siendo estos campos necesarios para el desarrollo de la práctica.

Hacemos el intersect de la capa de hurtos con la capa de comunas para que nos muestre información exclusivamente dentro del perímetro urbano de Cali.

hurtos=intersect(hurtos,comunas)

## Loading required namespace: rgeos

plot(hurtos)

plot(comunas)
points(hurtos,col="red", cex=0.1)

En la gráfica se pueden ver los puntos con marcadores en color rojo ubicados espacialmente en los sitios donde ocurrieron los 21641 hurtos en todas las modalidades.

# Construimos una lista de nombre "densidades_hurtos" con el fin de almacenar los 12 mapas que se obtendrán mediante el siguiente iterador.

densidades_hurtos=list()

# Lo que hace el iterador o ciclo es almacenar en la posición i de 1 a 12 los hurtos de cada uno de los meses del año donde el valor en la columna "DELITO" sea igual a "HURTO VEHICULOS", imprimiendo el mapa de las comunas y sobreponiendo los puntos georreferenciados de los hurtos en esa modalidad.
# Posteriormente saca las coordenadas de los hurtos por mes haciendo un extend en las comunas para extraer los valores de coordenadas mínimos y máximos y ponerlas dentro del patrón a analizar, donde la columna 1 será la coordenada X o longitud y la columna 2 será la coordenada Y o latitud.
# Luego imprime el patrón sacando el quadrat.test para rechazar hipótesis nula o que el patrón sea aleatorio sabiendo que los hurtos generalmente manejan una tendencia de patrón agregado o tipo clúster y se hace la estimación de la densidad con un valor de sigma = 500 el cual hace que se generen los resultados de una manera adecuada.
# Por último, cada resultado se convierte en formato ráster con mascara en el shp de comunas y se imprimen los resultados que serán almacenados en la variable i hasta que el ciclo recorra los 12 meses, quedando almacenado en la lista "densidades_hurtos".

for(i in 1:12){
  hurtos_mes=hurtos[hurtos$MES==i&hurtos$DELITO=="HURTO MOTOCICLETAS",]
  plot(comunas)
  points(hurtos_mes,col="blue",cex=1)
  
  require(spatstat)
  loc_hurto=coordinates(hurtos_mes)
  extent(comunas)
  patron_hurto_mes=ppp(x = loc_hurto[,1],y = loc_hurto[,2],
                       window=owin(c(1054098 ,1068492),c(860192.1,879441.5)))
  
  plot(patron_hurto_mes)
  
  quadrat.test(patron_hurto_mes)
  densidad_hurtos=density(patron_hurto_mes,500)
  plot(densidad_hurtos)
  
  densidad_hurtos=data.frame(densidad_hurtos)
  densidad_hurtos=rasterFromXYZ(densidad_hurtos[,1:3])
  densidad_hurtos=mask(densidad_hurtos,comunas)
  plot(densidad_hurtos)
  plot(comunas,add=T)
  
  densidades_hurtos[[i]]=densidad_hurtos
}

## Loading required package: spatstat

## Warning: package 'spatstat' was built under R version 4.0.5

## Loading required package: spatstat.data

## Warning: package 'spatstat.data' was built under R version 4.0.5

## Loading required package: spatstat.geom

## Warning: package 'spatstat.geom' was built under R version 4.0.5

## spatstat.geom 2.1-0

## 
## Attaching package: 'spatstat.geom'

## The following objects are masked from 'package:terra':
## 
##     area, convexhull, coords, perimeter, rescale, rotate, shift

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

## Loading required package: spatstat.core

## Warning: package 'spatstat.core' was built under R version 4.0.5

## Loading required package: nlme

## 
## Attaching package: 'nlme'

## The following object is masked from 'package:terra':
## 
##     collapse

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

## Loading required package: rpart

## spatstat.core 2.1-2

## 
## Attaching package: 'spatstat.core'

## The following object is masked from 'package:lattice':
## 
##     panel.histogram

## Loading required package: spatstat.linnet

## Warning: package 'spatstat.linnet' was built under R version 4.0.5

## spatstat.linnet 2.1-1

## 
## spatstat 2.1-0       (nickname: 'Comedic violence') 
## For an introduction to spatstat, type 'beginner'

## Warning: data contain duplicated points

## Warning: data contain duplicated points

## Warning: data contain duplicated points

## Warning: data contain duplicated points

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

## Warning: data contain duplicated points

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

## Warning: data contain duplicated points

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

## Warning: data contain duplicated points

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

## Warning: data contain duplicated points

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

## Warning: data contain duplicated points

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

## Warning: data contain duplicated points

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

## Warning: data contain duplicated points

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

## Warning: data contain duplicated points

Finalmente creamos el Stack de mapas con los nombres de cada mes para cada mapa y lo imprimimos, tanto en la escala de los datos crudos como en porcentaje para contar con una representación más entendible y estandarizada del fenómeno.

densidades_hurtos=stack(densidades_hurtos)
names(densidades_hurtos)=month.name
plot(densidades_hurtos)

levelplot(densidades_hurtos,par.settings=BuRdTheme)

densidades_hurtos2=densidades_hurtos
densidades_hurtos2[]=densidades_hurtos[]/max(densidades_hurtos[],na.rm = TRUE)*100
levelplot(densidades_hurtos2,par.settings=BuRdTheme)

table(hurtos$DELITO)

## 
##      HURTO A CELULARES         HURTO A FINCAS       HURTO A PERSONAS 
##                   2533                      3                   4102 
##       HURTO AUTOPARTES    HURTO DE BICICLETAS HURTO ESTABLECIMIENTOS 
##                    475                    115                   1638 
##     HURTO MOTOCICLETAS       HURTO RESIDENCIA        HURTO VEHICULOS 
##                   1445                   1203                   1044

Según la distribución espacial de los 12 mapas creados para cada mes en la modalidad “HURTO MOTOCICLETAS” se puede observar que las comunas 17 y 22 tienen una afectación mínima o nula para esta problemática, además, se nota que en los primeros 3 meses del año ocurrieron eventos marcados sobre las comunas, 16 y 20 para el mes de enero, comuna 3 y 6 en el mes de febrero, comunas 3, 5, 9, 10, 14 y 21 para el mes de marzo, así como en las comunas 19, 20 y 2 para el mes de mayo, 2 para el mes de julio y 3, 8, 9 y 14 para el mes de diciembre. Al ser una tendencia no representativa durante el año sobre algunas comunas de Cali, si bien se nota una mayor afectación sobre la parte centro de Cali, se plantea para el año 2018 aumentar el pie de fuerza pública en las comunas anteriormente mencionadas, promoviendo la instalación de cámaras y alarmas que estén conectadas al CAI más cercano en los sitios donde acurren más eventos de este tipo, asimismo, se propone aumentar el patrullaje sobre las calles principales de estas comunas.

Ahora hagamos el análisis sobre otra modalidad de hurto, esta vez trabajemos sobre “HURTO ESTABLECIMIENTOS”.

densidades_hurtos3=list()

for(i in 1:12){
  hurtos_mes2=hurtos[hurtos$MES==i&hurtos$DELITO=="HURTO ESTABLECIMIENTOS",]
  plot(comunas)
  points(hurtos_mes2,col="blue",cex=1)
  
  require(spatstat)
  loc_hurto2=coordinates(hurtos_mes2)
  extent(comunas)
  patron_hurto_mes2=ppp(x = loc_hurto2[,1],y = loc_hurto2[,2],
                       window=owin(c(1054098 ,1068492),c(860192.1,879441.5)))
  
  plot(patron_hurto_mes2)
  
  quadrat.test(patron_hurto_mes2)
  densidad_hurtos3=density(patron_hurto_mes2,500)
  plot(densidad_hurtos3)
  
  densidad_hurtos3=data.frame(densidad_hurtos3)
  densidad_hurtos3=rasterFromXYZ(densidad_hurtos3[,1:3])
  densidad_hurtos3=mask(densidad_hurtos3,comunas)
  plot(densidad_hurtos3)
  plot(comunas,add=T)
  
  densidades_hurtos3[[i]]=densidad_hurtos3
}

## Warning: data contain duplicated points

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

## Warning: data contain duplicated points

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

## Warning: data contain duplicated points

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

## Warning: data contain duplicated points

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

## Warning: data contain duplicated points

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

## Warning: data contain duplicated points

## Warning: data contain duplicated points

## Warning: data contain duplicated points

## Warning: data contain duplicated points

## Warning: data contain duplicated points

## Warning: data contain duplicated points

## Warning: data contain duplicated points

densidades_hurtos3=stack(densidades_hurtos3)
names(densidades_hurtos3)=month.name
plot(densidades_hurtos3)

levelplot(densidades_hurtos3,par.settings=BuRdTheme)

densidades_hurtos4=densidades_hurtos3
densidades_hurtos4[]=densidades_hurtos3[]/max(densidades_hurtos3[],na.rm = TRUE)*100
levelplot(densidades_hurtos4,par.settings=BuRdTheme)

Revisando los mapas resultantes de este análisis para los 12 meses del año, y a diferencia del resultado anterior (HURTOS MOTOCICLETAS) que no mostraba un patrón marcado sobre alguna comuna en específico, para este tipo de hurto, se presenta mayor afectación sobre las comunas 3 y 9 (centro de Cali) donde indiscutiblemente se debe tener una mayor y mejor comunicación entre las personas que atienden los establecimientos y la policía del cuadrante más cercano, si bien se espera que en todos los establecimientos existan cámaras de vigilancia, se requiere de una herramienta de alerta inmediata que posibilite el rápido apoyo de las autoridades a estos establecimientos, como alarmas silenciosas conectadas en tiempo real a los CAI o estaciones de policía.