library(leaflet.extras)
library(apcluster)
library(magrittr)
library(dplyr)   
library(leaflet)
library(rgdal)
library(rgeos)
library(geojsonio)
library(mapview)
library(contoureR)
library(geosphere)
library(cluster)
library(sp)
library(rgdal)
library(mapview)
library(RColorBrewer)
library(leaflet.extras)
library(cluster)
library(sp)
library(rgeos)
##           [,1]      [,2]
## [1,] -51.22374 -30.02819
## [2,] -51.20018 -30.01408
## [3,] -51.10942 -30.01665
## [4,] -51.23849 -30.10901
## [5,] -51.13676 -30.07403
## [6,] -51.22705 -30.10913
## [1] 12000     2
##   Length    Class     Mode 
##     2906 APResult       S4

0.1 Clusters

centroides = unique(apres@exemplars)
poly = data.frame()
centr_indice = 0
for (i in centroides){
  centr_indice = centr_indice + 1
  centr_lat=x2[i,1]
  centr_lon=x2[i,2]
  poly = rbind(poly, c(centr_lat, centr_lon, centr_indice))
}
names(poly) = c("Lat", "Lon", "Cluster")
#head(poly)
#dim(poly)
exemplars = poly
save(exemplars, file = "exemplars.rda")
predict.apcluster <- function(s, exemplars, newdata)
{
  simMat <- s(rbind(exemplars, newdata), sel=(1:nrow(newdata)) + nrow(exemplars))[1:nrow(exemplars), ]
  unname(apply(simMat, 2, which.max))
}
resultado <- list()
#filenames <- list.files(path = paste(workdirectory, "/TCC6/geo/consolidado", sep = "")) 
#filenames
#setwd(paste(workdirectory, "/TCC6/geo/consolidado", sep = "")) 
#data <- do.call("rbind", lapply(filenames, read.csv, header = TRUE, sep = ",")) 
load("x1-total.rda")
names(data)
## [1] "lat"      "long"     "ano"      "feridos"  "fatais"   "tipoAcid"
head(data)
tail(data)
acidentes = data
dim(acidentes)
## [1] 25934     6
acidentes <- cbind(acidentes$long, acidentes$lat, as.character(acidentes$ano), as.character(acidentes$tipo), as.character(acidentes$feridos), as.character(acidentes$fatais))
acidentes <- as.data.frame(acidentes)
acidentes$V1 <- as.numeric(as.character(acidentes$V1))
acidentes$V2 <- as.numeric(as.character(acidentes$V2))
acidentes <- acidentes[acidentes$V2 < 0, ]
acidentes <- subset(acidentes, !is.na(V1))
dim(acidentes)
## [1] 25934     6
head(acidentes)
#tail(acidentes)
setwd(paste(workdirectory, "/TCC6/", sep = "")) 
save(acidentes, file = "acidentesnc.rda")
acidentes$cluster = 0
#head(acidentes)
#tail(acidentes)
acidentesfim = acidentes[1:length(acidentes$cluster),]
#head(acidentesfim)
length(acidentesfim$V1)
## [1] 25934
aa= length(acidentesfim$V1)/1000
aa = aa+1
aa
## [1] 26.934
for(i in seq(from=1, to=length(acidentesfim$V1)-1000, by=1000)){
  #print(paste(i, (i+999), sep = " - "))
  inicio = i
  final = i+999
  print(paste(inicio, final, sep = " - "))
  resultado = predict.apcluster(negDistMat(r=2), x2[apres@exemplars, ],  acidentesfim[inicio:final, 1:2])
  acidentes$cluster[inicio:final] = resultado
}
## [1] "1 - 1000"
## [1] "1001 - 2000"
## [1] "2001 - 3000"
## [1] "3001 - 4000"
## [1] "4001 - 5000"
## [1] "5001 - 6000"
## [1] "6001 - 7000"
## [1] "7001 - 8000"
## [1] "8001 - 9000"
## [1] "9001 - 10000"
## [1] "10001 - 11000"
## [1] "11001 - 12000"
## [1] "12001 - 13000"
## [1] "13001 - 14000"
## [1] "14001 - 15000"
## [1] "15001 - 16000"
## [1] "16001 - 17000"
## [1] "17001 - 18000"
## [1] "18001 - 19000"
## [1] "19001 - 20000"
## [1] "20001 - 21000"
## [1] "21001 - 22000"
## [1] "22001 - 23000"
## [1] "23001 - 24000"
## [1] "24001 - 25000"
controle = length(acidentesfim$cluster)  - final
controle
## [1] 934
resultado = predict.apcluster(negDistMat(r=2), x2[apres@exemplars, ],  acidentesfim[(final + 1):length(acidentesfim$cluster), 1:2])
acidentes$cluster[(final + 1):length(acidentes$cluster)] = resultado
#length(resultado)
acidentes = acidentes[acidentes$V1 < -49, ]
acidentes = acidentes[acidentes$V2 < -29.954337664276498, ]
#head(acidentes)
#tail(acidentes)
rm(acidentesfim)

save(acidentes, file = "acidentes.rda")

Configuramos o diretório de trabalho e carregamos as bibliotecas necessárias.

#setwd("~/OneDrive/r-files/TCC6")
setwd("/Users/fagnersuteldemoura//OneDrive/r-files/TCC6/")
options(scipen = 999)

Carregamos o dataset.

load("acidentes.rda")
head(acidentes)
acidentes$V5 = as.numeric(as.character(acidentes$V5))
acidentes$V6 = as.numeric(as.character(acidentes$V6))
dim(acidentes)
## [1] 25933     7

Verificamso os clusters contidos no dataset.

clusters_encontrados = sort(unique(acidentes$cluster))
length(clusters_encontrados)
## [1] 2905

Criamos um dataframe “poly” e o povoamos com os dados que desejamos que estejam contidos no data frame de polígono espacial.

load("exemplars.rda")
poly = data.frame()
for (i in clusters_encontrados){
  temp = acidentes[acidentes$cluster == i,  ]
  ch1 = convexHullAM_Indexes(temp[,1],temp[,2], includeColinear=FALSE,zeroBased = FALSE)
  poligono = temp[ch1, 1:2 ]
  area <- geosphere::areaPolygon(x = poligono)
  medoide_lat=length(temp[temp$V5 > 0,5]) # numero de acidentes com vitimas
  medoide_lon=length(temp[temp$V6 > 0,6]) # numero de acidentes com fatais
  pol = temp
  coordinates(pol) = ~V1+V2
  centr_lat=gCentroid(pol, byid=FALSE)$x
  centr_lon=gCentroid(pol, byid=FALSE)$y
  centr_id=0#pam(temp, 1)$medoids[,6]
  exemplars_lat = exemplars[exemplars$Cluster == i, 1]
  exemplars_lon = exemplars[exemplars$Cluster == i, 2]
  for (ii in temp$V1) {
    polying = temp[ii,]
    #polying$area = area/1000000 #km2
    polying$area = area
    polying$ocorrencias = length(temp$cluster)
    polying$centroide_lat = centr_lat
    polying$centroide_lon = centr_lon
    polying$comFeridos = medoide_lat
    polying$medoide_lon = medoide_lon
    polying$exemplars_lat = exemplars_lat
    polying$exemplars_lon = exemplars_lon
    polying$centr_id = centr_id
    
  }
  poly = rbind(poly, polying)
}
head(poly)
tail(poly)
dim(poly)
## [1] 25921    16
dados = poly
names(dados)
##  [1] "V1"            "V2"            "V3"            "V4"           
##  [5] "V5"            "V6"            "cluster"       "area"         
##  [9] "ocorrencias"   "centroide_lat" "centroide_lon" "comFeridos"   
## [13] "medoide_lon"   "exemplars_lat" "exemplars_lon" "centr_id"
#dados$V5 = sum(dados$V5)
#dados$V6 = sum(dados$V6)
head(dados)

Renomeamos os atributos do data frame e criamos novas variáveis necessárias.

names(dados) = c("lat", "lon", "ano","logr",  "feridos", "Fatais","box_id", "area", "ocorrencias", "centroide_lat", "centroide_lon", "comFeridos", "comFatais", "exemplars_lat", "exemplars_lon","centr_id")
dados$group = dados$box_id
dados$cluster = dados$box_id
dados$id = (dados$ocorrencias / dados$area)
dados$box_id = dados$id
head(dados)
rm(poly)

Criamos o data frame temporário e enxuto.

dadostemp = dados
dadostemp$lat = NULL; dadostemp$lon = NULL; dadostemp$logr = NULL; dadostemp$segmento = NULL; dadostemp$box_id = NULL
dadostemp = unique(dadostemp)
dim(dadostemp)
## [1] 10992    15
head(dadostemp)
save(dadostemp, file = "clusters_centroids.rda")

Transformamos o data set original em um objeto da classe ‘SpatialPointsDataFrame’.

coordinates(dados)=c("lat","lon")
df = dados
class(df)
## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"
str(df)
## Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
##   ..@ data       :'data.frame':  25921 obs. of  17 variables:
##   .. ..$ ano          : Factor w/ 3 levels "2017","2018",..: 1 2 2 3 1 1 2 2 2 1 ...
##   .. ..$ logr         : Factor w/ 10 levels "ABALROAMENTO",..: 6 2 1 5 2 2 1 2 5 5 ...
##   .. ..$ feridos      : num [1:25921] 1 2 1 1 2 2 1 4 1 1 ...
##   .. ..$ Fatais       : num [1:25921] 0 0 0 0 0 0 0 0 0 0 ...
##   .. ..$ box_id       : num [1:25921] 0.104556 0.104556 0.104556 0.104556 0.000675 ...
##   .. ..$ area         : num [1:25921] 38.3 38.3 38.3 38.3 7411.7 ...
##   .. ..$ ocorrencias  : int [1:25921] 4 4 4 4 5 5 5 5 5 10 ...
##   .. ..$ centroide_lat: num [1:25921] -51.2 -51.2 -51.2 -51.2 -51.1 ...
##   .. ..$ centroide_lon: num [1:25921] -30 -30 -30 -30 -30 ...
##   .. ..$ comFeridos   : int [1:25921] 4 4 4 4 5 5 5 5 5 10 ...
##   .. ..$ comFatais    : int [1:25921] 0 0 0 0 0 0 0 0 0 0 ...
##   .. ..$ exemplars_lat: num [1:25921] -51.2 -51.2 -51.2 -51.2 -51.1 ...
##   .. ..$ exemplars_lon: num [1:25921] -30 -30 -30 -30 -30 ...
##   .. ..$ centr_id     : num [1:25921] 0 0 0 0 0 0 0 0 0 0 ...
##   .. ..$ group        : num [1:25921] 1 1 1 1 2 2 2 2 2 3 ...
##   .. ..$ cluster      : num [1:25921] 1 1 1 1 2 2 2 2 2 3 ...
##   .. ..$ id           : num [1:25921] 0.104556 0.104556 0.104556 0.104556 0.000675 ...
##   ..@ coords.nrs : int [1:2] 1 2
##   ..@ coords     : num [1:25921, 1:2] -51.2 -51.2 -51.2 -51.2 -51.1 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:25921] "8611" "16138" "20627" "24236" ...
##   .. .. ..$ : chr [1:2] "lat" "lon"
##   ..@ bbox       : num [1:2, 1:2] -51.3 -30.2 -51.1 -30
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:2] "lat" "lon"
##   .. .. ..$ : chr [1:2] "min" "max"
##   ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
##   .. .. ..@ projargs: chr NA

Recriamos uma função que fará o bind dos objetos ‘SpatialPointsDataFrame’ e ‘data.frame’, criando um objeto do tipo ‘SpatialPolygonsDataFrame’.

data <- data.frame(box_id=unique(df$box_id),row.names=unique(df$id))
head(data)
#dadostemp2 = dadostemp[unique(dadostemp$id), ]
dadostemp2 = dadostemp[!duplicated(dadostemp[15]),]
#head(dadostemp2, 15)
data = cbind(data, dadostemp2)
head(data)
dim(data)
## [1] 2360   16
dadosdf = data
points2polygons <- function(df,data) {
  get.grpPoly <- function(group,ID,df) {
    Polygon(coordinates(df[df$id==ID & df$group==group,]))
  }
  get.spPoly  <- function(ID,df) {
    Polygons(lapply(unique(df[df$id==ID,]$group),get.grpPoly,ID,df),ID)
  }
  spPolygons  <- SpatialPolygons(lapply(unique(df$id),get.spPoly,df))
  SpatialPolygonsDataFrame(spPolygons,match.ID=T,data=data)
}

Criamos o objeto ‘SpatialPolygonsDataFrame’.

## class       : SpatialPolygonsDataFrame 
## features    : 2360 
## extent      : -51.26576, -51.05432, -30.24464, -29.96328  (xmin, xmax, ymin, ymax)
## crs         : NA 
## variables   : 16
## names       :                 box_id,  ano, feridos, Fatais,             area, ocorrencias, centroide_lat, centroide_lon, comFeridos, comFatais, exemplars_lat, exemplars_lon, centr_id, group, cluster, ... 
## min values  : 0.00000524248874693897, 2017,       1,      0,                0,           1,  -51.26546062,  -30.23811053,          1,         0,    -51.265622,    -30.241536,        0,     1,       1, ... 
## max values  :                    Inf, 2019,       7,      2, 953745.490234852,          95,  -51.05784783,  -29.96830419,         95,         2,     -51.05782,    -29.968397,        0,     1,    2906, ...
A caption

A caption

Salvamos o objeto ‘SpatialPolygonsDataFrame’ em um banco de dados NoSql e não tabular e sim do tipo GeoJson.

library(rgdal)
rgdal::writeOGR(obj = spDF,
                dsn = "acidentespoly.json",
                layer = "acidentespoly",
                driver = "GeoJSON",
                overwrite_layer = TRUE)

Carregamos o arquivo GeoJson salvo previamente.

Acidentes <- geojsonio::geojson_read("acidentespoly.json", what = "sp")
head(Acidentes)
## class       : SpatialPolygonsDataFrame 
## features    : 6 
## extent      : -51.2648, -51.05466, -30.24464, -29.96328  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## variables   : 16
## names       :               box_id,  ano, feridos, Fatais,             area, ocorrencias, centroide_lat, centroide_lon, comFeridos, comFatais,    exemplars_lat,     exemplars_lon, centr_id, group, cluster, ... 
## min values  : 0.000674607979607541, 2017,       1,      0,                0,           1,  -51.17887028,  -30.12305895,          1,         0,       -51.178836, -30.1230589521574,        0,     1,       1, ... 
## max values  :    0.104555562195674, 2019,       2,      0, 7734.99870450795,          18,  -51.09730782,  -30.00363156,         18,         0, -51.097546100001,  -30.003539100001,        0,     1,       6, ...

Plotamos os dados com a biblioteca mapview.
Mapa de Clusteres de Acidentes.

2º

Mapa de área total.

Mapa de números de ocorrências de Acidentes.

## [1] 10992    15
##                                                                            
## 1 function (x, ..., degree = 1, coefs = NULL, raw = FALSE, simple = FALSE) 
## 2 {                                                                        
## 3     dots <- list(...)                                                    
## 4     if (nd <- length(dots)) {                                            
## 5         dots_deg <- nd == 1L && length(dots[[1L]]) == 1L                 
## 6         if (dots_deg)

0.2 Pontos Críticos

0.3 50º Percentil de pontos Críticos

0.4 75º Percentil de pontos Críticos

0.5 99º Percentil de pontos Críticos

0.6 Dimensão de Acidentes