Un recorrido por el tratamiento y análisis de datos espaciales con R
Ricardo Robayo Martínez
Sociedad Ecuatoriana de Estadística - R Weekend
1. Introducción
¿Qué es un dato espacial?
Un dato espacial es aquel que dispone de un componente de georreferenciación; característica o un atributo que lo hace único y por el cual podemos ubicarlo en algún lugar de la Tierra.
Dada la particularidad anterior, los datos espaciales son únicos y la estadística destina una rama para su estudio, la cual además de solventar el análisis también debe aboradar la forma con la que se puede representar y transformar la información goegráfica.
Originalmente R no disponia de librerías para procesar datos espaciales, no obstante, en las dos últimas décadas progresivamente expertos y entunciastas en la estadística espacial han trabajado para que R pueda procesar, visualizar y modelar este tipo de datos, puesto que estos necesitan de una gran capacidad computacional.
¿Dónde se encuentran los datos espaciales?
En las dos últimas décadas, cada minuto se generan miles de millones de datos, los cuales paulatinamente han incorporado atributos geográficos. Por ejemplo, cuando hacemos una compra, posteamos un tweet, hacemos una denuncia delictiva o si se dispone de sensores que miden el clima; todas estas acciones en su mayoría generan información que está sujeta a una dirección o a un par de coordenadas geográficas (latitud y longitud), a esto se le llama ubicación geográfica.
¿Cuál es su uso?
Por medio de la información georreferenciada se puede conocer el lugar específico o parcial en el que sucedio un evento de interés, y así ofrecer campañas publicitarias, desarrollar operativos policiales, conocer el clima de una área determinada para cosechar o sembrar, entre otras cosas más; por esto, la información geográfica es una herramienta potente para la toma de decisiones tanto en el sector público como privado.
¿En qué áreas se puede aplicar?
En Ecuador, los servicios financieros, la agricultura, la pesca, la ecología, la salud, la climatología y la criminología, entre otros muchos sectores más han y están aprovechado la información georreferenciada.
¿Por qué utilizar R para el manejo de datos espaciales?
Historicamente el proceso de datos espaciales estaba cubierto por software de paga y hardware sofisticados, mismos que no eran accesibles facilmente. Sin embargo, en la actualidad es posible descargar paquetes espaciales de alto rendimiento y ejecutarlo sobre una computadora que tenga al menos menos 4 Gb de RAM.
El primer software que permitió a gran esacala que el análisis de datos espaciales sea accecible fue QGIS; mismo que es un Sistema de Información Georreferenciada (SIG) de forma nativa y una interfaz de gráficas de usuario (GUI), lo cual desalienta la reprductividad en los trabajos de ciencia de datos o estadísticos, aunque esto puede ser solventado por medio de integración de código por medio de R.
R a diferencia de los SIG nativos, que prefieren un GUI, este software trabaja sobre un Interfaz de línea de comandos (CLI), con lo cual se permite procesos automáticos, reproducibles y personalizables.
Además de R, otros software como C++, Python o Java pueden manejar, procesar, visualizar y modelar datos espaciales, no obstante, R y Python permiten con mayor facilidad integrarse con ArcGis, QGIS e incluso con herramientas de Google.
Según Lovelace, Nowosad y Muenchow (2021), R en particular es superior a otros software dado que tienen algunas librerías únicas en su clase para modelar, además de la versatilidad para desarrollar estudios estadísticos. Sin embargo, si se desea mejorar la rapidez de compilación es aconcejable utilizar C++ o en su defecto se podría utilizar librerías que permiten acceder a C++ desde R (Rccp), y en el caso de aprendizaje profundo con datos espaciales puede que Python sea superior.
Finalmente, el utilizar R, tal como otros lenguajes de programación en general conlleva una gran variedad de ventajas pero también de desventajas. En el caso de R, la principal desventaja es que la curva de aprendizaje es pronunciada, no obstante, esto es opacado y considerado como irrelevante dado que permite desarrollar análisis y visualizaciones personalizables, transparentes y reproducibles.
2. Una breve historia de la estadística espacial
¡La estadística espacial está cerca!
¿Quién es Jonh Snow?
Es considerado como el padre de la epidemiología, ya que fue el primer médico y estadístico que utilizó el concepto de clusters espaciales al analizar datos del cólera en Inglaterra, puesto que todos los que tomaban agua cerca de una bomba general tenian mayor propención a tener cólera.
Waldo R. Tobler (1969) dice que: “Todos los lugares están relacionados, pero los lugares cercanos están más relacionados que los lugares lejanos.”
3. Datos espaciales
Los fenómenos espaciales generalmente pueden ser representados por medio de objetos discretos con límites claros (ej. mapas), o en fenómenos contínuos que se puede observar en todas partes pero con límites naturales (ej. imágenes satelitales).
En el caso de los objetos espaciales discretos, estos pueden referirse a un río, un país, una ciudad o un sitio de investigación cualquiera.
Este tipo de objetos espaciales se represetan por medio de datos vectoriales; los cuales al unirse entre varios constituyen geometrías, mismas que pueden representar las fronteras de los países del mundo, las provincias del Ecuador o puntos donde suceden determinados sucesos, que para caracterizarlos hacen refrencia a nombres de los paises, representan el tamaño de alguna población, la edad promedio, el índice de desigualdad; valores que son caracteristicas y que se les demoninan “atributos”.
4. Modelo de datos vectoriales
Representan al mundo por medio de puntos, líneas y polígonos; geometrías que están construidos por pares de coordenadas (x, y).
Puntos: son los objetos más simples de datos espaciales, y tienen un par de coordenadas (x,y), además de n atributos asociados.
Un punto puede represetar el lugar donde se realizó una compra y contener otros atributos más como: la hora, el día, el monto o el sector de una compra.
Adicionalmente, se puede agrupar a todos los puntos como conjunto de puntos, es decir, multipuntos, por ejemplo, se puede unir las cafeterías representadas por puntos, pasando a obtener un poligono de cafeterías, obteniendo información agregada.
Líneas: son las geometrías que resultan de la unión secuancial u ordenada de puntos o coordenadas (nodos), al unir varias de estas se obtienen polilíneas; con las cuales es posible representar uno o varios rios o carretera.
Polígonos: hace referencia al conjunto de polilíneas cerradas dado que ahora el último punto de la línea se conecta con el primer punto. Esta geometría es considerada como la más compleja compleja de todas, dado que es la unión de las anteriores. En el caso de que se dispongan de varios poligonos, estos pueden considerarse también como una geometría, lo cual se denomina multipoligonos, con lo que podemos construir el mapas.
5. Proyección de poligonos en el mapas
## Creamos el poligono donde y es la latitud y x es la longitud
POLIGONO <- data.frame(
x = c(-78.493957, -78.492162, -78.488182, -78.485296,-78.4865284, -78.4875284, -78.490855, -78.4921097),
y = c(-0.211751, -0.206055, -0.205250, -0.208694, -0.2093347, -0.2093347, -0.210689, -0.2101194))
## utilizamos librerías espaciales para contruir un poligono espacial (librería: sp)
PuntosPoligono1Poligono <- Polygon(POLIGONO)
## Construimos una lista de poligonos para generar una geometría
PuntosPoligono1PoligonoLista <- Polygons(list(PuntosPoligono1Poligono), "Poligono aletorio de HURTOS")
## Llegamos ha un SpatialPolygons
SpatialData <- SpatialPolygons(list(PuntosPoligono1PoligonoLista))
## Utilizamos la librería sf para estructurar los datos espaciales como un DataFrame
ShpDataFrame <- st_as_sf(SpatialData)
## Utilizamos la librería leaflet para proyectar los puntos en mapas
leaflet() %>%
addTiles() %>%
addPolygons(data = ShpDataFrame,
color = "green",
fillOpacity = 0.08) %>%
addMarkers(lng= -78.488945,
lat= -0.208713,
popup= "Se venden Poliburguers")6. ¿Qué hace que los datos vectoriales sean un objeto espacial?
Para que un dato vectorial pueda ser considerado como espacial necesita tener un Sistema de Referencia de Coordenadas (CRS). El CRS nos sirve para representar un mundo que es tridimensional en un espacio bidimencional.
Sin embargo, aún con esta transaformación, existen algunas propiedades de la superficie 3D que la 2D las pierde. Por ejemplo, el área, la distancia, la forma y la dirección se distorsionan al crear un mapa 2D.
Listado de sistemas de coordenadas
SRID - Spatial Reference System Identifier - Sistema de Referencia Espacial
EPSG - European Petroleum Survey Group - Grupo Europeo de Estudios del Petróleo
- Utilizaremos el WGS 84 - World Geodetic System 1984 - Sistema geodésico mundial 1984 - EPSG 4326.
https://epsg.io/4326 https://epsg.io/32617 https://epsg.org/ https://bookdown.org/martinmontaneb/CienciaDeDatosParaCuriosos/datos-espaciales-en-r.html https://spatialreference.org/ref/epsg/wgs-84/ Más detalle acerca de los CRS https://mappinggis.com/2016/04/los-codigos-epsg-srid-vinculacion-postgis/
7. Caso de práctico
Hablemos sobre la realidad delictiva del Ecuador en el 2017.
Ilustraciones de Ecuador por provincias y parroquias
repositorio_provincias <- "C:/Users/Usuario/Desktop/SEE/Curso SEE 1ra edicion/Curso Espacial 1ra edicion/Material/Shape File/Provincias_Cantones_Parroquias"
setwd(repositorio_provincias)
## librería para cargar y tratar datos espaciales de la misma forma que se hace con los dataframe
# library(sf)
##############################################
## Ilustraciones de Ecuador por provincias ###
##############################################
### Desagregacion espacial usada por el INEC a nivel de Provincias
Provincias <- read_sf("nxprovincias.shp")
head(Provincias)## Simple feature collection with 6 features and 7 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: 637401.6 ymin: 9598532 xmax: 886203.7 ymax: 10132520
## projected CRS: WGS 84 / UTM zone 17S
## # A tibble: 6 x 8
## DPA_PROVIN DPA_DESPRO DPA_VALOR DPA_ANIO REI_CODIGO REN_CODIGO PEE_CODIGO
## <chr> <chr> <int> <chr> <chr> <chr> <chr>
## 1 01 AZUAY 0 2012 05 01 593
## 2 02 BOLIVAR 0 2012 02 01 593
## 3 03 CAÑAR 0 2012 05 01 593
## 4 04 CARCHI 0 2012 04 01 593
## 5 05 COTOPAXI 0 2012 02 01 593
## 6 06 CHIMBORAZO 0 2012 02 01 593
## # ... with 1 more variable: geometry <MULTIPOLYGON [m]>### Modificacmos las proyecciones del shp a 4326 o WGS84, mismas que utiliza el GPS y nos permite proyectar los puntos tan solo utilizando geom_point() o geom_jitter
Provincias <- st_transform(Provincias,crs = 4326) ## funcion para cambiar la proyeccion
# Provincias <- st_transform(Provincias, "+init=epsg:4326")
head(Provincias)## Simple feature collection with 6 features and 7 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: -79.76355 ymin: -3.630165 xmax: -77.5312 ymax: 1.197772
## geographic CRS: WGS 84
## # A tibble: 6 x 8
## DPA_PROVIN DPA_DESPRO DPA_VALOR DPA_ANIO REI_CODIGO REN_CODIGO PEE_CODIGO
## <chr> <chr> <int> <chr> <chr> <chr> <chr>
## 1 01 AZUAY 0 2012 05 01 593
## 2 02 BOLIVAR 0 2012 02 01 593
## 3 03 CAÑAR 0 2012 05 01 593
## 4 04 CARCHI 0 2012 04 01 593
## 5 05 COTOPAXI 0 2012 02 01 593
## 6 06 CHIMBORAZO 0 2012 02 01 593
## # ... with 1 more variable: geometry <MULTIPOLYGON [°]>## Gráfica de provincias por medio de ggplot()
Provincias %>%
dplyr::filter(DPA_DESPRO != "GALAPAGOS") %>%
ggplot() +
geom_sf() + ## graficar mapa por medio de ggplot2
labs(y = "Latidud", x = "Longitud")
##############################################
## Ilustraciones de Ecuador por parroquias ###
##############################################
### Desagregacion espacial usada por el INEC a nivel de Parroquias
Parroquias <- read_sf("nxparroquias.shp") #1040
Parroquias <- st_transform(Parroquias, "+init=epsg:4326")
Parroquias %>%
dplyr::filter(DPA_DESPRO != "GALAPAGOS") %>%
ggplot() +
geom_sf() + labs(y = "Latidud", x = "Longitud")
Carga de datos delictivos
Al momento de cargar los datos delictivos que provienen de un xlsx, se puede observar que estos tienen un formato incorrecto, tanto para la latitud como para la longitud. Esto en primera instancia puede llevarnos a pensar que no sea adecuado utilizar el componente espacial, sin embargo, lo que primero debemos hacer es trabajar estas inconsistencias por medio de “depuración espacial”.
### Carga de datos delictivos espaciales
setwd("C:/Users/Usuario/Desktop/SEE/Curso SEE 1ra edicion/Curso Espacial 1ra edicion/Material/Data Points")
DELITOS <- readxl::read_xlsx("Delitos ECU.xlsx")
DELITOS <- DELITOS %>% dplyr::mutate(Lat.new = as.numeric(latitud),
lon.new = as.numeric(longitud))
## evaluamos como está la información
summary(DELITOS)## Fecha Hora latitud
## Min. :2014-01-01 00:00:00 Min. :1899-12-31 00:00:00 Length:140479
## 1st Qu.:2014-03-31 00:00:00 1st Qu.:1899-12-31 09:15:00 Class :character
## Median :2014-06-30 00:00:00 Median :1899-12-31 13:00:00 Mode :character
## Mean :2014-07-01 01:00:15 Mean :1899-12-31 13:02:30
## 3rd Qu.:2014-10-01 00:00:00 3rd Qu.:1899-12-31 18:00:00
## Max. :2014-12-31 00:00:00 Max. :1899-12-31 23:59:00
## NA's :53
## longitud subcircuito Provincia canton
## Min. :-2.210e+16 Length:140479 Length:140479 Length:140479
## 1st Qu.:-7.850e+14 Class :character Class :character Class :character
## Median :-7.921e+13 Mode :character Mode :character Mode :character
## Mean :-1.830e+15
## 3rd Qu.:-7.791e+13
## Max. : 7.813e+15
## NA's :12
## parroquia delito Lat.new
## Length:140479 Length:140479 Min. :-4.860e+16
## Class :character Class :character 1st Qu.:-2.210e+13
## Mode :character Mode :character Median :-2.169e+06
## Mean :-1.840e+15
## 3rd Qu.: 0.000e+00
## Max. : 1.880e+16
##
## lon.new
## Min. :-2.210e+16
## 1st Qu.:-7.850e+14
## Median :-7.921e+13
## Mean :-1.830e+15
## 3rd Qu.:-7.791e+13
## Max. : 7.813e+15
## NA's :12##############################################
## Se selecciona los mil primeros delitos ###
##############################################
muestra_errores <- DELITOS[1:1000,]
##############################################
## grafico general con puntos sin tratar ###
##############################################
Provincias %>%
dplyr::filter(DPA_DESPRO != "GALAPAGOS") %>%
ggplot() +
geom_sf() +
geom_jitter(data = muestra_errores, aes(x = lon.new,y = Lat.new)) +
labs(y = "Latidud", x = "Longitud", title = "Delitos representados en forma de puntos")
##########################################################
## Delitos que caen dentro de los limites del Ecuador ###
##########################################################
muestra_correctos <-
DELITOS %>%
dplyr::filter(
Lat.new > -5 & Lat.new < 2 &
lon.new > -81 & lon.new < -75
)
# total correctos
dim(muestra_correctos) ## existen muy pocos datos correctos sin tratar## [1] 6 11## [1] 140479 11######################################################################
## Proyeccion de delitos que podrian tener una ubicacion correcta ###
######################################################################
Provincias %>%
dplyr::filter(DPA_DESPRO != "GALAPAGOS") %>%
ggplot() +
geom_sf()+
geom_jitter(data = muestra_correctos, aes(x = lon.new,y = Lat.new)) +
labs(y = "Latidud", x = "Longitud")
######################################################################
## Procesar las coordenadas de longitud ###
######################################################################
## observamos los valores de coordenadas demasiado altos
summary(DELITOS$lon.new) ### dispersion de datos## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## -2.210e+16 -7.850e+14 -7.921e+13 -1.830e+15 -7.791e+13 7.813e+15 12### El formato que se presenta es: -7919419616321.419684 o 8541635181631651
### Procesamos la longitud al tener el siguiente formato: -79.4546 --> -XX,xxxxxx
DELITOS <- DELITOS %>%
dplyr::mutate(lon.new.1 = gsub("\\.","", lon.new), ## eliminamos el punto
lon.new.1 = gsub("-","", lon.new.1), ## eliminamos el signo negativo
lon.lenh = str_length(lon.new.1), ## tomamos el tamaño total del caracter
lon.new.2 = str_sub(lon.new.1,1,2), ## tomamos solo los dos primeros valores
lon.new.3 = str_sub(lon.new.1,3,lon.lenh), ## tomamos lo restante de la coordenada
lon.new = paste0(lon.new.2,".",lon.new.3)) %>% ## pegamos ambas partes
dplyr::select(-lon.new.1,-lon.lenh,-lon.new.2, -lon.new.3) ## eliminamos columnas intermedias
## finalmente se coloca el signo negativo y se transforma a numerico
DELITOS <- DELITOS %>%
dplyr::mutate(lon.new = paste0("-", lon.new),
lon.new = as.numeric(lon.new),
Lat.new = as.numeric(Lat.new))
summary(DELITOS$lon.new) ### dispersion de datos## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## -7.990e+16 -8.000e+01 -7.900e+01 -5.079e+12 -7.900e+01 -1.000e+00 12## se retiran los puntos que aun no logran pertenecer al Ecuador
# se pierde 280 observaciones con relacion a longitud
DELITOS <- DELITOS %>% dplyr::filter(lon.new >= -81 & lon.new <= -75)######################################################################
## Procesar las coordenadas de latitud (codigo simple) ###
######################################################################
### verificacion de que los delitos pertenecen a la provincia que describe la base
### formato de punto en proyeccion para el plano
point_analisis <- data.frame(y = DELITOS$Lat.new, x = DELITOS$lon.new)
# crear puntos espaciales para proyectarlos dentro de las areas geograficas
point_analisis <- st_as_sf(
point_analisis, # puntos en formato numeric
coords=c("x","y"), # hace referencia a los puntos del dataframe point_analisis
crs=st_crs(Provincias) # menciona el CRS que se utilizara en este caso el mismo de Provincias
)
## ahora los puntos son sf o SpatialDataFrame, entoncees se le puede unir las demas caracteristicas que se desee
point_analisis$ID.temp <- 1:length(point_analisis$geometry) # unimos un identificador
DELITOS$ID.temp <- 1:length(DELITOS$Fecha) # el identificador no sirve para ubicar a los elementos en ambas tablas
### provincias sobre las que se analizaran la pertenencia topologica
provincias.puntos <- unique(DELITOS$Provincia)
### funcion de pertenencia de puntos a poligonos provincia
correctos <- c()
for(j in 1:length(provincias.puntos))
{
# extraccion de provincias homologadas para ambos documentos
base.provincia <- DELITOS %>% dplyr::filter(Provincia == provincias.puntos[j])
poligono.provincia <- Provincias %>% dplyr::filter(DPA_DESPRO == provincias.puntos[j])
# extraer puntos por medio de ids
puntos.evaluar <- point_analisis %>% dplyr::filter(ID.temp %in% c(base.provincia$ID.temp))
aux.dentro.polig <-
st_intersects(puntos.evaluar, poligono.provincia, sparse = F) # si cae en el borde cae dentro
# st_within(p, a, sparse = FALSE)[, 1] # si cae completamente dentro del poligono
# st_disjoint() # opuesto no intersects
# st_touches(p, a, sparse = FALSE)[, 1] # si los poligonos se tocan # no aplica
dentro.polig <- puntos.evaluar[aux.dentro.polig,] #indexacion para ubicar posiciones
correctos <- c(correctos, dentro.polig$ID.temp) #analizamos pertenecia simple
}
geo_correctos <- DELITOS %>% dplyr::filter(ID.temp %in% c(correctos)) # filtrado de puntos correctos
paste0("Calidad de data: ", round(dim(geo_correctos)[1]/dim(DELITOS)[1],2)*100,"%")## [1] "Calidad de data: 44%"## exclusion de puntos correctos
DELITOS <- DELITOS %>% dplyr::filter(!(ID.temp %in% c(correctos)))
######################################################################
## Procesar las coordenadas de latitud (codigo avanzado) ###
######################################################################
### funcion de pertenencia de puntos a poligonos provincia
### la LATITUD esta construida 2.0151818 o -4.206114959 ---> X,xxxxxx o -X,xxxxxx
### la LATITUD esta mal construida 20151.818 o -420611495515
### segunda limpieza ubicar posibles ubicaciones afectadas por distintos posibles casos
geo.corregir.1 <- DELITOS %>% dplyr::mutate(limpia.lat =
gsub("\\.","",as.character(Lat.new)),
limpia.lat = gsub("-","",limpia.lat),
limpia.lat = as.numeric(limpia.lat),
ID.temp = 1:length(DELITOS$Hora) )
geo.corregir <- DELITOS %>% dplyr::mutate(ID.temp = 1:length(DELITOS$Hora) )
correctos3 <- c()
for(i in 1:10)
{
if(i== 1)
{
geo.corregir.1 <- geo.corregir.1 %>% dplyr::mutate(lat.1 = str_sub(limpia.lat,1,1),
lat.leng = str_length(limpia.lat),
lat.2 = str_sub(limpia.lat,2,lat.leng),
lat.final = paste0(lat.1, ".", lat.2))
point_analisis <- data.frame(y = as.numeric(geo.corregir.1$lat.final),x=geo.corregir.1$lon.new)
point_analisis <- st_as_sf(point_analisis,coords=c("x","y"),crs=st_crs(Provincias))
point_analisis$ID.temp <- geo.corregir.1$ID.temp
## base puntos
provincias.puntos <- unique(geo.corregir.1$Provincia)
for(j in 1:length(provincias.puntos))
{
# cat(j)
base.provincia <- geo.corregir.1 %>% dplyr::filter(Provincia == provincias.puntos[j])
poligono.provincia <- Provincias %>% dplyr::filter(DPA_DESPRO == provincias.puntos[j])
puntos.evaluar <- point_analisis %>% dplyr::filter(ID.temp %in% c(base.provincia$ID.temp))
aux.dentro.polig <- st_intersects(puntos.evaluar, poligono.provincia, sparse = F)
dentro.polig <- puntos.evaluar[aux.dentro.polig,]
correctos3 <- c(correctos3, dentro.polig$ID.temp)
}
geo.corregir.1a <- geo.corregir.1 %>% dplyr::filter(ID.temp %in% c(correctos3))
geo.corregir.1 <- geo.corregir.1 %>% dplyr::filter(!(ID.temp %in% c(correctos3)))
}
# cat(i)
if(i== 2)
{
geo.corregir.1 <- geo.corregir.1 %>% dplyr::mutate(lat.1 = str_sub(limpia.lat,1,1),
lat.leng = str_length(limpia.lat),
lat.2 = str_sub(limpia.lat,2,lat.leng),
lat.final = paste0(0,".",lat.1, lat.2))
point_analisis <- data.frame(y = as.numeric(geo.corregir.1$lat.final),x=geo.corregir.1$lon.new)
point_analisis <- st_as_sf(point_analisis,coords=c("x","y"),crs=st_crs(Provincias))
point_analisis$ID.temp <- geo.corregir.1$ID.temp
## base puntos
provincias.puntos <- unique(geo.corregir.1$Provincia)
for(j in 1:length(provincias.puntos))
{
# cat(j)
base.provincia <- geo.corregir.1 %>% dplyr::filter(Provincia == provincias.puntos[j])
poligono.provincia <- Provincias %>% dplyr::filter(DPA_DESPRO == provincias.puntos[j])
puntos.evaluar <- point_analisis %>% dplyr::filter(ID.temp %in% c(base.provincia$ID.temp))
aux.dentro.polig <- st_intersects(puntos.evaluar, poligono.provincia, sparse = F)
dentro.polig <- puntos.evaluar[aux.dentro.polig,]
correctos3 <- c(correctos3, dentro.polig$ID.temp)
}
geo.corregir.1b <- geo.corregir.1 %>% dplyr::filter(ID.temp %in% c(correctos3))
geo.corregir.1 <- geo.corregir.1 %>% dplyr::filter(!(ID.temp %in% c(correctos3)))
}
# cat(i)
if(i== 6)
{
geo.corregir.1 <- geo.corregir.1 %>% dplyr::mutate(lat.1 = str_sub(limpia.lat,1,1),
lat.leng = str_length(limpia.lat),
lat.2 = str_sub(limpia.lat,2,lat.leng),
lat.final = paste0("-",lat.1, ".", lat.2))
point_analisis <- data.frame(y = as.numeric(geo.corregir.1$lat.final),x=geo.corregir.1$lon.new)
point_analisis <- st_as_sf(point_analisis,coords=c("x","y"),crs=st_crs(Provincias))
point_analisis$ID.temp <- geo.corregir.1$ID.temp
## base puntos
provincias.puntos <- unique(geo.corregir.1$Provincia)
for(j in 1:length(provincias.puntos))
{
#cat(j)
base.provincia <- geo.corregir.1 %>% dplyr::filter(Provincia == provincias.puntos[j])
poligono.provincia <- Provincias %>% dplyr::filter(DPA_DESPRO == provincias.puntos[j])
puntos.evaluar <- point_analisis %>% dplyr::filter(ID.temp %in% c(base.provincia$ID.temp))
aux.dentro.polig <- st_intersects(puntos.evaluar, poligono.provincia, sparse = F)
dentro.polig <- puntos.evaluar[aux.dentro.polig,]
correctos3 <- c(correctos3, dentro.polig$ID.temp)
}
geo.corregir.1f <- geo.corregir.1 %>% dplyr::filter(ID.temp %in% c(correctos3))
geo.corregir.1 <- geo.corregir.1 %>% dplyr::filter(!(ID.temp %in% c(correctos3)))
}
#cat(i)
if(i== 7)
{
geo.corregir.1 <- geo.corregir.1 %>% dplyr::mutate(lat.1 = str_sub(limpia.lat,1,1),
lat.leng = str_length(limpia.lat),
lat.2 = str_sub(limpia.lat,2,lat.leng),
lat.final = paste0("-0.",lat.1, lat.2))
point_analisis <- data.frame(y = as.numeric(geo.corregir.1$lat.final),x=geo.corregir.1$lon.new)
point_analisis <- st_as_sf(point_analisis,coords=c("x","y"),crs=st_crs(Provincias))
point_analisis$ID.temp <- geo.corregir.1$ID.temp
## base puntos
provincias.puntos <- unique(geo.corregir.1$Provincia)
for(j in 1:length(provincias.puntos))
{
# cat(j)
base.provincia <- geo.corregir.1 %>% dplyr::filter(Provincia == provincias.puntos[j])
poligono.provincia <- Provincias %>% dplyr::filter(DPA_DESPRO == provincias.puntos[j])
puntos.evaluar <- point_analisis %>% dplyr::filter(ID.temp %in% c(base.provincia$ID.temp))
aux.dentro.polig <- st_intersects(puntos.evaluar, poligono.provincia, sparse = F)
dentro.polig <- puntos.evaluar[aux.dentro.polig,]
correctos3 <- c(correctos3, dentro.polig$ID.temp)
}
geo.corregir.1g <- geo.corregir.1 %>% dplyr::filter(ID.temp %in% c(correctos3))
geo.corregir.1 <- geo.corregir.1 %>% dplyr::filter(!(ID.temp %in% c(correctos3)))
}
}
### Confirmacion de que no existe duplicados en las revisiones anteriores
# if(dim(rbind(geo.corregir.1a %>% filter(ID.temp %in% c(geo.corregir.1b$ID.temp)),
# geo.corregir.1a %>% filter(ID.temp %in% c(geo.corregir.1f$ID.temp)),
# geo.corregir.1a %>% filter(ID.temp %in% c(geo.corregir.1g$ID.temp)),
#
# geo.corregir.1b %>% filter(ID.temp %in% c(geo.corregir.1f$ID.temp)),
# geo.corregir.1b %>% filter(ID.temp %in% c(geo.corregir.1g$ID.temp)),
#
# geo.corregir.1f %>% filter(ID.temp %in% c(geo.corregir.1g$ID.temp))))[1]==0)
# {print("NOTA: no se debe revisar problemas")}else{print("PELIGRO : existen cruces de valores// revisar")}
### Union de resultados ultima correccion
listo.corrcion.casos.posibles <- rbind(geo.corregir.1a,
geo.corregir.1b,
geo.corregir.1f,
geo.corregir.1g)
## corregimos la coordenada
listo.corrcion.casos.posibles$Lat.new <- listo.corrcion.casos.posibles$lat.final
## eliminamos variables intermedias para correccion
listo.corrcion.casos.posibles <-
listo.corrcion.casos.posibles %>%
dplyr::select(-limpia.lat,
-lat.1,
-lat.leng,
-lat.2,
-lat.final)
## eliminamos variables intermedias para correccion
geo.corregida <-
geo.corregir.1 %>%
dplyr::filter(!(ID.temp %in% c(listo.corrcion.casos.posibles$ID.temp))) %>%
dplyr::select(-limpia.lat,
-lat.1,
-lat.leng,
-lat.2,
-lat.final)
## Revisión final de resultados
geo.corregida <- rbind(geo.corregida, listo.corrcion.casos.posibles)
paste0("Porcentaje correctos: ",
round(
(dim(geo.corregir.1a)[1]
+dim(geo.corregir.1b)[1]
+dim(geo.corregir.1f)[1]
+dim(geo.corregir.1g)[1]
+ dim(geo_correctos)[1] ) / (dim(geo.corregida)[1]+dim(geo_correctos)[1]),2)*100,"%" )## [1] "Porcentaje correctos: 100%"######################################################################
## Nueva base correcta ###
######################################################################
Delitos_Nueva <- rbind(geo_correctos,listo.corrcion.casos.posibles)
## Depuracion de datos con NA
Delitos_Nueva <-
Delitos_Nueva %>%
dplyr::mutate(Lat.new = as.numeric(Lat.new),
lon.new = as.numeric(lon.new)) %>%
dplyr::filter(!is.na(Lat.new) &
!is.na(lon.new))######################################################################
## Grafico basico ###
######################################################################
Provincias %>%
dplyr::filter(DPA_DESPRO != "GALAPAGOS") %>%
ggplot() +
geom_sf()+
geom_jitter(data =
Delitos_Nueva %>%
dplyr::filter(delito %in% c("ROBO A PERSONAS",
"HURTO",
"ROBO DOMICILIOS")) # filtrar data de interes
,aes(x = lon.new,y = Lat.new,
color = delito), alpha = 0.7) 
######################################################################
## Nivel de subcircuito Pichincha ###
######################################################################
setwd("C:/Users/Usuario/Desktop/SEE/Curso SEE 1ra edicion/Curso Espacial 1ra edicion/Material/Shape File/Semplades")
subcircuitos <- read_sf("SUBCIRCUITOS.shp") #1871
subcircuitos <- st_transform(subcircuitos, "+init=epsg:4326")
subcircuitos %>% #nivel mas bajo o desagregado
dplyr::filter(PROVINCIA == "PICHINCHA") %>%
ggplot() +
geom_sf()+
geom_jitter(data =
Delitos_Nueva %>%
dplyr::filter(delito %in% c("ROBO A PERSONAS","HURTO","ROBO DOMICILIOS"),
Provincia == "PICHINCHA")
, aes(x = lon.new,y = Lat.new, color = delito), alpha = 0.5) +
theme_minimal() + ### tema del grafico
labs(title = "Desagregación de delitos por subcircuitos de Pichincha",
color="Tipo de delitos",
x="",
y="") + #### son los titulos
theme(
axis.title.y = element_text(size = rel(1.2), angle = 90),
axis.title.x = element_text(size = rel(1.3), angle = 0),
axis.title = element_text(size = rel(1.7)),
plot.title = element_text(size = rel(1.5)),
axis.text = element_text(size = rel(1.5)),
legend.text = element_text(size = rel(1.1)),
legend.title = element_text(size = rel(1.1)))
######################################################################
## Nivel de subcircuito Guayas ###
######################################################################
subcircuitos %>% #nivel mas bajo o desagregado
dplyr::filter(PROVINCIA == "GUAYAS") %>%
ggplot() +
geom_sf()+
geom_jitter(data =
Delitos_Nueva %>%
dplyr::filter(delito %in% c("ROBO A PERSONAS","HURTO","ROBO DOMICILIOS"),
Provincia == "GUAYAS")
, aes(x = lon.new,y = Lat.new, color = delito), alpha = 0.5) +
theme_minimal() + ### tema del grafico
labs(title = "Desagregación de delitos por subcircuitos de Guayas",
color="Tipo de delitos",
x="",
y="") + #### son los titulos
theme(
axis.title.y = element_text(size = rel(1.2), angle = 90),
axis.title.x = element_text(size = rel(1.3), angle = 0),
axis.title = element_text(size = rel(1.7)),
plot.title = element_text(size = rel(1.5)),
axis.text = element_text(size = rel(1.5)),
legend.text = element_text(size = rel(1.1)),
legend.title = element_text(size = rel(1.1)))
######################################################################
## Agregacion de datos ###
######################################################################
Delitos_Nueva_Prov <-
Delitos_Nueva %>%
dplyr::group_by(Provincia) %>%
dplyr::summarise(TotalDelitos = dplyr::n())
Provincias$Provincia <- Provincias$DPA_DESPRO
#### Unir bases total de casos delictivos por poligono
Provincias <-
Provincias %>%
dplyr::left_join(Delitos_Nueva_Prov,by = "Provincia")######################################################################
## Mapa compuesto ###
######################################################################
repositorio_provincias <- "C:\\Users\\Usuario\\Desktop\\SEE\\Curso SEE 1ra edicion\\Curso Espacial 1ra edicion\\Material\\Shape File\\Provincias_Cantones_Parroquias"
setwd(repositorio_provincias)
### Cargar el mapa del mundo para anexar los paises vecinos a Ecuador
mapa_mundo <- read_sf("Paises_Mundo.shp")
names(mapa_mundo)[1] <- "Pais"
### Se filtra los paises vecinos
mapa_mundo <- mapa_mundo %>% dplyr::filter(Pais %in% c("Perú","Colombia"))
### se obtienen los centriodes para colocar los nombres o valores en cada provincia
centroides <- st_coordinates(st_centroid(Provincias[,c("geometry")])) ## se obtienen los centriodes
### se construye la nueva data para colocar el numero de delitos
centroides_texto <- data.frame(DPA_DESPRO = Provincias %>% dplyr::pull(DPA_DESPRO),
centroides, TotalDelitos = Provincias$TotalDelitos)
### Grafico con varias caracteristicas
Provincias %>%
dplyr::filter(DPA_DESPRO != "GALAPAGOS") %>%
ggplot() +
### construccion de mapa y numero total de delitos para Ecuador
geom_sf(aes(fill = TotalDelitos)) +
### se coloca otra capa para los paises aledaños
geom_sf(data = mapa_mundo,
fill = "antiquewhite1") + # (capa nueva)
### se dispone de una paleta de colores con transformaciones log para abordar la dispersion de datos
scale_fill_viridis_c(trans = "log",
alpha = 0.7,
#option = 'inferno',
direction = -1) + #### formato de colores
theme_minimal() + ### tema del grafico
## se utiliza la base de centroides para colocar el numero total de delitos
geom_text(data = centroides_texto,
aes(x = X, y = Y, label = TotalDelitos),
size =5, col= "black", fontface = "bold") +
labs(title = "Mapa de los tres principales delitos por provincia",
fill="Número de delitos",
x="", y="")+ #### son los titulos
coord_sf(xlim = c(-82, -75), ylim = c(-5, 1.5), expand = FALSE) + ## limita el campo de posibles valores
## tamaño de barra de kilometros
annotation_scale(location = "bl", width_hint = 0.3) + # escala de mapa en kilometros
annotation_north_arrow(
location = "bl", which_north = "true",
pad_x = unit(0.4, "in"), pad_y = unit(0.5, "in"),
style = north_arrow_fancy_orienteering) +
### personalizacion
theme(
axis.title.y = element_text(size = rel(1.2), angle = 90),
axis.title.x = element_text(size = rel(1.3), angle = 0),
axis.title = element_text(size = rel(1.7)),
plot.title = element_text(size = rel(1.5)),
axis.text = element_text(size = rel(1.5)),
legend.text = element_text(size = rel(1.1)),
legend.title = element_text(size = rel(1.1)),
### color de fondo
panel.grid.major = element_line(color = gray(0.5), linetype = "dashed", size = 0.5), # enrejado
panel.background = element_rect(fill = "aliceblue") ## fondo de mapa
)
Mapas interactivos con leaflet
Mas detalles: < https://rstudio.github.io/leaflet/>
######################################################################
## Mapa interactivo leaflet ###
######################################################################
### llamado de bases de datos
repositorio_datos <- "C:\\Users\\Usuario\\Desktop\\SEE\\Curso SEE 1ra edicion\\Curso Espacial 1ra edicion\\Material\\Data Points"
setwd(repositorio_datos)
load("Base_Delitos.RData")
base_final_rural <- dplyr::left_join(Provincias,
base_final_rural %>%
dplyr::rename(DPA_PROVIN = cod_provincias))
# grafica con fill
pal <- colorNumeric("viridis", NULL)
# fill color
# forma mas eficiente de construir el html
names(base_final_rural)[25] <- "Delitos"
names(base_final_rural)[26] <- "Muertes"
base_final_rural %>%
dplyr::filter(PROVINCIAS != "GALAPAGOS") %>%
leaflet() %>%
addTiles() %>%
addPolygons(
stroke = FALSE,
smoothFactor = 0.3,
fillOpacity = 0.7, ## modificar
fillColor = ~pal((Delitos)),
label = ~paste0(PROVINCIAS, ": ", formatC(Delitos, big.mark = ","))) %>%
addLegend(pal = pal,
values = ~ Delitos,
opacity = 1.0)