Avaliação

Parâmetros de Análise

# dimensão dos hexagonos em m
cel_size = 75
# dimensão do buffer em m
buffer = 100
# distancia de vizinhanã entre centróides em m
DistanciaCentroides =cel_size*5
# unicio do pico avaliado
inicioPico = 6 # >=
# fim do pico avaliado
FimPico = 8    # <=
# velocidade mínima de fila admitida em km/h
MinDelay = 0
# velocidade mínima de fila admitida em km/h
MaxDelay = 600
set.seed(77)
#set.seed(3)
#set.seed(1)
#set.seed(1047)
# usuário do diretório
user = "fagne"
#user = "fsmoura"
require(sp)
## Carregando pacotes exigidos: sp
library(mapview)
library(rgeos)
## rgeos version: 0.6-4, (SVN revision 699)
##  GEOS runtime version: 3.11.2-CAPI-1.17.2 
##  Please note that rgeos will be retired during October 2023,
## plan transition to sf or terra functions using GEOS at your earliest convenience.
## See https://r-spatial.org/r/2023/05/15/evolution4.html for details.
##  GEOS using OverlayNG
##  Linking to sp version: 2.1-3 
##  Polygon checking: TRUE
library(sf)
## Linking to GEOS 3.11.2, GDAL 3.8.2, PROJ 9.3.1; sf_use_s2() is TRUE
library(raster)
library(DT)
library(tmap) #https://rpubs.com/leydetd/maps
# Função para titulo nas tabelas DT
js <- c(
  "function(settings){",
  "  var datatable = settings.oInstance.api();",
  "  var table = datatable.table().node();",
  "  var caption = 'ANOTHER CAPTION'",
  "  $(table).append('<caption style=\"caption-side: bottom\">' + caption + '</caption>');",
  "}"
)
####################
# ZONAS  DE TRAFEGO
# Carrega ZTs
#zts <- read_sf(paste("C:/Users/", user,"/OneDrive/EngenhariaDeTrafego/EDOM2023/P15/Rev03/20241218_P15 Relatorio Final REV02/P15_Mapas/1_Zoneamento/P15_zoneamento_n1.shp", sep=""))

Vias primárias e secundárias segundo o OSM

 # VIAS DE PRIMEIRO E SEGUNDO N?VEL
# Carrega as vias de interesse
bordeada <- read_sf(paste("C:/Users/", user,"/OneDrive/PPGEP-LASTRAN/Doutorado/R/OSM_PoA_primariasSecundarias/Bordeada", buffer,"m.shp", sep=""))
projection(bordeada)
## [1] "+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +vunits=m +no_defs"
bordeada = as(bordeada, "Spatial")
class(bordeada)
## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"
mapview(bordeada, color="red", col.regions="magenta",
        layer.name = c("Eixos Avaliados"))  #https://gis.stackexchange.com/questions/343815/custom-legend-mapview

Hexágonos

Criação de poligonos em formato de hexágonops para varredura da rede utilizando a função HexPoints2SpatialPolygons.

class(bordeada)
## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"
projection(bordeada)
## [1] "+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +vunits=m +no_defs"
bordeada
## class       : SpatialPolygonsDataFrame 
## features    : 1 
## extent      : -5706980, -5683045, -3513168, -3476968  (xmin, xmax, ymin, ymax)
## crs         : +proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +vunits=m +no_defs 
## variables   : 10
## names       :   osm_id,             name,   highway, waterway, aerialway, barrier, man_made, railway, z_order,                                                                                                                                                                        other_tags 
## value       : 14691002, Avenida Farrapos, secondary,       NA,        NA,      NA,       NA,      NA,       6, "bicycle"=>"yes","foot"=>"yes","lanes"=>"4","maxspeed"=>"60","oneway"=>"yes","oneway:bicycle"=>"no","sidewalk"=>"right","sidewalk:surface"=>"paving_stones","surface"=>"concrete"
HexPts <-spsample(bordeada, type="hexagonal", cellsize=cel_size)
HexPols <- HexPoints2SpatialPolygons(HexPts)
#mapview(zts, zcol="ZT1", layer.name = c("ZTs")) +  
 # mapview(HexPols, color="red", col.regions="magenta", layer.name = c("ZTs2"))
#clear_viewer_pane()
class(HexPols)
## [1] "SpatialPolygons"
## attr(,"package")
## [1] "sp"
projection(HexPols)
## [1] "+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +vunits=m +no_defs"

FILAS

Carrega os dados de filas

# FILAS
# CArrega os dados de filas
#load(paste("C:/Users/", user,"/OneDrive/PPGEP-LASTRAN/Doutorado/R/filasfull.Rda", sep=""))
load(paste("C:/Users/", user,"/OneDrive/PPGEP-LASTRAN/Doutorado/R/filas.Rda", sep=""))
max(filas$date)
## [1] "2023-12-29"
min(filas$date)
## [1] "2023-03-01"
filas <- filas[which(filas$city == "Porto Alegre"), ]
#filas = filas[filas$delay == "-1", ]
filas =  filas[(filas$delay >=MinDelay) & (filas$delay <= MaxDelay), ]
                                       #   filas = filas[filas$date < "2023-04-01",]
filas$identificador = 1:nrow(filas)
head(filas[ , c("lat", "lon")])
##              lat       lon
## 510000 -30.04026 -51.22590
## 201070 -30.03572 -51.19175
## 381300 -29.98100 -51.12687
## 481060 -30.00623 -51.12775
## 491060 -30.03864 -51.10956
## 511030 -30.02053 -51.21466

Caracterização das filas

library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:raster':
## 
##     intersect, select, union
## The following objects are masked from 'package:rgeos':
## 
##     intersect, setdiff, symdiff, union
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
numero = 1
datatable(head(filas, 10), caption = htmltools::tags$caption(
  style = 'caption-side: top; text-align: left;',
  paste('Tabela', numero), htmltools::em(iconv(paste('Filas no período entre', inicioPico, "h e ",FimPico, "h.", sep = ""), from = "", to = "UTF-8"))
))
numero = numero +1
glimpse(filas)
## Rows: 266,432
## Columns: 23
## $ country       <chr> "BR", "BR", "BR", "BR", "BR", "BR", "BR", "BR", "BR", "B…
## $ city          <chr> "Porto Alegre", "Porto Alegre", "Porto Alegre", "Porto A…
## $ level         <int> 2, 3, 3, 3, 4, 3, 3, 1, 2, 3, 2, 3, 3, 3, 2, 2, 3, 3, 3,…
## $ line          <list> [<data.frame[2 x 2]>], [<data.frame[4 x 2]>], [<data.fr…
## $ speedKMH      <dbl> 6.29, 6.14, 6.39, 15.63, 5.19, 7.71, 5.69, 28.74, 22.23,…
## $ length        <int> 226, 233, 326, 595, 275, 365, 224, 1574, 933, 485, 952, …
## $ turnType      <chr> "NONE", "NONE", "NONE", "NONE", "NONE", "NONE", "NONE", …
## $ type          <chr> "NONE", "NONE", "NONE", "NONE", "NONE", "NONE", "NONE", …
## $ uuid          <int> 630554630, 630556216, 630552304, 632032998, 631887650, 6…
## $ endNode       <chr> "", "", "Av. Alcides Maia", "Av. Gen. Raphael Zippin", "…
## $ speed         <dbl> 1.747222, 1.705556, 1.775000, 4.341667, 1.441667, 2.1416…
## $ segments      <list> [<data.frame[1 x 0]>], [<data.frame[2 x 0]>], [<data.fr…
## $ roadType      <int> 1, 1, 1, 2, 1, 2, 1, 2, 2, 2, 2, 2, 1, 2, 1, 1, 1, 2, 1,…
## $ delay         <int> 72, 103, 126, 91, 129, 132, 104, 76, 71, 166, 70, 116, 7…
## $ street        <chr> "R. da República", "R. Barão de Ubá", "Av. Souza Melo", …
## $ pubMillis     <dttm> 2023-03-01 00:00:01, 2023-03-01 00:00:01, 2023-03-01 00…
## $ hora          <dbl> 1.157407e-05, 1.157407e-05, 0.000000e+00, 4.098380e-02, …
## $ dia           <chr> "Wednesday", "Wednesday", "Wednesday", "Wednesday", "Wed…
## $ date          <date> 2023-03-01, 2023-03-01, 2023-03-01, 2023-03-01, 2023-03…
## $ lon           <dbl> -51.22590, -51.19175, -51.12687, -51.12775, -51.10956, -…
## $ lat           <dbl> -30.04026, -30.03572, -29.98100, -30.00623, -30.03864, -…
## $ delay60       <dbl> 120, 120, 180, 120, 180, 180, 120, 120, 120, 180, 120, 1…
## $ identificador <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 1…
coordinates(filas) = ~lon+lat
projection(filas) = "+proj=longlat +datum=WGS84 +no_defs"
class(filas)
## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"
head(filas@data)
##        country         city level
## 510000      BR Porto Alegre     2
## 201070      BR Porto Alegre     3
## 381300      BR Porto Alegre     3
## 481060      BR Porto Alegre     3
## 491060      BR Porto Alegre     4
## 511030      BR Porto Alegre     3
##                                                                                                                                                                                                                                                    line
## 510000                                                                                                                                                                                                       -51.22408, -51.22590, -30.03899, -30.04026
## 201070                                                                                                                                                           -51.18938, -51.18986, -51.19146, -51.19175, -30.03558, -30.03566, -30.03582, -30.03572
## 381300                                                                                                                                     -51.12350, -51.12503, -51.12595, -51.12670, -51.12687, -29.98084, -29.98090, -29.98095, -29.98099, -29.98100
## 481060                                                                                                                                                           -51.12857, -51.12788, -51.12779, -51.12775, -30.01155, -30.00713, -30.00652, -30.00623
## 491060 -51.10773, -51.10816, -51.10849, -51.10863, -51.10884, -51.10900, -51.10903, -51.10907, -51.10910, -51.10912, -51.10956, -30.04024, -30.04005, -30.03994, -30.03987, -30.03969, -30.03948, -30.03936, -30.03919, -30.03891, -30.03881, -30.03864
## 511030                                                                                                               -51.21633, -51.21620, -51.21579, -51.21561, -51.21556, -51.21466, -30.02349, -30.02330, -30.02269, -30.02240, -30.02230, -30.02053
##        speedKMH length turnType type      uuid                 endNode    speed
## 510000     6.29    226     NONE NONE 630554630                         1.747222
## 201070     6.14    233     NONE NONE 630556216                         1.705556
## 381300     6.39    326     NONE NONE 630552304        Av. Alcides Maia 1.775000
## 481060    15.63    595     NONE NONE 632032998 Av. Gen. Raphael Zippin 4.341667
## 491060     5.19    275     NONE NONE 631887650             R. Campinas 1.441667
## 511030     7.71    365     NONE NONE 631940708  R. Dr. Júlio Olszewski 2.141667
##        segments roadType delay                   street           pubMillis
## 510000     NULL        1    72          R. da República 2023-03-01 00:00:01
## 201070     NULL        1   103          R. Barão de Ubá 2023-03-01 00:00:01
## 381300     NULL        1   126           Av. Souza Melo 2023-03-01 00:00:00
## 481060     NULL        2    91             Av. Sertório 2023-03-01 00:59:01
## 491060     NULL        1   129                          2023-03-01 00:53:04
## 511030     NULL        2   132 R. Voluntários da Pátria 2023-03-01 00:55:27
##                hora       dia       date delay60 identificador
## 510000 1.157407e-05 Wednesday 2023-03-01     120             1
## 201070 1.157407e-05 Wednesday 2023-03-01     120             2
## 381300 0.000000e+00 Wednesday 2023-03-01     180             3
## 481060 4.098380e-02 Wednesday 2023-03-01     120             4
## 491060 3.685185e-02 Wednesday 2023-03-01     180             5
## 511030 3.850694e-02 Wednesday 2023-03-01     180             6

Qualificação da base de dados de filas

filas$hora = as.integer(substring(filas$pubMillis, 12, 13))
                    #filas@data[, 2:20] = NULL
filas = filas[(filas$hora  >= inicioPico) & (filas$hora  <= FimPico),]

#filas 16-18h = 103753

#filas 6-8h =  101169

Resumo das filas

resumo_filas_periodo = filas@data %>%
  summarise(Velocidade = mean(speedKMH), 
            VelocidadeSD = sd(speedKMH), 
            Atraso = mean(delay), 
            AtrasoSD = sd(delay), 
            Comprimento = mean(length), 
            ComprimentoSD = sd(length), 
            n = n())

resumo_filas_periodo
##   Velocidade VelocidadeSD   Atraso AtrasoSD Comprimento ComprimentoSD     n
## 1   10.31239     6.106539 148.6476  93.3187     606.977       532.382 51019
save(resumo_filas_periodo, file = paste("resumo_filas_periodo", MinDelay, MaxDelay, inicioPico,FimPico, ".Rda",sep = "_"))

Resumo das filas por data e hora

resumo_filas = filas@data %>%
  group_by(date, hora) %>%
  summarise(Velocidade = mean(speedKMH), 
            Atraso = mean(delay), 
            Comprimento = mean(length), 
            n = n())
## `summarise()` has grouped output by 'date'. You can override using the
## `.groups` argument.
# ATENÇÃO, considerando apenas horas com trinta ou mais ocorrências
resumo_filas = resumo_filas[resumo_filas$n >= 30,]
datatable(head(resumo_filas, 10), caption = htmltools::tags$caption(
  style = 'caption-side: top; text-align: left;',
  paste('Tabela', numero), htmltools::em(iconv(paste('Resumo de filas no período entre', inicioPico, "h e ",FimPico, "h.", sep = ""), from = "", to = "UTF-8"))
))
numero = numero +1
#resumo_filas = resumo_filas[(resumo_filas$hora  >= 0) & (resumo_filas$hora  < 24),]
nrow(resumo_filas)
## [1] 497
#resumo_filas = resumo_filas[(resumo_filas$hora  > 13) & (resumo_filas$hora  < 19),]


resumo_filas$hora = paste(resumo_filas$hora,":00:00", sep = "")
resumo_filas$Data =  as.POSIXlt(paste(resumo_filas$date, resumo_filas$hora, sep = " "))
attributes(resumo_filas$Data[1])
## $names
##  [1] "sec"    "min"    "hour"   "mday"   "mon"    "year"   "wday"   "yday"  
##  [9] "isdst"  "zone"   "gmtoff"
## 
## $balanced
## [1] TRUE
## 
## $class
## [1] "POSIXlt" "POSIXt" 
## 
## $tzone
## [1] ""    "-03" "-02"
hours <- sapply( seq_len( nrow( resumo_filas ) ) , function(x) resumo_filas$Data[x]$hour )
resumo_filas$hours <- hours

Criação de Time Series

library(zoo)
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
resumo_filas$date = NULL
resumo_filas$hora = NULL

#tt_resumo_filas = zoo(resumo_filas[,c("Velocidade", "n")], order.by=resumo_filas[,"Data"])

tt_resumo_filas = zoo(resumo_filas)
head(tt_resumo_filas)
##   Velocidade Atraso   Comprimento n   Data                hours
## 1 13.333590  140.7179  871.0256    39 2023-03-01 06:00:00 6    
## 2  9.643782  182.7949  687.5577   156 2023-03-01 07:00:00 7    
## 3 10.614607  132.4270  552.2584    89 2023-03-01 08:00:00 8    
## 4  9.630000  169.3333  646.3272   162 2023-03-02 07:00:00 7    
## 5  9.993942  133.2019  520.3462   104 2023-03-02 08:00:00 8    
## 6 12.270000  177.2727  926.0000    33 2023-03-03 06:00:00 6
class(resumo_filas)
## [1] "tbl_df"     "tbl"        "data.frame"
tryCatch(
  expr = {
    plot(tt_resumo_filas)
  },
  error = function(e){ 
    # (Optional)
    # Do this if an error is caught...
  }
  )
## Warning in xy.coords(x, y, xlabel, ylabel, log): NAs introduzidos por coerção
## Warning in min(x): nenhum argumento não faltante para min; retornando Inf
## Warning in max(x): nenhum argumento não faltante para max; retornando -Inf

## NULL
resumo_filas$ms = resumo_filas$Velocidade/3.6
resumo_filas$VelTipica = (resumo_filas$Comprimento/((resumo_filas$Comprimento/resumo_filas$ms)-resumo_filas$Atraso))*3.6
resumo_filas$DeltaVelocidade =  resumo_filas$VelTipica - resumo_filas$Velocidade

library(zoo)

tt_resumo_filas <- zoo(
  resumo_filas[, c("Velocidade", "Atraso", "VelTipica","Comprimento", "n")],
  order.by = as.POSIXct(resumo_filas$Data)
)

# Preparar dados
library(tidyr)
## 
## Attaching package: 'tidyr'
## The following object is masked from 'package:raster':
## 
##     extract
df_plot <- fortify.zoo(tt_resumo_filas) %>%
  mutate(Data = Index) %>%
  select(Data, Velocidade, VelTipica, Atraso, Comprimento) %>%
  pivot_longer(cols = -Data, names_to = "Variável", values_to = "Valor")

# Separar variáveis
df_primario <- df_plot %>% filter(Variável %in% c("Atraso", "Comprimento"))
df_secundario <- df_plot %>% filter(Variável %in% c("Velocidade", "VelTipica"))

# Definir limites personalizados
limite_esquerdo <- max(df_primario$Valor, na.rm = TRUE) * 1.05
limite_direito <- max(df_secundario$Valor, na.rm = TRUE) * 1.05
proporcao <- limite_esquerdo / limite_direito

# Reescalar as velocidades manualmente para caber no gráfico
df_secundario <- df_secundario %>%
  mutate(Valor_esc = Valor * proporcao)

# Plotar gráfico com dois eixos reais
library(ggplot2)
ggplot() +
  # Atraso e comprimento (eixo principal)
  geom_line(data = df_primario,
            aes(x = Data, y = Valor, color = Variável),
            size = 1) +
  
  # Velocidade e velocidade típica (eixo secundário, reescalado)
  geom_line(data = df_secundario,
            aes(x = Data, y = Valor_esc, color = Variável),
            size = 1) +
  
  scale_y_continuous(
    name = "Atraso e Comprimento",
    sec.axis = sec_axis(~ . / proporcao,
                        name = "Velocidade e Velocidade típica (km/h)")
  ) +
  labs(title = "Comparação Temporal das Variáveis de Mobilidade",
       x = "Tempo", y = "Valor", color = "Variável") +
  theme_minimal()
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Correlações

Traffic congestion level (0 = free flow 5 = blocked) RoadType
Value Type
1 Streets
2 Primary Street
3 Freeways
4 Ramps
5 Trails
6 Primary
7 Secondary
8, 14 4X4 Trails
15 Ferry crossing
9 Walkway
10 Pedestrian
11 Exit
16 Stairway
17 Private road
18 Railroads
19 Runway/Taxiway
20 Parking lot road
21 Service road

Hierarquia Funcional das Vias
Nível Tipo_de_via Funcao_principal
1️⃣ Freeways Rodovias expressas, alta velocidade, sem cruzamentos
2️⃣ Ramps Alcas de acesso as rodovias com controle de trafego
3️⃣ Primary Vias arteriais principais que conectam regioes
4️⃣ Primary Street Avenidas urbanas com alto fluxo e semaforos
5️⃣ Secondary Vias coletoras que distribuem trafego local
6️⃣ Streets Ruas locais com acesso direto a residencias e comercios
7️⃣ Trails Trilhas para pedestres ou ciclistas, nao veiculares 
paste("Consolidacao",inicioPico,FimPico,cel_size, DistanciaCentroides, ".Rda",sep = "_")
## [1] "Consolidacao_6_8_75_375_.Rda"
filas$VelTipica = (filas$length/((filas$length/filas$speed)-filas$delay))*3.6
filas$DeltaVelocidade =  filas$VelTipica - filas$speedKMH
filas$hierarquia = ifelse(filas$roadType==1, 6,
                   ifelse(filas$roadType==2, 4,
                   ifelse(filas$roadType==3, 1,
                   ifelse(filas$roadType==4, 2,
                   ifelse(filas$roadType==5, 7,
                   ifelse(filas$roadType==6, 3,
                   ifelse(filas$roadType==7, 5,filas$roadType)))))))


library(ggcorrplot)
#para_correlacao = resumo_filas[,c("Velocidade","Atraso","Comprimento", "n")]
para_correlacao = filas@data[,c("speedKMH","length", "delay","level", "hierarquia", "DeltaVelocidade")]

# matriz de correla??o
corr <- round(cor(para_correlacao), 1)
corr
##                 speedKMH length delay level hierarquia DeltaVelocidade
## speedKMH             1.0    0.7  -0.2  -0.7       -0.4             0.1
## length               0.7    1.0   0.4  -0.3       -0.4             0.2
## delay               -0.2    0.4   1.0   0.4       -0.2             0.2
## level               -0.7   -0.3   0.4   1.0        0.0             0.1
## hierarquia          -0.4   -0.4  -0.2   0.0        1.0            -0.3
## DeltaVelocidade      0.1    0.2   0.2   0.1       -0.3             1.0
datatable(corr)
# Obtendo os índices do triângulo inferior
indices <- which(lower.tri(corr), arr.ind = TRUE)
# Criando um dataframe inicial
df_cor <- data.frame(
  Linha = rownames(corr)[indices[, 1]],
  Coluna = colnames(corr)[indices[, 2]],
  Valor = corr[lower.tri(corr)]
)
# Transformando o formato para ter as combinações como colunas
library(tidyr)
df_final_correlacoes <- df_cor %>%
  pivot_wider(names_from = c(Linha, Coluna), values_from = Valor)

print(df_final_correlacoes)
## # A tibble: 1 × 15
##   length_speedKMH delay_speedKMH level_speedKMH hierarquia_speedKMH
##             <dbl>          <dbl>          <dbl>               <dbl>
## 1             0.7           -0.2           -0.7                -0.4
## # ℹ 11 more variables: DeltaVelocidade_speedKMH <dbl>, delay_length <dbl>,
## #   level_length <dbl>, hierarquia_length <dbl>, DeltaVelocidade_length <dbl>,
## #   level_delay <dbl>, hierarquia_delay <dbl>, DeltaVelocidade_delay <dbl>,
## #   hierarquia_level <dbl>, DeltaVelocidade_level <dbl>,
## #   DeltaVelocidade_hierarquia <dbl>
# p=valores
p.mat <- cor_pmat(para_correlacao)
p.mat
##                      speedKMH length delay         level   hierarquia
## speedKMH         0.000000e+00      0     0  0.000000e+00 0.000000e+00
## length           0.000000e+00      0     0  0.000000e+00 0.000000e+00
## delay            0.000000e+00      0     0  0.000000e+00 0.000000e+00
## level            0.000000e+00      0     0  0.000000e+00 1.041528e-18
## hierarquia       0.000000e+00      0     0  1.041528e-18 0.000000e+00
## DeltaVelocidade 1.040624e-140      0     0 7.692471e-207 0.000000e+00
##                 DeltaVelocidade
## speedKMH          1.040624e-140
## length             0.000000e+00
## delay              0.000000e+00
## level             7.692471e-207
## hierarquia         0.000000e+00
## DeltaVelocidade    0.000000e+00
ggcorrplot(corr, hc.order = TRUE, type = "lower",
           lab = TRUE)+ 
  labs(title = "Correlações",
       subtitle = paste("Período entre ", inicioPico,"h e ",FimPico,"h,", sep = ""))

HexPols_sf = st_as_sf(HexPols)
HexPols_sf = st_transform(HexPols_sf, 3857)
class(filas) 
## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"
class(HexPols_sf)
## [1] "sf"         "data.frame"
mapview(HexPols_sf, col.regions="magenta", color="magenta",
        layer.name = c("Hexágonos"))
## Warning in cbind(`Feature ID` = fid, mat): number of rows of result is not a
## multiple of vector length (arg 1)
filas_sf = st_as_sf(filas)

projection(HexPols_sf)
## [1] "+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs"
st_crs(filas_sf)   <- "+proj=longlat +datum=WGS84 +no_defs"
projection(HexPols_sf)
## [1] "+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs"
filas_sf = st_transform(filas_sf, 3857)
projection(filas_sf)
## [1] "+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs"
nrow(HexPols_sf)
## [1] 10192
length(HexPols)
## [1] 10192
HexPols$Ids = 1:length(HexPols)
HexPols_sf$Ids = 1:length(HexPols_sf$geometry)
head(HexPols)
##   Ids
## 1   1
## 2   2
## 3   3
## 4   4
## 5   5
## 6   6
tail(HexPols)
##         Ids
## 10187 10187
## 10188 10188
## 10189 10189
## 10190 10190
## 10191 10191
## 10192 10192
head(HexPols_sf)
## Simple feature collection with 6 features and 1 field
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -5686230 ymin: -3534680 xmax: -5685780 ymax: -3534593
## Projected CRS: WGS 84 / Pseudo-Mercator
##                         geometry Ids
## 1 POLYGON ((-5686230 -3534615...   1
## 2 POLYGON ((-5686155 -3534615...   2
## 3 POLYGON ((-5686080 -3534615...   3
## 4 POLYGON ((-5686005 -3534615...   4
## 5 POLYGON ((-5685930 -3534615...   5
## 6 POLYGON ((-5685855 -3534615...   6
tail(HexPols_sf)
## Simple feature collection with 6 features and 1 field
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -5689830 ymin: -3498581 xmax: -5689568 ymax: -3498298
## Projected CRS: WGS 84 / Pseudo-Mercator
##                             geometry   Ids
## 10187 POLYGON ((-5689793 -3498516... 10187
## 10188 POLYGON ((-5689830 -3498450... 10188
## 10189 POLYGON ((-5689755 -3498450... 10189
## 10190 POLYGON ((-5689718 -3498385... 10190
## 10191 POLYGON ((-5689643 -3498385... 10191
## 10192 POLYGON ((-5689680 -3498320... 10192
length(HexPols_sf$Ids)
## [1] 10192
length(unique(HexPols_sf$Ids))
## [1] 10192
res = st_intersection(HexPols_sf, filas_sf)
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries
length(res$Ids)
## [1] 40659
length(unique(res$Ids))
## [1] 1765
nrow(res)/nrow(filas_sf)*100
## [1] 79.69384
dim(res)
## [1] 40659    26
head(res)
## Simple feature collection with 6 features and 25 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -5703072 ymin: -3521219 xmax: -5688605 ymax: -3507162
## Projected CRS: WGS 84 / Pseudo-Mercator
##       Ids country         city level
## 4294 4294      BR Porto Alegre     2
## 7041 7041      BR Porto Alegre     3
## 2641 2641      BR Porto Alegre     4
## 5521 5521      BR Porto Alegre     4
## 1967 1967      BR Porto Alegre     4
## 3313 3313      BR Porto Alegre     2
##                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                line
## 4294 -51.11581, -51.11632, -51.11706, -51.11792, -51.11815, -51.11853, -51.11875, -51.12003, -51.12047, -51.12067, -51.12097, -51.12239, -51.12327, -51.12355, -51.12391, -51.12441, -51.12468, -51.12509, -51.12615, -51.12646, -51.12677, -51.12694, -51.12826, -51.12844, -51.12861, -51.12908, -51.12938, -51.13028, -51.13124, -51.13156, -51.13214, -51.13266, -51.13328, -51.13378, -51.13427, -51.13710, -51.13784, -51.14002, -51.14205, -51.14251, -51.14294, -51.14339, -51.14374, -51.14584, -51.14614, -51.14670, -51.14851, -51.14955, -30.07895, -30.07885, -30.07875, -30.07868, -30.07864, -30.07854, -30.07843, -30.07761, -30.07740, -30.07732, -30.07723, -30.07687, -30.07667, -30.07663, -30.07660, -30.07656, -30.07657, -30.07660, -30.07681, -30.07684, -30.07685, -30.07685, -30.07685, -30.07684, -30.07682, -30.07674, -30.07667, -30.07644, -30.07616, -30.07608, -30.07596, -30.07586, -30.07568, -30.07549, -30.07525, -30.07376, -30.07325, -30.07154, -30.06992, -30.06957, -30.06930, -30.06906, -30.06893, -30.06825, -30.06815, -30.06795, -30.06727, -30.06670
## 7041                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                           -51.16331, -51.16360, -51.16557, -51.16611, -51.16807, -51.16880, -51.16954, -51.16984, -51.17040, -30.02801, -30.02802, -30.02809, -30.02810, -30.02816, -30.02818, -30.02818, -30.02816, -30.02810
## 2641                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                       -51.17193, -51.17210, -51.17257, -51.17304, -51.17346, -51.17406, -51.17453, -51.17501, -51.17552, -51.17564, -51.17576, -51.17599, -51.17632, -51.17673, -51.17705, -30.11635, -30.11634, -30.11630, -30.11623, -30.11620, -30.11613, -30.11604, -30.11598, -30.11593, -30.11591, -30.11586, -30.11573, -30.11554, -30.11534, -30.11514
## 5521                                                                                                                                                                                                                                                                                                                                                                                                                                                         -51.09168, -51.09209, -51.09268, -51.09291, -51.09399, -51.09432, -51.09460, -51.09477, -51.09495, -51.09610, -51.09635, -51.09704, -51.09777, -51.09801, -51.09826, -51.09848, -51.09892, -51.09920, -51.09933, -51.09953, -51.09978, -51.10026, -51.10043, -51.10061, -51.10083, -51.10102, -51.10127, -51.10161, -30.04682, -30.04705, -30.04736, -30.04745, -30.04773, -30.04778, -30.04779, -30.04778, -30.04775, -30.04740, -30.04731, -30.04701, -30.04667, -30.04657, -30.04644, -30.04633, -30.04610, -30.04597, -30.04594, -30.04590, -30.04589, -30.04600, -30.04602, -30.04603, -30.04602, -30.04600, -30.04593, -30.04581
## 1967                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                   -51.17312, -51.17295, -51.17287, -51.17282, -51.17231, -30.13898, -30.13888, -30.13878, -30.13860, -30.13736
## 3313                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                       -51.22914, -51.22923, -51.22943, -51.22954, -51.22964, -51.23010, -51.23021, -51.23036, -51.23087, -51.23095, -51.23098, -51.23101, -51.23129, -51.23155, -51.23157, -30.10616, -30.10605, -30.10575, -30.10556, -30.10534, -30.10393, -30.10368, -30.10338, -30.10249, -30.10226, -30.10213, -30.10199, -30.10005, -30.09818, -30.09805
##      speedKMH length turnType type      uuid         endNode    speed segments
## 4294    19.47   3640     NONE NONE 651986012                 5.408333     NULL
## 7041    12.76    685     NONE NONE 652831980                 3.544444     NULL
## 2641     5.63    520     NONE NONE 652186202 Est. Costa Gama 1.563889     NULL
## 5521     8.38   1035     NONE NONE 650312522            <NA> 2.327778     NULL
## 1967     3.98    200     NONE NONE 652334058 Est. Costa Gama 1.105556     NULL
## 3313    20.14    939     NONE NONE 653016672                 5.594444     NULL
##      roadType delay                        street           pubMillis hora
## 4294        2   423           Av. Bento Gonçalves 2023-03-01 06:40:12    6
## 7041        2   138          Av. Dr. Nilo Peçanha 2023-03-01 06:53:25    6
## 2641        2   283 Est. Afonso Lourenço Mariante 2023-03-01 06:43:44    6
## 5521        2   366            Av. Protásio Alves 2023-03-01 06:12:05    6
## 1967        2   157             Est. Gedeon Leite 2023-03-01 06:46:10    6
## 3313        2    95              Av. da Cavalhada 2023-03-01 06:55:48    6
##            dia       date delay60 identificador VelTipica DeltaVelocidade
## 4294 Wednesday 2023-03-01     480            20  52.40857        32.93857
## 7041 Wednesday 2023-03-01     180            21  44.62526        31.86526
## 2641 Wednesday 2023-03-01     300            22  37.81479        32.18479
## 5521 Wednesday 2023-03-01     420            23  47.38645        39.00645
## 1967 Wednesday 2023-03-01     180            24  30.11982        26.13982
## 3313 Wednesday 2023-03-01     120            25  46.40533        26.26533
##      hierarquia                  geometry
## 4294          4 POINT (-5693942 -3512127)
## 7041          4 POINT (-5696262 -3507162)
## 2641          4 POINT (-5697003 -3518359)
## 5521          4 POINT (-5688605 -3509440)
## 1967          4 POINT (-5696475 -3521219)
## 3313          4 POINT (-5703072 -3516160)
pontos_por_hex = as.data.frame(table(res$Id))

res_sub = res[res$Ids == 7886, ]
#mapview(res_sub)
library(tidyr)
res_summ = res %>%
  group_by(Ids) %>%
  summarise(speed = round(mean(speedKMH),2), 
            length = round(mean(length),2),
            roadType = mean(roadType),
            delay = round(mean(delay),2),
            Contagem = n())
res_summ_road_type = res %>%
  count(roadType)
res_summ_road_type$geometry = NULL
names(res_summ_road_type) = c("roadType", "Freq")

datatable(head(res_summ, 10), caption = htmltools::tags$caption(
  style = 'caption-side: top; text-align: left;',
  paste('Tabela', numero), htmltools::em('Resumo Polígonos 24h.')
))
numero = numero +1
#mapview(HexPols_sf[HexPols$Ids==9242,]) + mapview(res[res$Ids==9242,], col.regions="magenta")
datatable(res_summ[2234:2250,])
restemp = res
restemp$geometry = NULL

# https://link.springer.com/article/10.1140/epjs/s11734-021-00424-2
# (0 = free flow, 1 = Low, 2 = Medium, 3 = High, 4 = Very High, 5 = closed

res_summ_table = restemp %>%
  group_by(Ids, level) %>%
  summarise(Contagem = n(), 
            .groups = "drop")%>%
  pivot_wider(names_from = level, values_from = Contagem, values_fill = 0) 

names(res_summ_table) = c("Ids",paste("L", names(res_summ_table[,2:ncol(res_summ_table)]), sep = "_"))

Lista_para_renomear <- list(
  L_0 = "Free",
  L_1 = "Low",
  L_2 = 'Medium',
  L_3 = "High",
  L_4 = "Very High",
  L_5 = "Closed"
)

rename_df_from_ls <- function(df, l){
  l <- l[names(l) %in% names(df)]
  if (length(l)){
    rename_with(df, ~unlist(l), all_of(names(l)))
  } else {
    df
  }
}

res_summ_table = res_summ_table %>% 
  rename_df_from_ls(Lista_para_renomear)

rm(restemp)

datatable(res_summ_table[2234:2250,], caption = htmltools::tags$caption(
  style = 'caption-side: top; text-align: left;',
  paste('Tabela', numero), htmltools::em('Filas no período de 24h com  nível de criticidade.')
))
numero = numero +1
length(HexPols)
## [1] 10192
length(unique(res_summ$Ids))
## [1] 1765
100*(length(unique(res_summ$Ids))/length(HexPols_sf$Ids))
## [1] 17.3175
length(unique(res_summ_table$Ids))
## [1] 1765
res_summ = merge(x = res_summ, y = res_summ_table, by = "Ids", all.x = TRUE)
res_summ$Eventos = res_summ$`Very High`+res_summ$Medium+res_summ$High+res_summ$Low
glimpse(res_summ)
## Rows: 1,765
## Columns: 12
## $ Ids         <int> 31, 155, 161, 211, 335, 394, 413, 427, 433, 441, 458, 463,…
## $ speed       <dbl> 24.93, 23.09, 28.25, 13.50, 13.60, 10.20, 7.18, 9.52, 21.4…
## $ length      <dbl> 1994.0, 831.0, 1105.0, 914.0, 909.0, 2435.0, 406.0, 381.0,…
## $ roadType    <dbl> 2.0, 2.0, 7.0, 7.0, 1.0, 1.0, 2.0, 2.0, 2.0, 2.0, 2.0, 1.0…
## $ delay       <dbl> 114.0, 67.0, 62.0, 174.0, 79.0, 370.0, 151.0, 118.0, 67.0,…
## $ Contagem    <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 4, 2, 1, 1, 1, 1, 1, 1…
## $ Medium      <int> 1, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0…
## $ Low         <int> 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 2, 1, 0, 0, 0, 1, 0…
## $ High        <int> 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 3, 0, 0, 1, 1, 0, 0, 1…
## $ `Very High` <int> 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0…
## $ geometry    <GEOMETRY [m]> POINT (-5686869 -3534223), POINT (-5688316 -35328…
## $ Eventos     <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 4, 2, 1, 1, 1, 1, 1, 1…
save(res_summ, file = "res_summ.Rda")
save(res_summ_table, file = "res_summ_table.Rda")
# left join HexPols_sf res_summ
a = as.data.frame(HexPols_sf$Ids)
names(a) = "Ids"
b = res_summ
b$geometry = NULL
binds  <- left_join(a, b, by=c('Ids'='Ids'))
HexPols_sf$Eventos = binds$Eventos
HexPols_sf$speed = binds$speed
HexPols_sf$length = binds$length
HexPols_sf$r
## NULL
length(HexPols)
## [1] 10192
Dfhigh <- res_summ[order(-res_summ$Contagem),]
Dfhigh = Dfhigh$Ids[1:10]
deletar = res[res$Ids %in% Dfhigh, c(1:4,6:12,14:22)]
library(RColorBrewer)
mapview(HexPols_sf[c(Dfhigh), ], zcol = "Eventos",alpha.regions=1,
        col.regions=brewer.pal(9, "YlOrRd"),
        layer.name = paste("Polígonos mais críticos no período entre ", inicioPico, "h e ",FimPico, "h.", sep = ""))
## Warning: Found less unique colors (9) than unique zcol values (10)! 
## Interpolating color vector to match number of zcol values.
#mapview(deletar, cex=5, col.regions="red", color="magenta", alpha.regions=1, layer.name = c("Polígonos mais críticos período de 24h"))
library(leaflet)
library(leaflet.extras)
res_spatial = filas[filas$identificador %in% res$identificador, ]
class(res_spatial)
## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"
leaflet(res_spatial) %>% 
  addProviderTiles("CartoDB.Positron") %>% 
  #addCircles(color = "green") %>% 
  addHeatmap(blur= 1, max = 1, radius = 2)
#HexPols_sf  <- left_join(HexPols_sf, as.data.frame(res_summ), by=c('Ids'='Ids'))
library(RColorBrewer)
#HexPols_sf$Eventos = ifelse(is.na(HexPols_sf$Eventos), 0, HexPols_sf$Eventos)

#mapview(HexPols_sf, color="magenta", zcol="Eventos", col.regions=brewer.pal(9, "YlOrRd"), alpha.regions=1, alpha=0, na.color="#333333",         layer.name = c("Polígonos com eventos no período de 24h"))
mapview(HexPols_sf[c(Dfhigh), ], color="magenta", zcol="Eventos", col.regions=brewer.pal(9, "YlOrRd"), alpha.regions=1, alpha=0, na.color="#FFFFFF",
        layer.name = paste("Dez polígonos com mais eventos no período entre ", inicioPico, "h e ",FimPico, "h.", sep = ""))
## Warning: Found less unique colors (9) than unique zcol values (10)! 
## Interpolating color vector to match number of zcol values.
class(HexPols)
## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"
#mapview(HexPols)
head(HexPols)
##   Ids
## 1   1
## 2   2
## 3   3
## 4   4
## 5   5
## 6   6
head(HexPols)
##   Ids
## 1   1
## 2   2
## 3   3
## 4   4
## 5   5
## 6   6
HexPols@data  <- left_join(HexPols@data, as.data.frame(res_summ), by=c('Ids'='Ids'))
#mapview(HexPols, zcol="Eventos", alpha.regions=1, alpha=0)
mapview(HexPols, zcol="Eventos", col.regions=brewer.pal(9, "YlOrRd"), at = seq(0, max(HexPols$Eventos, na.rm = 1), 50), legend = TRUE, alpha.regions=1, stroke = FALSE)
## Warning in st_as_sf.Spatial(x): column "geometry" will be overwritten by
## geometry column
## Warning: Found less unique colors (9) than unique zcol values (11)! 
## Interpolating color vector to match number of zcol values.
HexPols$Eventos = ifelse(is.na(HexPols$Contagem), 0, HexPols$Contagem)

Análise de correlação espacial

### analise de correlação espacial 
require(maptools)      
## Carregando pacotes exigidos: maptools
## Please note that 'maptools' will be retired during October 2023,
## plan transition at your earliest convenience (see
## https://r-spatial.org/r/2023/05/15/evolution4.html and earlier blogs
## for guidance);some functionality will be moved to 'sp'.
##  Checking rgeos availability: TRUE
## 
## Attaching package: 'maptools'
## The following object is masked from 'package:sp':
## 
##     sp2Mondrian
require(sp)            
require(spdep)        
## Carregando pacotes exigidos: spdep
## Carregando pacotes exigidos: spData
## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`
require(classInt)     
## Carregando pacotes exigidos: classInt
require(RColorBrewer)
## Visualizando um mapa do atributo: Quantis
# HexPols = HexPols[HexPols$Eventos >0, ]
INT1 <- classIntervals(HexPols$Eventos, n=4, style="quantile")
CORES.1 <- c(rev(brewer.pal(3, "Blues")), brewer.pal(3, "Reds"))[-1] 
COL1 <- findColours(INT1, CORES.1)
# plot(HexPols, col=COL1, axes = FALSE, type = "n", border = NA)
# title("Ocorrencias: Quantis")
# TB1 <- attr(COL1, "table")
# legtext <- paste(names(TB1))
# legend("bottomright", fill=attr(COL1, "palette"), legend=legtext, bty="n",cex=0.8)
## Calculando a matriz de vizinhan?as
ccods = coordinates(HexPols)
points = cbind(ccods[,1],ccods[,2])
dnb = dnearneigh(points,0,DistanciaCentroides) # distancia entre centroides
W.Bin= nb2mat(neighbours = dnb, style = "B") #Converte a lista de vizinhos em uma matriz binária (1 para vizinhos, 0 para não vizinhos).
W.Normal= nb2mat(neighbours = dnb, style = "W") #Quando nb2mat() usa style="W", ele gera uma matriz de vizinhança ponderada, onde os pesos são atribuídos de forma inversamente proporcional ao número de vizinhos de cada polígono. Em outras palavras, cada entrada da matriz representa 1 / número de vizinhos do polígono correspondente. Isso significa que:
# - Se um polígono tem 3 vizinhos, cada um recebe um peso de 1/3.
# - Se um polígono tem 5 vizinhos, cada um recebe um peso de 1/5.
# - Quanto mais vizinhos um polígono tiver, menor será o peso atribuído a cada um.
## Calculando o ?ndice de Moran Global

# Pelo teste de Normalidade
moran.test(HexPols$Eventos,listw=nb2listw(dnb, style = "W"), randomisation= FALSE) 
## 
##  Moran I test under normality
## 
## data:  HexPols$Eventos  
## weights: nb2listw(dnb, style = "W")    
## 
## Moran I statistic standard deviate = 8.6367, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      2.661497e-02     -9.812580e-05      9.566464e-06
# Pelo teste de Permuta??o
moran.test(HexPols$Eventos,listw=nb2listw(dnb, style = "W"), randomisation= TRUE)
## 
##  Moran I test under randomisation
## 
## data:  HexPols$Eventos  
## weights: nb2listw(dnb, style = "W")    
## 
## Moran I statistic standard deviate = 8.6989, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      2.661497e-02     -9.812580e-05      9.430215e-06
# Por simula??o de Monte-Carlo 
moran.mc(HexPols$Eventos, listw=nb2listw(dnb, style = "W"), nsim=999)
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  HexPols$Eventos 
## weights: nb2listw(dnb, style = "W")  
## number of simulations + 1: 1000 
## 
## statistic = 0.026615, observed rank = 1000, p-value = 0.001
## alternative hypothesis: greater
## Moran Local
ShapePB.mloc <- localmoran(HexPols$Eventos, listw=nb2listw(dnb, style="W")) 
head(ShapePB.mloc)
##           Ii          E.Ii      Var.Ii      Z.Ii Pr(z != E(Ii))
## 1 0.02888918 -2.834773e-06 0.002405062 0.5891348      0.5557709
## 2 0.02888918 -2.834773e-06 0.002219839 0.6132210      0.5397302
## 3 0.02888918 -2.834773e-06 0.002405062 0.5891348      0.5557709
## 4 0.02888918 -2.834773e-06 0.002405062 0.5891348      0.5557709
## 5 0.02888918 -2.834773e-06 0.002623961 0.5640259      0.5727365
## 6 0.02888918 -2.834773e-06 0.002623961 0.5640259      0.5727365
## Mapa das probabilidades (Signific?ncias do I de Moral Local)
INT4 <- classIntervals(ShapePB.mloc[,5], style="fixed", 
                       fixedBreaks=c(0,0.01, 0.05, 0.10))
## Warning in classIntervals(ShapePB.mloc[, 5], style = "fixed", fixedBreaks =
## c(0, : variable range greater than fixedBreaks
CORES.4 <- c(rev(brewer.pal(3, "Reds")), brewer.pal(3, "Blues"))
COL4 <- findColours(INT4, CORES.4)
# plot(HexPols, col=COL4, axes = FALSE, type = "n", border = NA)
# title("P-valores do I de Moran Local")
# 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)
## LISA Map (+/+) , (-,-), (+,-), (-,+)

# Montando matrix W de vizinhan?a
ShapeCG.nb1.mat <- nb2mat(dnb)

# Incid?ncia de dengue padronizada
Dengue_SD <- scale(HexPols$Eventos)

# M?dia das incid?ncias de dengue padronizada
Dengue_W <- ShapeCG.nb1.mat %*% Dengue_SD
# Diagrama de espalhamento de Moran
plot(Dengue_SD, Dengue_W,xlab="Z",ylab="WZ")
abline(v=0, h=0)
title("Diagrama de Espalhamento de Moran")

# Montando o Mapa LISA
Q <- vector(mode = "numeric", length = nrow(ShapePB.mloc))
Q[(Dengue_SD>0  & Dengue_W > 0)] <- 1            
Q[(Dengue_SD<0  & Dengue_W < 0)] <- 2
Q[(Dengue_SD>=0 & Dengue_W < 0)] <- 3
Q[(Dengue_SD<0  & Dengue_W >= 0)]<- 4
signif=0.05
            #signif=1
HexPols$Q = Q
Q[ShapePB.mloc[,5]>signif]<-5
df = data.frame()

CORES.5 <- c("blue", "green" , "red", "yellow", "gray", rgb(0.95,0.95,0.95))
HexPols$cores5Q = CORES.5[Q]
HexPols$cores5 =  ifelse(HexPols$cores5Q == "blue", "A-A", ifelse(HexPols$cores5Q == "green", "B-B", 
                                                                   ifelse(HexPols$ShapeCG == "red", "A-B", ifelse(HexPols$cores5Q == "yellow", "B-A", "N-S"))))
df = as.data.frame(cbind(Q, CORES.5[Q]))
names(df)= c("Q", "cores5Q")
df$Incid = HexPols$Eventos
df$pvalor = ShapePB.mloc[,5]
 # plot(HexPols, col=CORES.5[Q], axes = FALSE, type = "n", border = NA)
 # title("Mapa LISA")
 # legend("bottomright", c("Q1(+/+)", "Q2(-/-)", "Q3(+/-)", "Q4(-/+)","NS"), fill=CORES.5)
HexPols$pvalor = ShapePB.mloc[,5]
library(mapview)
yourNewPolygon<- HexPols
yourNewPolygon@data$cores5 = ifelse(is.na(yourNewPolygon@data$cores5),"NS",yourNewPolygon@data$cores5)
mapview(yourNewPolygon, zcol = "cores5", col.regions=c("red", "orange", "green", "yellow", "grey"), layer.name = c("Moran's Index"))
## Warning in st_as_sf.Spatial(x): column "geometry" will be overwritten by
## geometry column
# mapview(yourNewPolygon, zcol = "cores5Q", col.regions=c("red", "orange", "green", "yellow", "grey"),
#         layer.name = c("Moran's Index"))
mapview(yourNewPolygon[yourNewPolygon$Q==1,], zcol = "cores5", col.regions=c("red", "orange", "green", "yellow", "grey"),
        layer.name = c("Moran's Index"))
## Warning in st_as_sf.Spatial(x): column "geometry" will be overwritten by
## geometry column
# mapview(yourNewPolygon[yourNewPolygon$Q==1,], zcol = "cores5Q", col.regions=c("red", "orange", "green", "yellow", "grey"),
#         layer.name = c("Moran's Index"))
unique(HexPols$cores5Q)
## [1] "gray"   "yellow" "blue"
# Numero de congestionamentos no ano
Filas = length(filas)
Filas
## [1] 51019
# Numero de congestionamentos nas vias monitoradas
CongestionamentosViasMonitoradas= sum(yourNewPolygon$Eventos)
CongestionamentosViasMonitoradas
## [1] 40659
# Percentual de congestionamentos nas vias monitoradas
PercentualCongestionamentosViasMonitoradas = sum(yourNewPolygon$Eventos)/length(filas)*100
PercentualCongestionamentosViasMonitoradas
## [1] 79.69384
blue = yourNewPolygon[yourNewPolygon$cores5Q == "blue", c("Eventos")]
# N?mero de poligonos nas vias monitoradas
PolihonosMonitorados = length(yourNewPolygon$Eventos)
PolihonosMonitorados
## [1] 10192
# N?mero de  poligonos com congestionamentos nas vias monitoradas
PoligonosCongestionados =  length(res_summ$Ids)
PoligonosCongestionados
## [1] 1765
# N?mero de  poligonos comcongestionamentos significativos nas vias monitoradas
PoligonosCongestionadosSignificativos = nrow(yourNewPolygon[yourNewPolygon$cores5Q=="blue", ])
PoligonosCongestionadosSignificativos
## [1] 151
# Percentual de  poligonos com congestionamentos significativos nas vias monitoradas
PercentualPoligonosCongestionadosSignificativos = (nrow(yourNewPolygon[yourNewPolygon$cores5Q=="blue", ])/length(yourNewPolygon$Eventos))*100
PercentualPoligonosCongestionadosSignificativos
## [1] 1.481554
# Percentual de congestionamentos registrados nos pol?gonos significativos em rela??o aos comngestionamentos das vias monitoradas
PercentCongNosPoligonosSignificativos =  sum(blue$Eventos)/sum(yourNewPolygon$Eventos)*100
PercentCongNosPoligonosSignificativos
## [1] 20.54158
# Percentual de congestionamentos registrados nos pol?gonos significativos em rela??o aos comngestionamentos totais
PercentCongNosPoligonosSignificativosRelacaoAoTotal = sum(blue$Eventos)/length(filas)*100
PercentCongNosPoligonosSignificativosRelacaoAoTotal
## [1] 16.37037
periodo = c(inicioPico,FimPico,cel_size, DistanciaCentroides, Filas,
            CongestionamentosViasMonitoradas, PercentualCongestionamentosViasMonitoradas,
            PolihonosMonitorados, PoligonosCongestionados, PoligonosCongestionadosSignificativos,
            PercentualPoligonosCongestionadosSignificativos, PercentCongNosPoligonosSignificativos,
            PercentCongNosPoligonosSignificativosRelacaoAoTotal, resumo_filas_periodo, df_final_correlacoes, res_summ_road_type)
periodo
## [[1]]
## [1] 6
## 
## [[2]]
## [1] 8
## 
## [[3]]
## [1] 75
## 
## [[4]]
## [1] 375
## 
## [[5]]
## [1] 51019
## 
## [[6]]
## [1] 40659
## 
## [[7]]
## [1] 79.69384
## 
## [[8]]
## [1] 10192
## 
## [[9]]
## [1] 1765
## 
## [[10]]
## [1] 151
## 
## [[11]]
## [1] 1.481554
## 
## [[12]]
## [1] 20.54158
## 
## [[13]]
## [1] 16.37037
## 
## $Velocidade
## [1] 10.31239
## 
## $VelocidadeSD
## [1] 6.106539
## 
## $Atraso
## [1] 148.6476
## 
## $AtrasoSD
## [1] 93.3187
## 
## $Comprimento
## [1] 606.977
## 
## $ComprimentoSD
## [1] 532.382
## 
## $n
## [1] 51019
## 
## $length_speedKMH
## [1] 0.7
## 
## $delay_speedKMH
## [1] -0.2
## 
## $level_speedKMH
## [1] -0.7
## 
## $hierarquia_speedKMH
## [1] -0.4
## 
## $DeltaVelocidade_speedKMH
## [1] 0.1
## 
## $delay_length
## [1] 0.4
## 
## $level_length
## [1] -0.3
## 
## $hierarquia_length
## [1] -0.4
## 
## $DeltaVelocidade_length
## [1] 0.2
## 
## $level_delay
## [1] 0.4
## 
## $hierarquia_delay
## [1] -0.2
## 
## $DeltaVelocidade_delay
## [1] 0.2
## 
## $hierarquia_level
## [1] 0
## 
## $DeltaVelocidade_level
## [1] 0.1
## 
## $DeltaVelocidade_hierarquia
## [1] -0.3
## 
## $roadType
## [1]  1  2  3  4  6  7 20
## 
## $Freq
## [1] 11005 29373     2   230    44     2     3
save(periodo, file = paste("Consolidacao", MinDelay, MaxDelay, inicioPico,FimPico,cel_size, DistanciaCentroides, ".Rda",sep = "_"))
significativos = yourNewPolygon[yourNewPolygon$Q==1,]
# mapview(significativos, zcol = "Eventos", col.regions=rev(c("red", "orange", "green", "yellow", "grey")),
#         layer.name = c("Ocorrências"), alpha.regions = 1)
# mapview(significativos, zcol = "speed", col.regions=(c("red", "orange", "green", "yellow", "grey")),
#         layer.name = c("Velocidade Média"), alpha.regions = 1)
# mapview(significativos, zcol = "length", col.regions=rev(c("red", "orange", "green", "yellow", "grey")),
#         layer.name = c("Fila Média"), alpha.regions = 1)
mapview(significativos, zcol = "delay", col.regions=rev(c("red", "orange", "green", "yellow", "grey")),
        layer.name = c("Atraso Médio"), alpha.regions = 1)
## Warning in st_as_sf.Spatial(x): column "geometry" will be overwritten by
## geometry column
## Warning: Found less unique colors (5) than unique zcol values (508)! 
## Interpolating color vector to match number of zcol values.
# Definição dos pesos
pesos <- c(1, 3, 5, 10) # Medium, Low, High, Very High
dados = significativos@data[ , c("Low", "Medium", "High", "Very High")]
significativos@data$Score <- rowSums(significativos@data[ , c("Low", "Medium", "High", "Very High")] * pesos)
# mapview(significativos, zcol = "Score", col.regions=rev(c("red", "orange", "green", "yellow", "grey")),
#         layer.name = c("Score"), alpha.regions = 1)
significativos@data$Score_Normalizado <- (significativos@data$Score / max(significativos@data$Score)) * 100
mapview(significativos, zcol = "Score_Normalizado", col.regions=rev(c("red", "orange", "green", "yellow", "grey")),
        layer.name = c("Score Normalizado"), alpha.regions = 1)
## Warning in st_as_sf.Spatial(x): column "geometry" will be overwritten by
## geometry column
## Warning: Found less unique colors (5) than unique zcol values (281)! 
## Interpolating color vector to match number of zcol values.