Clusters de Acidentes por Affinity Propagation

Clusters de Acidentes por Affinity Propagation: N=22531 = (30%) | APCluster Parameter: q=0.999 | 1169 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-99930.rda")
load("data/AZURE/apres2-20000-99930.rda")
names(x2) = c("LONGITUDE", "LATITUDE" )
head(x2)

dim(x1)

## [1] 75102     2

dim(x2)

## [1] 22531     2

Treino

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

A caption

summary(apres)

##   Length    Class     Mode 
##     1169 APResult       S4

#save(apres, file = "data/apres2-20000-99930.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] 1169    3

exemplars = poly
#save(exemplars, file = "data/exemplars-20000-99930.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-99930.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] 68842.89

median(poly$area)

## [1] 57594.13

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% 
##   2152.367  36172.491  57594.131  83866.550 262411.335

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

## [1] 7754   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    : 7754 
## extent      : -51.26617, -51.07949, -30.24132, -29.96812  (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 B CONJ RES ALTO PETROPOLIS,                     AC HAMMOND,     0, Cruzamento, ABALROAMENTO,       0 R CEL CLAUDINO, 01/01/2015,     DOMINGO, 00:00,          0,   1,      1, 2166.48207807541,         6, ... 
## max values  : 683124,        VDT JOSE EDUARDO UTZIG, VDT DOS ACORIANOS-ACESSO OESTE, 90100, Logradouro,   TOMBAMENTO, VDT JOSE EDUARDO UTZIG, 31/12/2020, TERCA-FEIRA, 23:59,         23,  13,   1169, 125491.769590304,       381, ...

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    : 993 
## extent      : -51.26617, -51.07949, -30.24132, -29.96812  (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  :      1, 601013, AC DOS HIBISCOS CONJ RES JD MEDIANEIRA,     0, Cruzamento, ABALROAMENTO,                    10 R JULIEN FERREIRA, 01/01/2015,     DOMINGO, 00:05,          0,   1,        1, 2166.48207807541,         6, ... 
## max values  :   1169, 683079,             TRVS ENG REGIS BITENCOURTH, 90010, Logradouro,   TOMBAMENTO, R VISCONDE DE MACAE & R FERNANDO CORTEZ, 31/12/2020, TERCA-FEIRA, 23:50,         23,  13,     1169, 125491.769590304,       381, ...

class(spDF)

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

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

## [1] 993  18

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

## [1] 993  24

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

A caption

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

Acidentes por Cluster

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

## [1] 993  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% 
## 12691.43

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,] 483993.3 6674640
## [2,] 481304.0 6678471
## [3,] 481955.4 6671521
## [4,] 484521.5 6663398
## [5,] 488078.3 6679704
## [6,] 481095.8 6679538

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]  2  5  2  1  1  5  2  1  1  2  2  3  1  4  1  6  2  6  3  3  5  1  5  5  5
##  [26]  5  2  2  3  5  3  3  4  1  2  1  4  7  4  1  4  3  3  1  1  5  3  1  1  5
##  [51]  2  1  1  2  4  1  1  4  3  1  5  4  5  1  5  3  3  3  2  5  3  1  3  2  8
##  [76]  3  8  2  5  1  7  7  1  2  3  5  1  4  3  3  3  5  7  4  6  4  1  4  4  5
## [101]  1  2  6  1  1  3  5  1  4  2  2  1  3  2  2  1  3  2  2  1  4  6  2  1  1
## [126]  1 10  4  2  5  4  3  3  5  5  4  1  4  4  2  1  1  2  2  6  1  2  1  3  2
## [151]  4  2  4  3  1  3  9  3  1  3  1  1  1  1  1  8  2  3  8  5  4  9  2  5  1
## [176]  1  1  5  6  5  3  3  3  4  5  1  5  6  1  1  5  1  3  1  3  9  9 10  3  1
## [201]  4  3  3  5  1  1  5  3  5  2  4  4  6  1  5  1  3  7  3  6  7  4  4  2  3
## [226]  1  4  1  3  3  4  6  4  7  2  1  3  5  5  2  2  4  4  2  1  2  5  1  3  6
## [251]  1  4  4  1  8  2  1  1  1  4  2  5  2  3  5  5  1  2  4  1  1  1  6  4  3
## [276]  1  3  3  4  1  3  4  3  7  5  3  2  8  4  6 11  5  2  2  3  4  1  1  3  5
## [301]  7  1  8  2  1  1  5  3  9  1  3  5  4  5  1  1  1  1  1  3  1  2  1  4  4
## [326]  1  4  3  6  1  2  1  4  2  6  1  1  2  4  1  1  2  5  4  5  2  2  1  1  3
## [351]  1  6  3  1  5  4  1  7  1  6  3  4  2  4  3  1  5  1  1  4  1  2  1  3  7
## [376]  1  2  1  1  6  3  3  4  2  1  7  1  5  1  1  7  1  4  7  1  6  2  1  1  1
## [401]  6  2  1  1  6  7 11  1  8  6  1  4  4  4  4  2  5  5  1 11  5  1  5  4  1
## [426]  1  4  8  2  1  2  5 11  1  1  2  5  3  1  6  6  4  5  1  4  2  3  1  1  6
## [451]  1  5  4  3  2  2  5  5  1  6  1  6  1  1  7  3  7  4  3  5  1  1  6  3  2
## [476]  6  4  1  4  6  5  1  2  3  5  1  3  4  2  5  1  4  3 10  2  1  1  3  1  4
## [501]  4  1  1  5  5  3  4  4  7  4  5  5  1  2  1  3  1  7  1  6  1  8  1  1  2
## [526]  3  4  5  4  8  3  1  2  1  1  6  6  1  3  7  3  2  2  2  7  1  3  6  1  1
## [551]  1  4  5  4  3  4  5  1  1  2  1  3  5  1  1  1  4  1  1  2  5  4  3  5  4
## [576]  1  2  1  1  7  5  3  4  3  4  6  4  2  1  8  3  1  4  1  5  1  4  3  1  1
## [601]  6  6  1  5  4  4  6  1  1  1  4  1  2  2  5  5  1  6  3  1  5  2  1  3  1
## [626]  2  2  1  5  8  3  3  3  5  1  3  2  6  2  4  3  4  4  2  1  1  3  3  1  7
## [651]  2  4 10  3  2  3  5  3  2  3  5  1  5  3  5  3  3  2  1  4  2  1  2  4  2
## [676]  2  3  4  4  5  1  1  1  4  2  2  3  1  2  2  1  1  2  5  5  2  2  1  1  4
## [701]  6  2  2  5  2  3  2  3  5  4  1  5  1  4  4  1  1  1  1  6  2  3  4  4  1
## [726]  4  3  1  6  3  2  5  5  1  5  1  3  1  4  1  3  1  3  4  5  2  2  5  5  2
## [751]  7  1  3  2  5  5  4  6  6  2  3  2  2  1  2  4  1  5  4  1  4  4  1  6  1
## [776]  3 11  4  4  1  5  1  1  3  1  3  1  5  6  4  1  2  1  5  4  4  2  4  1  4
## [801]  5  1  1  1  1  2  1  1  1  2  1  1  1  5  2  4  2  1  2  5  3  3  5  1  3
## [826]  1  4  4  7  5  1  4  1  5  3  4  2  5  5  4  1  1  1  5  3  3  3  1 10  1
## [851]  1  1  1  5  5  1  1  1  1  1  4  2  1  8  4  2  5  5  4  1  5  2  3  2  4
## [876]  6  2  5  4  4  2  2  1  5  1  3  4  2  4  3 10  5  6  1  1  4  1  5  3  5
## [901]  3  1  3  2  5  3  1  1  2  1  3  1  2  2  2  1  4  3  1  5  3  7  3  5  1
## [926]  1  1  4  3  2  1  3  9  2  2  5  1  2  3  2  3  2  1  1  1  5  1  5  5  2
## [951]  4  1  6  5  2  3  1  1  2  2  7  3  1  1  2  4  3  3  5  1  6  7  5  1  2
## [976]  3  1  1  5  5  4  5  2  4  4  7  5  2  3  1  2  3  6

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

##   [1]   9  13  22  36  40  45  48  49  52  53  56  57  72  80  87 104 105 112
##  [19] 116 125 142 146 148 155 163 164 175 190 192 194 205 214 216 226 228 236
##  [37] 245 248 251 254 257 258 259 270 272 276 315 317 318 319 321 326 330 337
##  [55] 340 341 348 349 351 354 357 366 368 371 376 385 395 400 404 411 419 425
##  [73] 435 444 449 451 461 463 482 486 491 497 499 503 515 517 532 535 538 546
##  [91] 550 551 561 564 569 576 578 589 594 596 599 603 610 612 617 623 625 628
## [109] 635 646 649 669 681 682 688 698 699 716 718 719 725 734 738 740 764 767
## [127] 770 780 782 783 787 793 802 803 805 809 811 818 824 831 842 843 852 853
## [145] 856 858 859 860 863 870 885 894 895 902 912 916 925 926 927 937 943 947
## [163] 958 970 977 978 990

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

## [1] 993

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

## [1] 826

length(dnb)

## [1] 993

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

## [1] 826   2

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

##          [,1]    [,2]
## [1,] 483993.3 6674640
## [2,] 481304.0 6678471
## [3,] 481955.4 6671521
## [4,] 484521.5 6663398
## [5,] 488078.3 6679704
## [6,] 481095.8 6679538

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

## Neighbour list object:
## Number of regions: 826 
## Number of nonzero links: 2948 
## Percentage nonzero weights: 0.4320832 
## Average number of links: 3.569007

length(dnb)

## [1] 826

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,]  644  461  547  487
## [2,]  679  604  816  516
## [3,]    8  314  324  138
## [4,]  393  354   61  429
## [5,]  582  775  589   17
## [6,]  598  223  408  336

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 = 9.587, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.2894795262     -0.0012121212      0.0009193941

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 = 9.6083, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.2894795262     -0.0012121212      0.0009153254

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.28948, 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.56754, 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.53809, 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 = 8.5562, p-value < 2.2e-16
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic       Expectation          Variance 
##      0.7299927079      1.0000000000      0.0009958277

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.0138, p-value = 5.562e-16
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic       Expectation          Variance 
##       0.729992708       1.000000000       0.001135213

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.72999, 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 = 15.126, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Global G statistic        Expectation           Variance 
##       6.020140e-03       4.326069e-03       1.254379e-08

Getis e Ord Local

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

##   [1]  0.3152738988  1.0763827938 -0.3469702677 -0.7948980813 -0.5099470616
##   [6]  0.5005046775 -0.2295330513 -0.4892876344 -1.5718969382  0.9803489422
##  [11] -0.5673449887 -0.8503490776 -1.1011240415  1.0115745178 -0.2160208251
##  [16]  4.5429194922 -0.2606004301 -0.7633882857  0.0133565525 -1.2239818059
##  [21]  0.1836399529  2.8098201493  1.0748159847  1.5126179846 -0.1311873502
##  [26] -1.6750219401 -1.0974787848 -0.5551976871 -0.2366920294  0.8098940308
##  [31]  0.3660671785  0.5177242719 -0.8506447653  1.1615133444  1.4209176952
##  [36] -0.3193289966  0.1297264727 -0.3535043770  1.1597957153  0.8955999330
##  [41]  0.0865485721  0.3009329366 -0.9971365529 -0.5183616774 -1.3158804646
##  [46] -0.1187582706 -0.5199837159 -0.2857288380  3.7826777869 -0.6751227490
##  [51]  0.9223955630  1.1398259323 -0.0646118148 -1.0716929283 -0.8132904637
##  [56] -0.6220097114  0.5188809916 -0.4059180419 -0.3201860483 -0.0344594934
##  [61] -1.3983340031  2.5308811282  2.0590215430  1.6190664837  0.5459077411
##  [66]  3.8201617678  2.8867217687  0.4569605065  1.3838353882  0.6328026695
##  [71] -0.6022352772  3.0307793013 -0.3646338273  1.3062703540  0.5875404638
##  [76] -0.2499770675  0.3529222106  1.3297134439  2.9537121726 -0.3860295324
##  [81]  1.1789991136 -0.5508023959  2.1794205124 -1.1547534510  0.6916343231
##  [86] -0.5902560633 -0.4923360866  1.2176797195  0.5533514585  0.2116177564
##  [91] -1.1602672956 -0.5642473685 -0.5212165213 -0.7367245937 -0.2578179796
##  [96] -0.1890554039 -0.3619041852 -0.8842211088  1.0088826234 -0.6323738418
## [101] -0.7345338470  0.8647928617  0.8646693809  0.0877000899  2.2798459131
## [106] -0.8158860510  3.8489328883  2.7226448925 -1.6162352849  1.0854029604
## [111] -1.5809574791 -0.1856814988  1.1892380533  0.7767896562 -0.5561815898
## [116]  0.4092514182  1.6902655180 -0.5414269026 -1.6355109545 -0.4063100149
## [121] -0.8973982415 -0.1032629365 -0.7196026850  2.4290661265 -0.0713453110
## [126] -1.0383128309  2.4643898331 -1.0839651965 -0.7954693411 -0.3381787841
## [131]  0.1539357016 -0.6494873172  3.3004775041 -0.3072686070 -0.0397580836
## [136]  0.0349224749 -0.7356467694 -0.6939149934 -1.1008528295  2.2603729014
## [141] -1.1672677088 -1.5913316898  2.2123437444  0.1356460240  3.2996222668
## [146]  4.2474673483  2.3041787703  3.0277412924  0.3473276062 -1.1639372334
## [151]  2.8371845534  0.1156328979  1.3168825434 -0.2839574131  0.2806217953
## [156] -0.5415188647  0.1779263610  2.0962317162  1.3838052884  0.0072477365
## [161]  0.1076576319  2.1983961718  2.4925326146 -1.6390305286 -0.4714134767
## [166]  5.2364242014  1.7941929571  3.2672206444 -1.6602071913 -1.2224341484
## [171]  1.1311776120 -0.6252639151  0.7395098163  1.3490714988 -0.4473013549
## [176]  1.3864073239  0.1041522925 -0.0395164596  4.0330483892 -1.2174520389
## [181]  0.7215598094  1.3639795966  0.3950747322  0.9165384025  0.1211746902
## [186] -0.6134931335  1.0672035707 -0.0636948499 -0.8995372007  0.6594109447
## [191] -0.1143734740  1.0106154222 -0.8657817882 -0.6966351514 -0.3175323990
## [196] -0.3700594425 -1.7456755849  2.5190064370  1.4106828737 -0.9960195114
## [201] -0.9181781240  0.9014116039  5.0156876521  0.9216043162 -0.7801643379
## [206]  1.0793228408  1.0096910263  0.4036284324 -1.2556541581  0.9955512101
## [211] -1.3901347103 -0.1596872619  0.8296916402 -0.7938730443 -0.5066026505
## [216] -0.5773267986 -0.4090484351  1.5115381077  1.7709105017 -0.7366257837
## [221]  0.1868794190  0.9542985053 -0.8499857034 -1.0614158337 -1.0809075972
## [226] -0.5020553455 -1.1024096461 -0.8804840549  0.7039388278  1.4297809932
## [231] -1.0377558058  0.3347153552 -0.4859406021 -1.4265189233 -0.2930643228
## [236] -1.3176964309  0.1408025908  0.5346547770 -0.5648336660 -0.3317128248
## [241]  0.1864428814  2.7472956337 -0.2965042498  2.5915788614  1.7517393158
## [246] -0.1764365729  0.6337090848 -1.1391071955  0.8117779758 -0.0793771522
## [251]  1.4858686651 -1.4879616139 -0.1904774099 -1.8943243602  0.8725779684
## [256]  0.7119182841  2.7166166524 -1.1408497410 -0.8983871043  1.7300914588
## [261]  1.2485750652 -1.4487119431  1.7253121526 -0.4499441905  0.8710871535
## [266] -0.3948462660  1.5652891027 -1.3981035437 -0.8565802382 -0.4711113862
## [271] -0.3024240705 -1.3872122627 -0.1594212282  0.1093337304 -0.9529563132
## [276] -1.4617591460  0.4170022784  0.3300651528 -0.5508684059 -0.6871841839
## [281] -0.1161238376  1.3326331398 -0.6117962317 -0.1163758966 -0.6249035326
## [286] -0.3044671785 -0.8965726647  0.1131485064  2.4008730879 -1.2544074037
## [291]  0.0283571206  0.2471157082  2.6514513039 -0.4937060920  0.4764814610
## [296] -0.9987040208  3.3449157711 -0.8980351716  1.5401879650  0.8593660712
## [301] -0.9184220642 -0.5909757120 -1.5707587105 -0.9805141198  0.9787715341
## [306] -0.8971923406  1.2731409818 -1.1542062279 -1.0814597887  1.8505530225
## [311]  0.9730548621 -0.1740050171 -1.0620733535 -0.4891218161 -0.8876241585
## [316]  1.2972167477  0.6454871269  3.3994045855 -0.0455574225  1.3333256310
## [321] -0.5108429475  0.9293485378 -1.3055856556 -0.3874537446 -0.0007864011
## [326] -0.3892585584  0.0536675443  1.2584866383  0.1487431639 -0.9811948678
## [331]  0.7522834370  0.1626837048  1.3035893739 -0.2892252005 -0.8971131537
## [336] -0.3092087266  4.5413309386  4.4345618620  0.3893188937  2.1952869547
## [341]  0.7796205832 -0.4585355045 -0.3999830132 -1.1535202904  1.5386667125
## [346] -1.1967800676 -1.4707985878  0.5458486146  2.1915565203  0.1152540107
## [351] -0.0416410076 -0.0806817175 -0.2238039962 -0.3274817909 -1.3076316982
## [356]  4.3312917597  0.8339935884 -1.0807095277 -1.4449585847  0.7008915292
## [361]  3.1376139426 -0.2460190225 -0.7514957309  1.1297784376 -1.0371326873
## [366] -0.9398159821  3.7591464229  0.3985260396  0.7350428544  0.4451845686
## [371]  0.5704825050 -0.0005238405 -0.4221182041  0.3658859612  2.2492760367
## [376] -0.2080013747  1.7645279785 -0.0919373985  0.4176242748 -0.7073824393
## [381]  2.4268285081  0.6298136059 -0.9576948946  1.6076777669  1.6934092612
## [386] -0.4498693983  0.4548399645 -0.8385827747  1.3024413028  3.8754481205
## [391] -1.1080353691  0.6652182287 -0.4289360656  0.6102013562 -0.2179884832
## [396] -1.0349897953 -0.3180012403  0.8682710261  0.0117862061 -0.3878566754
## [401]  1.3889042847  2.3017723820  2.9768637785 -0.2318612079 -1.6367174563
## [406] -0.8246455909 -1.3536429517  0.7328508250 -0.1444744577 -0.9230075796
## [411]  2.4705490212 -0.2597851907  3.3396605440 -0.8805800463 -0.9797937872
## [416] -0.8613440919  2.1016262161  1.8287460567 -0.6133519176  2.0152100226
## [421] -0.1565594048 -0.5044769907  1.6994109112  0.9875142477  3.2769370948
## [426]  0.7946633930  0.4492171495 -0.0758323444 -0.3273183241 -0.2889628924
## [431] -0.3666405826  1.2338773271 -0.8148622517  0.2886143776 -0.3868978822
## [436]  0.5610750412 -0.7563877717 -0.8567462951 -0.6354286205 -0.6011009891
## [441] -0.1029779893  3.0141031160  1.7551840615  2.7778074000  0.7317415139
## [446]  0.4591471368 -1.0599946324  0.7073637908 -0.1079088454 -1.2356570468
## [451]  0.6931090637 -1.1646379759 -0.7232075414 -0.8221176301 -1.7294939563
## [456] -0.2147496777 -0.1635508712  1.4154447058 -1.4680928090 -0.5235529697
## [461] -0.1092596353  0.0376575471 -0.8855593697 -0.6811820327 -1.6934736173
## [466] -0.6328829183 -0.1434869760  0.6182325173 -1.1204095035  0.0281610444
## [471] -1.1636156114 -1.0393935511 -0.1651417315 -1.3453422099 -0.4771881092
## [476] -0.1356243428 -0.9678691513  0.0571471311 -0.7405121249  0.8986065200
## [481]  0.4452752430 -0.7538088073  3.3011289545  1.5260643611 -0.2378112165
## [486]  0.6506886834 -0.5160924480  0.6224539538  2.2449948966  0.8620488593
## [491]  5.1724689074  2.3687914073 -0.7406485622 -0.4885056336 -1.1619514522
## [496] -0.3490161524 -1.1661897158  0.1776169519 -1.0405017854  2.4848858066
## [501]  0.2792551450  1.1155863891  2.0406939958  3.3443189515  0.9429593120
## [506]  0.5507211120  0.2025742997  1.3269751830 -0.9228372575 -0.7785394762
## [511]  3.7593619614  1.8360675495  0.1062203612 -1.1892496424 -0.8165837643
## [516]  1.4693393615 -1.0691577385 -0.0718629323  2.7083451263 -0.6361874168
## [521] -0.5681829897  1.2294571958  0.0111273267 -1.0027401989  0.7054339847
## [526]  0.0722295314 -1.4265779569  1.6701799021  3.8394133892  0.4021069940
## [531]  1.8163855345 -0.0575436785 -1.0520496879  0.2416946859 -1.6295903727
## [536] -0.6101652913 -0.2117264358  2.8942485682  3.3646425038 -1.1229071976
## [541] -0.4409027938  0.4227232789  0.7877012249 -0.3909820295  5.2247588816
## [546] -0.1476577238  0.2025964387 -0.0416429857 -0.5029981150  0.7947259313
## [551] -0.5919158945 -1.7474782241 -0.5439599737  1.1134032548  0.5056911312
## [556] -0.8608941002 -0.5506679418  0.0379000509 -1.8905757745 -1.2231416091
## [561] -0.4336195147  0.5183541080 -0.2876947332 -1.6147330289  0.9915792485
## [566]  0.3953513581 -0.7449954430  2.1008211423 -0.6135349008 -1.1527954211
## [571]  0.4898001093 -0.3183509149 -1.4015288147  0.0426582004  1.2388653219
## [576] -0.8160927178 -0.1431459410 -0.3896944093 -0.0749117195  0.2742900019
## [581]  1.0007050793 -0.7648354295 -0.5190949143 -0.2830067893  1.1660404899
## [586] -0.4180399937  2.5735162460 -0.3054152195  0.7402866452 -0.9821123592
## [591] -0.7647187436  1.1474896353 -0.0427320923 -0.6331248425  0.8780362219
## [596] -0.4688703867  1.7652394678  0.1552692387 -0.2465334335  0.9090550177
## [601] -0.4304865397  1.2119321341  0.4075505547  0.0811441349 -1.3380301526
## [606] -0.3870155448  0.2845321271  1.3039607405 -0.9565432406 -1.4437003373
## [611]  1.4141007996 -0.1655744493  1.1754644164 -1.0201288420 -1.5909834396
## [616]  0.5634497208 -0.5541609714 -0.9391442867 -0.8831189854  1.0400208886
## [621]  1.7543183040  0.6935643699 -1.4431415564  1.1288847749  0.7636962586
## [626]  1.7565280854  1.1821760187  0.7522770225  0.9165532867 -0.3892289169
## [631]  0.4003023256  0.7318770884  1.3124996154  0.9234166889  3.0446087447
## [636]  0.6619098371 -1.0720785094 -0.7221182402 -0.6494875452  1.4537775128
## [641] -0.9104508183 -1.0048040046 -0.9782002497 -0.6324719991  0.2434027613
## [646]  0.1825191157  4.3970856287 -0.2238736428 -0.2246089572  3.4832137223
## [651] -1.1551688682 -1.8476240155  0.7847511809 -1.3792163890 -1.0199314042
## [656]  0.1756866899 -0.0464121647  1.8659283105  1.1517122547 -0.5914652183
## [661] -1.7888953437  0.2827901544 -1.4394485978  0.1610961973  0.3729477196
## [666] -1.0138867966 -0.3058251517  1.8578797612  1.6433883637  1.5898858278
## [671] -1.2277437030 -1.0201794936  0.2656499932 -0.8800519780 -0.5713012350
## [676] -1.2043792156  0.1916261675 -1.0085038402  0.9799038726 -0.1152470845
## [681]  1.1954913675  0.9670178372  1.5907677012 -0.4138266387  0.6364584136
## [686] -0.8967214650 -0.3066069353 -0.8280716792 -1.1017545751  1.9642332384
## [691]  3.3488114743 -0.2509653107 -0.7145775660 -1.2336659810  0.0227684462
## [696]  0.0612929716  1.3397893505  0.1437537653  1.3231038780  0.2841071958
## [701] -0.7961625240  1.0217828769 -1.6984034597 -0.3309943767  0.3983502373
## [706] -0.4709211498  3.5279868283 -0.7935785845 -0.3476241529  0.8070391729
## [711]  0.1937721564 -0.7344941756  0.7401022978 -0.7052148012  2.6344176570
## [716] -0.2860071152 -1.1533315578  1.0419935538 -0.6596079508  0.7742193314
## [721] -0.1397929510  1.0087481692 -0.5313457850  0.1127449341 -0.4506525997
## [726]  0.2695241656  1.7011843664  1.2680068651  1.1587652391  0.3453975185
## [731] -0.8529633890 -0.9251117324 -0.7752308274  0.7115035028  1.2587117832
## [736]  0.1760635852 -0.7516019640  0.5996409921  0.3989484170  2.1160986973
## [741]  0.9012574319  0.3779830683  2.4879583924 -0.9387347004  0.1666125626
## [746]  0.9407819241 -0.0988361174 -0.5880550677 -1.8761010437  0.7196592198
## [751] -0.3503847722  0.3633667716 -0.3467314898 -0.6725690487  0.2139180242
## [756] -1.2049122769 -0.7665132935 -1.3421257017 -1.4730640206  0.6770006524
## [761]  0.3468653081  0.1506834013 -0.6525567973 -0.1989702784  0.1974825813
## [766] -0.6972306209  2.4983559610 -0.2593008359  1.3657944834 -0.0850058212
## [771]  0.8643926184 -0.1641810510  0.1637083412  3.2432752495 -0.6322892235
## [776]  0.2583426875 -0.1106951115  0.1305629880  0.5754286595  0.7193856686
## [781] -1.6024906736  0.8915728450 -0.7153129879 -0.4704757395 -0.5308677645
## [786]  0.9790776705 -0.8509279595  0.0874061357 -0.6657153764  1.4044306694
## [791]  0.7238168953 -0.9870848311  1.8154722827 -1.8274297638 -0.8579950999
## [796] -1.1374850048 -0.0731183380  0.9165386409 -1.0472705581 -1.3261460984
## [801] -1.3262730389 -0.0150343715  3.1244742598 -0.8393026819  1.1752943297
## [806]  0.2287062154  0.6677943306  0.1779799234  0.3098929111 -1.4280300412
## [811] -1.1546531272 -0.3545829379  1.7243308537 -1.4466537278  2.0586145006
## [816]  0.7008997148 -0.4753942237  0.1212515808 -0.0990498583  3.2611704141
## [821] -0.0071398481 -1.1968740229 -1.0377558058  1.5125758412 -1.1918640846
## [826] -0.2989984131
## 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.18703090 -8.461689e-04 0.3487482858 -0.3152739      0.7525537
## 2 -0.01969822 -2.036736e-06 0.0003348348 -1.0763828      0.2817561
## 3  0.11971161 -2.900637e-04 0.1196162276  0.3469703      0.7286137
## 4  0.34021687 -2.221120e-04 0.1834237878  0.7948981      0.4266728
## 5  0.25964944 -3.147266e-04 0.2598823446  0.5099471      0.6100886
## 6 -0.08214210 -1.632157e-04 0.0268279633 -0.5005047      0.6167198

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

head(CORES.5[Q])

## [1] "gray" "gray" "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"))