Clusters de Acidentes por Affinity Propagation

Clusters de Acidentes por Affinity Propagation: N=25009 = (33%) | APCluster Parameter: q=0.999 | 1210 Objetos

library(leaflet.extras)

## Loading required package: leaflet

library(apcluster)

## 
## Attaching package: 'apcluster'

## The following object is masked from 'package:stats':
## 
##     heatmap

library(magrittr)
library(dplyr)   

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(leaflet)
library(rgdal)

## Loading required package: sp

## Please note that rgdal will be retired by the end of 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## 
## rgdal: version: 1.5-27, (SVN revision 1148)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.2.1, released 2020/12/29
## Path to GDAL shared files: C:/Users/fagne/Documents/R/win-library/4.1/rgdal/gdal
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 7.2.1, January 1st, 2021, [PJ_VERSION: 721]
## Path to PROJ shared files: C:/Users/fagne/Documents/R/win-library/4.1/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.4-5
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.
## Overwritten PROJ_LIB was C:/Users/fagne/Documents/R/win-library/4.1/rgdal/proj

library(rgeos)

## rgeos version: 0.5-8, (SVN revision 679)
##  GEOS runtime version: 3.9.1-CAPI-1.14.2 
##  Please note that rgeos will be retired by the end of 2023,
## plan transition to sf functions using GEOS at your earliest convenience.
##  GEOS using OverlayNG
##  Linking to sp version: 1.4-5 
##  Polygon checking: TRUE

library(geojsonio)

## Registered S3 method overwritten by 'geojsonsf':
##   method        from   
##   print.geojson geojson

## 
## Attaching package: 'geojsonio'

## The following object is masked from 'package:base':
## 
##     pretty

library(mapview)
library(contoureR)

## Loading required package: geometry

library(geosphere)

## 
## Attaching package: 'geosphere'

## The following object is masked from 'package:geojsonio':
## 
##     centroid

library(cluster)
library(sp)
library(rgdal)
library(RColorBrewer)
library(sp)
library(foreign)

Carga de Dados

dados = read.dbf("/Users/fagne/OneDrive/r-files/CIET/acidentes2020/_Base/acidentes_2014a2020_WGS84.dbf")
dados = dados[dados$ANO > 2014, ]
sort(unique(dados$ANO))

## [1] 2015 2016 2017 2018 2019 2020

anos = length(unique(dados$ANO))
anos

## [1] 6

class(dados)

## [1] "data.frame"

x2 <- cbind(dados$LONGITUDE, dados$LATITUDE)
x2 <- x2[complete.cases(x2), ]
dim(x2)

## [1] 75102     2

head(x2)

##           [,1]      [,2]
## [1,] -51.01213 -30.19741
## [2,] -51.03732 -30.21369
## [3,] -51.05466 -30.23513
## [4,] -51.05470 -30.23517
## [5,] -51.05470 -30.23517
## [6,] -51.05470 -30.23517

Preparação

x1 <- x2
#x2 <- x2[sample(nrow(x2), round(nrow(dados)*0.30, 0)), ]
load("data/AZURE/x2-20000-99933.rda")
load("data/AZURE/apres2-20000-99933.rda")
names(x2) = c("LONGITUDE", "LATITUDE" )
head(x2)

dim(x1)

## [1] 75102     2

dim(x2)

## [1] 25009     2

Treino

#apres <- apcluster(negDistMat(r=2), x2, q=0.999)
plot(apres, x2)

A caption

summary(apres)

##   Length    Class     Mode 
##     1210 APResult       S4

#save(apres, file = "data/apres2-20000-99933.rda")

Obtenção de Centróides

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)

## [1] 1210    3

exemplars = poly
#save(exemplars, file = "data/exemplars-20000-99933.rda")

Classificação Global

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

dados$cluster = 0
for(i in seq(from=1, to=length(dados$ID)-1000, by=1000)){
  inicio = i
  final = i+999
  resultado = predict.apcluster(negDistMat(r=2), x2[apres@exemplars, ],  dados[inicio:final, 2:3])
  dados$cluster[inicio:final] = resultado
}
controle = length(dados$cluster)  - final

resultado = predict.apcluster(negDistMat(r=2), x2[apres@exemplars, ],  dados[(final + 1):length(dados$cluster), 2:3])
dados$cluster[(final + 1):length(dados$cluster)] = resultado
head(dados)

tail(dados)

#save(dados, file = "data/acidentes-20000-99933.rda")

Acidentes por Cluster

pal <- colorFactor(
  palette = 'Dark2',
  domain = dados$cluster
)

leaflet(dados) %>%
  addTiles(group="Mapa") %>% 
  addCircles(group="Acidentes", ~LONGITUDE, ~LATITUDE, weight = 0.1, radius=7, color=~pal(cluster),
             stroke = TRUE, fillOpacity = 0.8, popup=~paste("Cluster Nº: ", cluster,  
            "<br>Ano: ", ANO, "<br>Tipo: ", TIPO_ACID, "<br>Local: ", LOG1,  "<br>UPS: ", UPS,   sep = " ")) %>% 
  addLegend(group="Legenda", "topright", colors= "", labels=paste("Classificados em meio a ", summary(apres)[1], "Clusters"), title="Acidentes em Porto Alegre") %>% 
  addLayersControl(overlayGroups = c("Mapa", "Acidentes", "Legenda"),
                   options = layersControlOptions(collapsed = FALSE)) %>% 
  addProviderTiles(providers$CartoDB.DarkMatter)

A caption

Enriquecimento Informacional

rm(apres)
clusters_encontrados = sort(unique(dados$cluster))
#clusters_encontrados
parq = dados
poly = data.frame()
for (i in clusters_encontrados){
  temp = parq[parq$"cluster" == i,  ]
  ch1 = convexHullAM_Indexes(temp[,2],temp[,3], includeColinear=FALSE,zeroBased = FALSE)
  #print(i)
  #print(ch1)
  poligono = temp[ch1, 2:3 ]
  area <- geosphere::areaPolygon(x = poligono)
  acidentes = nrow(temp)
  pol = temp
  coordinates(pol) = ~LONGITUDE+LATITUDE
  centr_lat=gCentroid(pol, byid=FALSE)$x
  centr_lon=gCentroid(pol, byid=FALSE)$y
  if(nrow(temp) >= anos) {
    for (ii in ch1) {
    polying = temp[ii,]
    polying$area = area
    polying$acidentes = acidentes
    polying$centroide_lat = centr_lat
    polying$centroide_lon = centr_lon
    poly = rbind(poly, polying)
    }  
  }
}
head(poly)

tail(poly)

mean(poly$area)

## [1] 69104.99

median(poly$area)

## [1] 50608.18

minimoquantil = quantile(poly$area, probs = 0.01)
maximoquantil = quantile(poly$area, probs = 0.90)
quantile(poly$area, probs = c(0.01, 0.25, 0.5,0.75,0.99))

##        1%       25%       50%       75%       99% 
##   2385.92  33832.86  50608.18  78682.14 350499.98

poly = poly[(poly$area < maximoquantil) & (poly$area > minimoquantil), ]
dim(poly)

## [1] 8041   48

class(poly)

## [1] "data.frame"

pol = poly

Combinando Informações

dados = poly[,c(1:9, 11:13,38,41,44:48)]
head(dados)

names(dados) = c("ID","lat", "lon", "log1", "Log2", "Pred", "Local", "Tipo", "Via", "Data", "Dia", "Hora", 
                 "Fx_horaria","UPS", "box_id", "Area", "Acidentes", "CentLon", "CentLat")
dados$id = (dados$box_id * 11)
dados$group = dados$id
head(dados)

dadostemp = dados[, c(15:21)]
coordinates(dados)=c("lat","lon")
df = dados
df

## class       : SpatialPointsDataFrame 
## features    : 8041 
## extent      : -51.26266, -51.07949, -30.24132, -29.96443  (xmin, xmax, ymin, ymax)
## crs         : NA 
## variables   : 19
## names       :     ID,               log1,                       Log2,  Pred,      Local,         Tipo,                          Via,       Data,         Dia,  Hora, Fx_horaria, UPS, box_id,             Area, Acidentes, ... 
## min values  : 601001, AC A AV A J RENNER, AC SEIS TREVO ASSIS BRASIL,     0, Cruzamento, ABALROAMENTO, 0 AV PRESIDENTE JOAO GOULART, 01/01/2016,     DOMINGO, 00:00,          0,   1,      2, 2771.38786965981,         6, ... 
## max values  : 683124,   VDT OTAVIO ROCHA,   VDT ABDIAS DO NASCIMENTO, 90100, Logradouro,   TOMBAMENTO,             VDT OTAVIO ROCHA, 31/12/2020, TERCA-FEIRA, 23:58,         23,  13,   1210, 122405.631346516,       357, ...

data <- data.frame(box_id=unique(df$box_id),row.names=unique(df$id))
head(data)

dadostemp2 = dados[!duplicated(dados$id),]
data = as.data.frame(cbind(data, dadostemp2@data))

Criação de Polígonos

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 SpatialPolygonsDataFrame
data$Log2 = NULL
spDF <- points2polygons(df,data)
spDF

## class       : SpatialPolygonsDataFrame 
## features    : 1035 
## extent      : -51.26266, -51.07949, -30.24132, -29.96443  (xmin, xmax, ymin, ymax)
## crs         : NA 
## variables   : 19
## names       : box_id,     ID,                                   log1,  Pred,      Local,         Tipo,                 Via,       Data,         Dia,  Hora, Fx_horaria, UPS, box_id.1,             Area, Acidentes, ... 
## min values  :      2, 601013, AC DOS HIBISCOS CONJ RES JD MEDIANEIRA,     0, Cruzamento, ABALROAMENTO,    0 R SAO BENEDITO, 01/01/2017,     DOMINGO, 00:01,          0,   1,        2, 2771.38786965981,         6, ... 
## max values  :   1210, 683079,             TRVS ENG REGIS BITENCOURTH, 90010, Logradouro,   TOMBAMENTO, TREV CARRION JUNIOR, 31/12/2020, TERCA-FEIRA, 23:55,         23,  13,     1210, 122405.631346516,       357, ...

class(spDF)

## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"

spDF@data$group = 1
spDF@data$box_id = NULL
dim(spDF@data)

## [1] 1035   18

dadostemp = unique(dadostemp)
spDF@data = merge(spDF@data, dadostemp, by = "box_id")
dim(spDF@data)

## [1] 1035   24

plot(spDF,col=spDF$box_id+1)

A caption

library(rgdal)
rgdal::writeOGR(obj = spDF,
                dsn = "data/myParq-99933.json",
                layer = "myParq",
                driver = "GeoJSON",
                overwrite_layer = TRUE)

Acidentes por Cluster

#carregamos os dados SpatialPolygonsDataFrame
parqs <- geojsonio::geojson_read("data/myParq-99933.json", what = "sp")
#Verificamos o objeto
dim(parqs)

## [1] 1035   24

library(raster)

## 
## Attaching package: 'raster'

## The following object is masked from 'package:dplyr':
## 
##     select

## The following object is masked from 'package:magrittr':
## 
##     extract

projection(parqs)

## [1] "+proj=longlat +datum=WGS84 +no_defs"

library(mapview)
mapviewPalette(name = "Viridis")
library(RColorBrewer)
mapview(parqs, zcol = "Acidentes.x", col.regions=brewer.pal(9, "YlOrRd"))

Acidentes por m2

parqs@data$densidade = (parqs@data$Acidentes.x/parqs@data$Area.x)*1000000
mapview(parqs, zcol = "densidade", col.regions=brewer.pal(9, "YlOrRd"))

Acidentes por poligono

hist(parqs@data$Acidentes.x, col = "magenta")

A caption

Acidentes por KM2

hist(parqs@data$densidade, col = "orange")

A caption

Locais mais densos

quantile(parqs$densidade, probs = 0.99)

##      99% 
## 14500.56

temp = parqs[parqs$densidade > quantile(parqs$densidade, probs = 0.99), ]
mapview(temp, zcol = "densidade")

Calculando a matriz de vizinhanças

projection(parqs)

## [1] "+proj=longlat +datum=WGS84 +no_defs"

parqstemp = parqs
require(sf)

## Loading required package: sf

## Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1

shape <- read_sf(dsn = ".", layer = "mercator_32722_2014_2019")
projection(shape)

## [1] "+proj=utm +zone=22 +south +datum=WGS84 +units=m +no_defs"

projection(parqstemp) = projection(shape)

ccods = coordinates(parqs)
temps = as.data.frame(ccods)
cord.dec = SpatialPoints(cbind(temps$V1, temps$V2), proj4string=CRS("+proj=longlat"))
cord.UTM <- spTransform(cord.dec, CRS("+init=epsg:32722"))
ccods = as.data.frame(cord.UTM)
points = cbind(ccods[,1],ccods[,2])
head(points)

##          [,1]    [,2]
## [1,] 480144.0 6678092
## [2,] 479606.2 6675988
## [3,] 477867.8 6675924
## [4,] 485161.6 6677051
## [5,] 483806.5 6679915
## [6,] 481709.1 6680126

library(spdep)

## Loading required package: spData

## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`

distNeighbors = 400
dnb = dnearneigh(points,0,distNeighbors)
class(dnb)

## [1] "nb"

subsets = as.data.frame(matrix(dnb))
class(subsets)

## [1] "data.frame"

subsets = subsets$V1
lengths(subsets)

##    [1]  8  8  5  1  5  4  1  3  3  6  6  1  1  4  1  2  4  8  6  4  1  1  5  2
##   [25]  2  3  3  3  1  7  7  1  3  4  3  1  5 10  8  4  6  7  5  3  3  2  3  2
##   [49]  1  7  3  5  4  2  1  4  2  4  4  5  6  2  1  1  4  3  3  6  6  7  6  7
##   [73]  3  1  2  2  5  3  1  4  6  2  4  5  5  4  1  1  2  4  1  2  4  3  2  8
##   [97]  5  2  4  7  4  3  1  8  4  5  5  4  3  3  1  1  4  6  6  1  1  3  2  3
##  [121]  1  2  2  5  2  1  3  5  1  1  4  5  4  4  4  2  3  2  2  2  4  1  8  3
##  [145]  6  6  1  6  4  6  6  6  6  6  5  3  1  1  5  3  4  2  7  2  5  4  1  1
##  [169]  1  1  3  6  4  2  1  5  3  3  5  1  4  3  1  2  1  7  2  5  4  3  7  8
##  [193]  4  1  1  2  3  1  2  6  6  5  6  2  3  4  5  4  1  8  4  2  3  2  6  1
##  [217]  3  6  3  5  5  5  2  1  8  2 11  1  7  5  4  2  5  6  6  3  5 10  2  3
##  [241]  6  3  2  2  3  3  4  3  1  3  3  1  2  3  5  8  5  1  5  6  4  1  2  3
##  [265]  4  1  1  4  1  1  9  5  2  2  5  4  3  1  2  6  3  1  1 10  2  5  2  6
##  [289]  3  1  3  2  1  5  1  1  1  3  7  6  4  1  4  5  1  1  4  9  2  3  6  2
##  [313]  4  1  7  3  3  7  7  5  3  5  4  1  1  4  3  4  6  1  5  7  5  1  2  1
##  [337]  3  1  2  1  1  1  1  3  2  4  3  1  3  5  4  4  3  4  6  1  3  2  2  3
##  [361]  2  1  3  7  2  5  2  2  4  2  4  2  1  6  4  3  3  7  4  4  4  2  3  3
##  [385]  4  6  1  2  4  4  1  2  6  4  3  6  2  3  2  4  2  1  1  3  1  2  4  2
##  [409]  7  3  1  1  8  1  2  1  4  6  1  4  3  4  5  5  7  4  3  5  4  6  5  1
##  [433]  4  3  4  6  7  4  1  4  1  6  5  3  6  6  6  4  2  4  1  1  3  3  7  5
##  [457]  5  1  3  5  1  5  1  5  4  1  1  2  2  5  5  2  5  3  4  1  5  1  1  5
##  [481]  1  4  3  5  5  3  7  2  1  1  2  6  6  1  2  3  5  5  4  1  1  5  2  2
##  [505]  1  1  1  5  2  1  5  1  3  3  5  3  3  9  1  2  3  1  3  5  1  2  1  4
##  [529]  1  2  3  3  4  4  4  1  3  2  5  6  1  7  7  4  1  1  9  5  1  4  3  6
##  [553]  4  1  5  1  2  4  1  4  6  5  1  4  9  5  3  1  1  2  5  3  4  5  4  1
##  [577]  1  1  2  5  2  1  2  5  6  9  4  2  5  3  4  1  1  2  2  1  3  3  2  4
##  [601]  3  4  3  3  1  3  5  4  1  2 10  6  1  6  1  1  8  1  1  1  4  3  4  3
##  [625]  5  1  1  3  5  4  2  6  2  4  1  8  3  1  1  2  2  2  3  4  2  3  4  1
##  [649]  3  1  6  1  3  5  1  4  1  1  7  3  7  6  1  1  4  3  2  2  4  1  5  5
##  [673]  2  9  1  4  1  1  4  4  4  7  1  1  6  3  7 10  7  1  1  1  3  6  3  1
##  [697]  4  4  1  7  3  3  3  1  4  1  2  1  1  3  4  9  3  3  1  6  1  3  1  4
##  [721]  7  2  5  3  3  8  3  2  4  1  2  1  3  1  4  6  1  3  5  5  3  1  3  6
##  [745]  7  1  2  5  6  2  1  2  2  6  7  1  2  3  3  6  4  6  1  4  7  1  1  2
##  [769]  1  6  5  6  6  4  3  6  4  1  5  6  1  1  1  7  5  6  1  1  7  1  6  2
##  [793]  2  7  6  3  7  6  4  5  6  2  4  1  9  1  2  2  1  1  4  5  6  2  2  2
##  [817]  3  7  2  1  2  6  1  3  3  1  1  2  2  2  3  1  4  1  2  4  1  1  5  5
##  [841]  2  3  2  3  1  6  1  9  4  4  4  1  3  1  5  3  4  8  1  1  5  6  3  3
##  [865]  4  4  3  7  6  3  5  6  3  2  2  2  1  1  4  2  1  2  2  3  4  3  4  1
##  [889]  2  2  2  2  4  3  2  3  3  4  2  2  1  2  2  1  4  1  3  8  3  1  8  1
##  [913]  6  1  7  2  1  1  3  5  2  3  1  6  1  2  4  5  5  2  2  1  2  4  4  2
##  [937]  3  1  3  5  4  1  2  4  1  1  1  7  2  4  6  1  2  7  3  2  7  1  8  3
##  [961]  1  2  2  5  3  2  1  1  1  2  4  2  2  5  2  1  7  2  4  2  1  3  4  1
##  [985]  5  6  1  7  5  2  2  4  6  2  2  8  5  3  3  2  1  3  7  5  3  1  3  3
## [1009]  2  1  3  1  3  1  3  1  2  2  1  2  2  6  2  2  5  3  1  1  3  1  6  3
## [1033]  5  1  4

parqs$n = 1
sub = which(subsets == '0')
sub

##   [1]   13   15   21   22   29   32   36   49   63   64   88  103  112  116  121
##  [16]  126  142  147  157  158  167  168  169  175  183  194  209  216  249  252
##  [31]  258  262  269  282  290  293  302  314  336  340  341  342  348  356  362
##  [46]  373  402  403  416  419  432  439  461  467  476  479  481  489  490  494
##  [61]  506  507  522  525  536  541  545  549  554  559  563  569  577  592  593
##  [76]  596  609  615  618  620  627  639  650  663  664  677  678  691  692  696
##  [91]  704  706  709  719  732  734  737  751  756  763  767  769  778  787  790
## [106]  809  810  823  827  832  845  854  888  904  906  917  918  938  942  945
## [121]  968  976  981 1001 1006 1010 1012 1019 1034

parqs$n[sub] = 0
length(parqs)

## [1] 1035

parqs = parqs[parqs$n > 0,]
length(parqs)

## [1] 906

length(dnb)

## [1] 1035

#dim(ccods)
ccods = ccods[-sub, ]
dim(ccods)

## [1] 906   2

points = cbind(ccods[,1],ccods[,2])
head(points)

##          [,1]    [,2]
## [1,] 480144.0 6678092
## [2,] 479606.2 6675988
## [3,] 477867.8 6675924
## [4,] 485161.6 6677051
## [5,] 483806.5 6679915
## [6,] 481709.1 6680126

#dnb = dnearneigh(points,0,2000)
dnb = dnearneigh(points,0,distNeighbors)
dnb

## Neighbour list object:
## Number of regions: 906 
## Number of nonzero links: 3322 
## Percentage nonzero weights: 0.4047093 
## Average number of links: 3.666667

length(dnb)

## [1] 906

Matriz de Vizinhanca

Matriz Binária

W.Bin= nb2mat(neighbours = dnb, style = "B")
#parqs <- parqs[!sub,]

Matriz Normalizada

W.Normal= nb2mat(neighbours = dnb, style = "W")
#head(W.Normal)

KNN

vizinhos_4 <- knearneigh(points, k = 4)
class(vizinhos_4)

## [1] "knn"

head(vizinhos_4$nn)

##      [,1] [,2] [,3] [,4]
## [1,]  576  568  750  828
## [2,]  668  247  201   61
## [3,]  627   73  361  451
## [4,]  223  882  298  186
## [5,]  282  383  157  404
## [6,]  128  474  879  898

vizinhanca_4 <- knn2nb(vizinhos_4)
class(vizinhanca_4)

## [1] "nb"

Preparação para Análisis Global e Local

mv_simpl = st_as_sf(parqs)
plot(mv_simpl)

A caption

class(mv_simpl)

## [1] "sf"         "data.frame"

library(dplyr)
mv_simpl =  mv_simpl %>% dplyr::select(Acidentes.y)
#mv_simpl <- st_simplify(mv_simpl, preserveTopology = FALSE,                  dTolerance = 1)
class(mv_simpl)

## [1] "sf"         "data.frame"

mapview::mapview(mv_simpl)

sf::sf_use_s2(FALSE)#trips and tiks

## Spherical geometry (s2) switched off

mv_simpl = st_as_sf(mv_simpl)
vizinhanca_neig <- poly2nb(mv_simpl)

## although coordinates are longitude/latitude, st_intersects assumes that they are planar

ShapeNEIG = parqs
ShapeNEIG$vizinhos = card(vizinhanca_neig)
ShapeNEIG <- subset(ShapeNEIG, parqs$vizinhos != 0)
#vizinhanca2neig <- poly2nb(ShapeNEIG)

Calculando o Índice de Moran Global

Os índices de autocorrelção espacial global calculados pelos testes de normalidade e permutação.

Pelo teste de Normalidade

moran.test(parqs$Acidentes.y,listw=nb2listw(dnb, style = "W"), randomisation= FALSE)

## 
##  Moran I test under normality
## 
## data:  parqs$Acidentes.y  
## weights: nb2listw(dnb, style = "W")    
## 
## Moran I statistic standard deviate = 10.593, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.3014634960     -0.0011049724      0.0008157776

Pelo teste de Permutação ou Teste de pseudo-significˆancia

 moran.test(parqs$Acidentes.y,listw=nb2listw(dnb, style = "W"), randomisation= TRUE)

## 
##  Moran I test under randomisation
## 
## data:  parqs$Acidentes.y  
## weights: nb2listw(dnb, style = "W")    
## 
## Moran I statistic standard deviate = 10.617, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.3014634960     -0.0011049724      0.0008122204

Por simulação de Monte-Carlo

moran.mc(parqs$Acidentes.y, listw=nb2listw(dnb, style = "W"), nsim=999)

## 
##  Monte-Carlo simulation of Moran I
## 
## data:  parqs$Acidentes.y 
## weights: nb2listw(dnb, style = "W")  
## number of simulations + 1: 1000 
## 
## statistic = 0.30146, observed rank = 1000, p-value = 0.001
## alternative hypothesis: greater

Pelo teste de Permutação

Diferente dos demais testes globais o teste para o EBI é exclusivo para taxas e tem-se apenas a opção de teste da permutação

EBImoran.mc(parqs$Acidentes.y,parqs$Area.y,
            nb2listw(dnb, style="B", zero.policy=TRUE), nsim=999, zero.policy=TRUE)

## The default for subtract_mean_in_numerator set TRUE from February 2016

## 
##  Monte-Carlo simulation of Empirical Bayes Index (mean subtracted)
## 
## data:  cases: parqs$Acidentes.y, risk population: parqs$Area.y
## weights: nb2listw(dnb, style = "B", zero.policy = TRUE)
## number of simulations + 1: 1000
## 
## statistic = 0.38776, observed rank = 1000, p-value = 0.001
## alternative hypothesis: greater

Por simulação de Monte-Carlo

shapeCG.p=parqs$Acidentes.y/parqs$Area.y
moran.mc(shapeCG.p, nb2listw(dnb, style="B", zero.policy=TRUE),
         nsim=999, zero.policy=TRUE)

## 
##  Monte-Carlo simulation of Moran I
## 
## data:  shapeCG.p 
## weights: nb2listw(dnb, style = "B", zero.policy = TRUE)  
## number of simulations + 1: 1000 
## 
## statistic = 0.36979, observed rank = 1000, p-value = 0.001
## alternative hypothesis: greater

Calculando a Estatística C de Geary Global

Pelo teste de Normalidade

geary.test(parqs$Acidentes.y, listw=nb2listw(dnb, style = "W"), randomisation= FALSE)

## 
##  Geary C test under normality
## 
## data:  parqs$Acidentes.y 
## weights: nb2listw(dnb, style = "W") 
## 
## Geary C statistic standard deviate = 9.3384, p-value < 2.2e-16
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic       Expectation          Variance 
##      0.7229430619      1.0000000000      0.0008802269

Pelo teste de Permutação

geary.test(parqs$Acidentes.y, listw=nb2listw(dnb, style = "W"), randomisation=TRUE)

## 
##  Geary C test under randomisation
## 
## data:  parqs$Acidentes.y 
## weights: nb2listw(dnb, style = "W") 
## 
## Geary C statistic standard deviate = 8.7296, p-value < 2.2e-16
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic       Expectation          Variance 
##       0.722943062       1.000000000       0.001007274

Por simulação de Monte-Carlo

geary.mc(parqs$Acidentes.y, listw=nb2listw(dnb, style = "W"),nsim=999)

## 
##  Monte-Carlo simulation of Geary C
## 
## data:  parqs$Acidentes.y 
## weights: nb2listw(dnb, style = "W") 
## number of simulations + 1: 1000 
## 
## statistic = 0.72294, observed rank = 1, p-value = 0.001
## alternative hypothesis: greater

Calculando Índice de Getis e Ord Global

Getis-Ord é um indicador que mede a concentração local de uma variável de atributo distribuída espacialmente

globalG.test(parqs$Acidentes.y, nb2listw(dnb, style="B"))

## 
##  Getis-Ord global G statistic
## 
## data:  parqs$Acidentes.y 
## weights: nb2listw(dnb, style = "B") 
## 
## standard deviate = 18.064, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Global G statistic        Expectation           Variance 
##       5.843736e-03       4.051565e-03       9.843596e-09

Getis e Ord Local

localG(parqs$Acidentes.y, nb2listw(dnb, style="B"), zero.policy=NULL, spChk=NULL, return_internals=FALSE)

##   [1]  1.258265248  3.445637632  0.618404321 -0.621905575  1.696036059
##   [6]  0.354917680 -1.233785662  2.149439942 -0.902415303  0.127161958
##  [11]  1.548927648 -1.081755165 -0.101950312  0.042593542  0.656535319
##  [16]  3.174283715  1.104527106 -0.712263287 -0.919598484 -0.262988828
##  [21] -1.174901185 -0.512257131  1.533505499 -1.094168652  2.023650136
##  [26]  0.284187870 -0.665854520  2.489357350  1.321487922 -0.732348723
##  [31]  3.995368912  3.689701014 -0.461461715 -0.431322453  0.087702352
##  [36] -0.584850911 -0.563221706 -0.980830565 -0.913036522  0.190990084
##  [41]  0.244165552  3.376439615 -1.342660372  0.055675363 -0.223859400
##  [46] -1.390875205 -0.339872882  2.390309035  0.056958038  0.702634408
##  [51] -0.987224579  1.000612721  1.861331847 -0.157571075  0.006678866
##  [56] -1.734574307  3.494243315  0.458420555  0.928193781  1.257555765
##  [61]  3.809031857  5.449338640 -1.506191239 -0.689352540 -1.127974036
##  [66]  1.553587843  1.377982761 -0.225347704 -0.645758995  0.345569626
##  [71]  1.234107499 -1.004898833  2.465899268 -0.939456345 -1.183903096
##  [76]  0.449585653  2.399339120 -1.298007609  1.232387911 -0.929234539
##  [81] -0.742099530 -0.681549028  3.167538226  0.627889780  2.193890345
##  [86]  0.133951463 -1.236043714 -0.833263814  4.356242184  0.102042768
##  [91]  1.533407251  0.495857886  1.049488209 -0.164968127 -1.229798324
##  [96] -1.028920103 -0.601895179 -0.036475629 -0.732952245 -0.356327273
## [101]  3.714414770  1.961064954 -0.818937685 -0.551259322 -0.127284567
## [106]  0.215342748  1.318582449  0.522060875  0.417655225 -1.374460557
## [111] -0.992629577 -0.335799273 -0.689267920 -1.429170332  1.107091273
## [116]  2.151039719 -1.289350669  0.093656471  0.701713934  0.533071109
## [121]  0.113378307 -0.666442024 -0.171929180 -1.437048935  0.754895755
## [126]  3.464086365  1.284001034  1.034325910  3.963887630  0.696120003
## [131]  0.189439536  0.323960945  0.347228162  0.786370525 -0.471297126
## [136]  1.938390241  1.335877874  1.109163521 -0.087503624  1.494321673
## [141]  0.943037988  0.642929425  2.255162558 -1.174336031  1.637637063
## [146] -0.800453040 -0.732820116  1.196420326  1.183557744 -1.358629917
## [151] -0.403064127 -1.356664080  3.934504366 -1.040718754  0.251681897
## [156] -0.274775821  0.711226499  2.366748837 -1.036406458 -0.252984882
## [161]  2.382264444 -0.064946055  2.828487155 -0.595270874  0.452189881
## [166]  2.666691376  1.801862159  3.206221421 -1.059990194 -0.204690201
## [171] -1.142546244 -1.102012254 -0.589736486  1.215239981  0.927183470
## [176] -1.162559728  3.804958982 -0.681959592 -1.608160879  1.933548977
## [181]  0.543750702  0.835073963  3.039746710 -0.694097068 -0.144146767
## [186]  0.768112136 -1.529775219 -0.363177297 -1.256444951  0.840758475
## [191] -0.739406641  1.663473740  1.307621601  0.517417170 -0.881255046
## [196] -1.428919657  2.830588728  0.352693925  3.907036277 -1.319919551
## [201]  3.575392040  0.885542888  0.342253857  0.105222661 -0.104061145
## [206]  1.732044051  1.304363805 -0.364891826  0.243995037  4.241693918
## [211] -0.558463896 -1.129488577  0.997525642  3.433590734 -1.113290004
## [216]  1.180036168 -1.457782604  1.559347743  0.671133042 -0.775997758
## [221] -0.615996268  0.025731165  1.133515498  1.471022630  1.000755613
## [226]  4.153300300  0.475488281  1.002820159  2.431189376 -1.094999026
## [231] -1.283069487 -1.319109632 -0.854799553 -0.753846775 -0.992436884
## [236]  1.567291835 -0.950225018  3.594787857 -0.915883379  1.997748209
## [241] -0.805379077  1.284781912 -1.313176401  0.163378375 -0.885762327
## [246] -0.898129787  2.746713184 -0.161883208  0.508356482  5.105416757
## [251] -0.850371170  3.642191452  0.797222957 -0.176154433  0.644006358
## [256] -0.816891069  1.234222160 -0.988273547 -0.057513317 -0.798222691
## [261] -0.992937727 -0.775529388  1.005589573  0.829859005  0.112979594
## [266]  0.353901834 -0.869861471 -0.296914602 -1.123494435 -0.611776850
## [271]  3.890173477 -1.298664616  0.523612036  1.142132752 -1.018063494
## [276] -0.301103737  1.001140640  0.352066568 -0.062492426 -0.083921760
## [281]  0.443084313  1.805901747 -0.462163386 -0.403799853  0.081837921
## [286] -0.819086288 -1.168815052 -1.214161161  0.930829656 -1.357175044
## [291] -0.354885033  0.247264961  0.299828479  2.427721572  1.014659160
## [296] -1.103397303 -1.420677563 -0.373772029 -0.710694278 -0.680053869
## [301] -1.122505948 -1.644726646 -0.481432254 -0.590868880 -1.056518666
## [306]  1.069437428 -0.733063796  0.777050267 -1.670362831 -0.477726982
## [311]  0.276509551 -0.567942407 -0.639722709 -1.127682884 -1.298133543
## [316] -0.138784921  0.180380033  0.415384286  3.634372777 -1.159453919
## [321] -0.135266905 -1.235045429 -1.298070235  3.274315418 -1.668909089
## [326] -1.215798841  0.226198385  1.098330955  0.972584394  0.027355781
## [331] -0.362812186  1.066539451 -0.094777237  0.343562356  1.130576093
## [336]  0.071551358 -0.024578868 -0.262702990  0.844967604 -1.523348320
## [341] -0.187912594 -1.681531951  0.842843488 -0.069650684 -0.797302354
## [346] -1.375277968  1.812085121 -1.748038979  0.792564747  2.432572042
## [351] -1.406058129 -0.112850578  1.979704683 -0.593623358 -0.898413007
## [356]  0.176337423 -0.732249451 -0.881823630  0.069514004 -0.929137902
## [361]  1.615311153 -0.616266081 -1.341843256 -1.296966555  2.740541287
## [366] -0.971077337  0.964304415 -0.246442911  0.109899495 -1.041743169
## [371] -0.287173943  0.875632233  1.450226081  1.509232935  1.748134053
## [376]  2.386088240 -1.031243596  2.041098262  1.255770284  0.575647519
## [381]  0.045239386 -0.558976375  0.131992312 -0.321561096  1.024417893
## [386] -0.493016153  0.279265145  0.636473993  0.051818495  1.693830097
## [391]  1.238327473 -1.280236667  2.224896045  1.623133407  1.193448517
## [396] -0.931258757 -0.342769335  1.312026883 -0.144429879 -1.343005415
## [401] -0.676957830  1.473229329 -0.878960074  2.173187959 -0.226336677
## [406] -1.124251860  0.163612137  1.183750699  3.828020222 -0.254148089
## [411] -0.275004569  1.116737324 -0.078914241 -0.681020871 -1.790967457
## [416] -0.525830700  0.318374697 -1.145423838  1.202171994 -1.509199184
## [421]  1.501987444  0.221457319 -1.189961526 -0.295176343 -0.361883091
## [426] -0.289272379 -0.869020982 -0.774382627  0.064425891  0.373910754
## [431] -1.004527527 -0.650257743  1.578263613  3.828614368  1.770713847
## [436] -0.954246060  1.344591945  0.539146521 -0.081536842 -0.646409622
## [441] -0.774735191 -1.094476740  0.917868152 -1.064444881  0.115710573
## [446] -0.422602828  0.810158035 -0.601172004  0.070669753  1.921092842
## [451]  0.831446669 -1.482792162 -0.032254807  0.869982925  0.390062958
## [456]  2.915156075 -0.145832358 -0.204581752  0.026436747 -1.204269988
## [461] -0.780882835 -0.098047345 -0.405299574  0.953180105 -1.079919494
## [466] -0.098194226 -0.502969467  0.126026393  0.212696136 -0.431903534
## [471]  0.777589639  0.091139098  0.811354818  1.012331570 -1.210758639
## [476] -0.521036008  3.594085277  1.103305152 -0.906072305  2.883336673
## [481]  1.265921148  0.560208928 -1.119869438  1.158044932 -0.029174054
## [486] -0.947696663 -0.993931166 -0.405925172 -0.486586105 -0.114682545
## [491]  1.451313136  0.465194736 -0.147688309  4.201652413  1.522684134
## [496] -1.168585376 -0.885058086 -0.848818692  0.387736246 -0.072048960
## [501] -0.494057640  0.565151728 -0.877608719 -0.341296036 -0.492081647
## [506] -1.496204193 -1.640107884 -0.420251657  1.942825881 -1.342006496
## [511]  1.675126870  2.451012975  1.792276868 -0.355186839 -0.373425170
## [516]  2.403136295 -0.049586725 -0.276306581  0.658483888  0.149955666
## [521] -0.803947556 -1.443127365 -1.360232853  2.106044715 -0.564240633
## [526]  0.384767692  0.140718572  0.416572596 -0.928567012 -0.250094797
## [531]  1.315021738  1.736381359 -1.604039657  4.702136586 -0.158770710
## [536] -0.752844214  0.826957754 -1.080449493  2.946509755 -1.233073740
## [541]  0.772700017 -0.919154500 -1.116886119 -0.113551046  1.127800877
## [546] -0.994410564 -1.320683630 -0.009262718 -0.703008385 -0.927701136
## [551]  1.890702435  0.658750669 -1.247551256 -0.427405771  0.825418562
## [556] -1.534828876 -0.732907815 -1.529470548  1.518515243 -1.021607547
## [561] -1.698746895  0.180627319 -1.264707330 -0.199915950  0.573643177
## [566] -0.188356375 -0.777178411  1.620798200 -0.448728978 -0.137297964
## [571]  2.735829770 -0.122638983  0.508463341 -0.709142999 -0.428844927
## [576]  2.553957605  1.231908031 -0.299356513  0.321849000  0.440862366
## [581] -0.740569008 -1.097461706 -0.930197484 -0.165237319 -0.710539205
## [586]  0.479276724  2.227207216 -1.437193150  3.097187351 -1.081340410
## [591]  0.333401774 -1.358629917  1.366178241 -0.681712802  0.035617532
## [596]  0.093871168 -0.253085541 -0.203092922 -0.602962945  1.312131784
## [601]  4.382751894  3.376180162 -0.949332855 -0.500531744  0.634733239
## [606]  0.341219244  0.222292071  1.147635609  1.247078290  1.637893468
## [611] -0.199415303 -0.225218449  1.370323708  1.444770514 -0.958698029
## [616] -0.752868288  1.935691418  3.906770705  3.022838577 -0.767005194
## [621]  0.051360734 -0.558515117  2.091995230 -0.253293228 -1.496626576
## [626] -0.973771570  2.122486685 -0.127388368  0.767628055  0.936962270
## [631] -0.125202308  2.694333946  0.465018153  0.579166627  0.788604403
## [636] -0.319096553  0.472628660 -0.527752646  0.606693520  1.281813483
## [641] -1.255790312 -0.520080988  1.236279757  0.603675174  0.311847109
## [646] -1.329611760  2.288414637  0.178264442 -1.056997651 -1.066958530
## [651]  0.708605211  1.200434943  0.181902528 -0.713056917  0.536786102
## [656] -0.194372170  0.713467344 -0.911059209 -0.915954509  0.038256938
## [661]  0.502125774 -1.083361911 -1.512482613 -0.942720537  2.570132871
## [666] -0.123276843 -0.405191393  3.264632064 -0.022075213  1.368289939
## [671]  1.237602731 -1.485824340 -0.162430379  2.299863454 -0.324267451
## [676] -0.674134601  1.595994288  0.267997869 -1.449165043 -1.320700555
## [681] -0.843516852  1.113642845  4.408200555 -0.384005795  2.934818207
## [686]  1.915740694 -0.745228542 -0.188776127  0.944765286  0.995554973
## [691] -0.075489058  1.169528550  1.052133230 -1.216204560  0.895287144
## [696]  1.740008839 -0.389805239  0.551917228 -0.993402851  3.056354333
## [701] -0.971978373 -0.098527853 -0.358409008  0.202155352  0.913354707
## [706]  0.769112732 -1.097514728 -0.618921535 -1.545652566  0.138594315
## [711]  2.254021543  0.749257154 -0.623768949 -1.390025081  1.878889040
## [716] -0.527853442 -0.086033793 -0.297408015  1.873668947  2.502747276
## [721] -0.142700730 -0.827166872 -0.834278722 -1.015107341  1.718650030
## [726]  1.576561031 -0.579840966 -0.819011127  0.261801926  0.639400947
## [731] -0.186079025 -0.414394223 -0.943915153 -0.413450632  0.901948265
## [736]  1.444010265  5.017390358 -0.995742307 -0.508198305  1.736571575
## [741] -0.993982335 -0.149992229  0.975152594 -0.074677966  0.174107871
## [746]  4.160158997 -0.949408349 -1.102655690 -0.203715132  2.654043485
## [751]  0.229611918 -0.261961514  3.782921567  0.440528923 -0.288223356
## [756]  2.963542752 -0.068060857 -0.013100835  1.888397191  1.193262955
## [761]  0.087482457  0.395819919  0.025788767  1.642103516  0.073221010
## [766] -1.059752664 -1.641787358  0.703227982 -1.363257570 -1.420798537
## [771] -1.051304868 -0.840009740  1.204028026 -0.390599134  1.899661246
## [776]  1.688330925  0.101950374 -0.018090291  0.749285612 -0.093245165
## [781] -0.627611513 -1.065484600 -0.387507635 -1.067055383 -1.269108056
## [786] -1.329041627 -0.865913699 -0.796584624  0.411234621 -0.237022961
## [791]  0.369261835  0.216899416  3.509171483 -0.817115230 -0.319005865
## [796] -0.075566755  0.725186705  3.080665722 -0.296441310  0.003280153
## [801]  0.521341971 -0.563014805  0.017067370  0.594885554  2.921576390
## [806] -0.123224317  1.456764890 -1.037933309 -0.636121784 -0.300279207
## [811]  1.615685587  0.699661795  1.151976129 -0.806345348  0.681263239
## [816]  0.164127381  0.352430401  1.520713086  1.242222362 -0.073798192
## [821]  0.754454205 -0.378723753  0.047639388  0.011557300 -0.823563361
## [826] -0.949295854 -1.059468350  1.884819328  0.087756407  0.276419543
## [831]  0.703423744 -0.428577334 -0.603571045  1.437334092 -1.003137877
## [836] -0.605347731 -0.434105495 -0.231369818  1.651643217 -0.602338902
## [841] -0.798222691 -1.774972791  1.256577658 -0.526053361  0.213248101
## [846] -0.836570953 -0.601690475 -0.797573159  0.025286139  1.736995597
## [851] -0.220146615  0.995200288 -1.520474929  0.071260776  0.184022200
## [856] -0.420536036 -0.376002820  0.026327997 -1.092141101 -0.635725372
## [861]  0.356143186  3.214188731  4.620576778 -0.819979506  2.416508941
## [866] -1.022352899 -1.467688667  0.503327704  1.960439818  1.613656524
## [871] -1.726217420 -0.189860319  5.145276700  1.130673631  0.591429272
## [876] -1.030926978  1.502706845 -0.981084800  1.640902501  0.338961801
## [881] -1.180867731 -0.422366219 -0.463379509  1.087361617 -0.474866016
## [886] -0.804833917 -1.121949071  0.982932827 -0.688399069 -1.851164291
## [891] -0.790401827 -1.251966574 -1.159224165  1.089284764  0.056032143
## [896]  0.364429255 -0.149607219  0.578825728  0.159136269  0.681533088
## [901]  0.529234098 -1.145327731  0.409651696 -0.855361691  0.103315035
## [906] -1.650555205
## attr(,"gstari")
## [1] FALSE
## attr(,"call")
## localG(x = parqs$Acidentes.y, listw = nb2listw(dnb, style = "B"), 
##     zero.policy = NULL, spChk = NULL, return_internals = FALSE)
## attr(,"class")
## [1] "localG"

Moran Local

Todas as análises feitas até o momento foram de escala global. No entanto, é necessário que seja feita também uma análise local do estudo. Essa análise pode ser feita pelo índice local de autocorrelaçãoo espacial (LISA). Para isso é preciso calcular o índice de Moran local.

ShapePB.mloc <- localmoran(parqs$Acidentes.y, listw=nb2listw(dnb, style="W")) 
head(ShapePB.mloc)

##            Ii          E.Ii     Var.Ii       Z.Ii Pr(z != E(Ii))
## 1  0.83698866 -0.0039911804 0.44671113  1.2582652   0.2082958447
## 2  4.27062123 -0.0139545979 1.54623843  3.4456376   0.0005697138
## 3  0.51969477 -0.0039911804 0.71712823  0.6184043   0.5363088504
## 4 -0.34368470 -0.0003365414 0.30480393 -0.6219056   0.5340039594
## 5  0.48322437 -0.0004510259 0.08132758  1.6960361   0.0898790514
## 6  0.09364609 -0.0003105345 0.07008088  0.3549177   0.7226512629

Mapa das probabilidades (Signific?ncias do I de Moral Local)

Por meio dos valor-p do éndice de Moran local é possível construir um mapa de probabilidades.

library(classInt)
INT4 <- classIntervals(ShapePB.mloc[,5], style="fixed", 
                       fixedBreaks=c(0,0.01, 0.05, 0.10))
CORES.4 <- c(rev(brewer.pal(3, "Reds")), brewer.pal(3, "Blues"))
COL4 <- findColours(INT4, CORES.4)
parqs$COL = COL4  
parqs$p_valor = ifelse(parqs$COL == "#DE2D26", "[0,0.01)", ifelse(parqs$COL == "#EEE5E4", "[0.01,0.05)", "[0.05,0.1]"))
plot(parqs, col=COL4)
title("P-valores do I de Moran Local por Distäncia de Centróides")
TB4 <- attr(COL4, "table")
legtext <- paste(names(TB4))
legend("bottomright", fill=attr(COL4, "palette"), legend=legtext, 
       bty="n", cex=0.7, y.inter=0.7)

A caption

mapview(parqs, zcol = "p_valor", col.regions=c("red", "orange", "green"))

temp = parqs[parqs$p_valor != "[0.05,0.1]", ]
mapview(temp, zcol = "p_valor", col.regions=c("red", "orange"))

Montando matrix W de vizinhança

ShapeCG.nb1.mat <- nb2mat(dnb)

Incidência de acidentes padronizada

Acidentes_SD <- scale(parqs$Acidentes.y)

Média das incidências de acedentes padronizada

Acidentes_W <- ShapeCG.nb1.mat %*% Acidentes_SD

Diagrama de espalhamento de Moran

plot(Acidentes_SD, Acidentes_W,xlab="Z",ylab="WZ")
abline(v=0, h=0)
title("Diagrama de Espalhamento de Moran por Distancia de Centróides")

A caption

Q <- vector(mode = "numeric", length = nrow(ShapePB.mloc))
Q[(Acidentes_SD>0  & Acidentes_W > 0)] <- 1            
Q[(Acidentes_SD<0  & Acidentes_W < 0)] <- 2
Q[(Acidentes_SD>=0 & Acidentes_W < 0)] <- 3
Q[(Acidentes_SD<0  & Acidentes_W >= 0)]<- 4
signif=0.05
parqs$Q = Q

Mapa LISA

Q[ShapePB.mloc[,5]>signif]<-5
CORES.5 <- c("blue", "green" , "red", "yellow", "gray", rgb(0.95,0.95,0.95))
#CORES.5 <- c(1:5, rgb(0.95,0.95,0.95))
parqs$cores5Q = CORES.5[Q]
plot(parqs, col=CORES.5[Q])
title("Mapa LISA por Distancia Centroides")
legend("bottomright", c("Q1(+/+)", "Q2(-/-)", "Q3(+/-)", "Q4(-/+)","NS"), 
       fill=CORES.5)

A caption

CORES.5[Q][1:5]

## [1] "gray" "blue" "gray" "gray" "gray"

head(CORES.5[Q])

## [1] "gray" "blue" "gray" "gray" "gray" "gray"

#save(parqs, file = "parqsFinal-975.Rds")

parqs$cores5 =  ifelse(parqs$cores5Q == "blue", "A-A", ifelse(parqs$cores5Q == "green", "B-B", 
            ifelse(parqs$cores5Q == "red", "A-B", ifelse(parqs$cores5Q == "yellow", "B-A", "NA"))))
mapview(parqs, zcol = "cores5", col.regions=c("red", "orange", "green", "yellow", "grey"))

temp = parqs[parqs$cores5 == "A-A", ]
mapview(temp, zcol = "cores5", col.regions=c("red", "orange", "green", "yellow", "grey"))

mapview(temp, zcol = "Acidentes.x", col.regions=brewer.pal(9, "YlOrRd"))