exploring.Rmd

---
title: "Exploring OD2017 data"
output: github_document
---

```{r, include=FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  echo = FALSE
)

library(tidyverse)
library(piggyback)
library(sf)
library(tmap)
library(leaflet)
library(stplanr)
library(geobr)
library(od)
library(osrm)
# assume data is in the folder
 folder = "/home/lucas/"
# folder = "~/Dropbox/DATA/"
# folder = "~/hd/data/project-data/spshared/"

# pb_download(repo = "spstreets/OD2017")

# folder = "~/hd/data/project-data/spshared/"
# load(file.path(folder, "OD2017.rda"))
OD_2017_v1 = readRDS("OD_2017_v1.Rds")

```

I create this separated document to be the main reference on the exploratory data analysis. The data set is at the trip level and I take advantage of the survey design and use the expansion factor (for the trips) in all the graphs and analysis here. The survey implies that there are `r sum(OD_2017_v1$fe_via, na.rm=TRUE)` trips/day in the São Paulo Metro Region.

### Key variables

The data set contains a lot of socioeconomic variables. Since our interest relies on trips, flows and mode choices, the analyses so far are at the trip level. In this context, the variables that contain a the modes used in each trip are the main ones (`modo1`, `modo2`, `modo3`, and `modo4`). There is also the coordinates for each residence (`co_dom_x` and `co_dom_y`) and workplaces (`co_tr1_x` and `co_tr1_y`), and obviously the origin and destination zones and coordinates (e.g. `zona_o`, `co_o_x`, `co_o_y`) and euclidean distance. As you can see in the graph below, the vast majority of trips are made using just one mode. I use the main mode used in the trip (`modoprin`) in the remaining graphs. 

```{r, echo=FALSE, eval=FALSE}

OD_2017_v1 = OD_2017_v1 %>%
  rowwise() %>%
  mutate(
    n_modes = n_distinct(c(modo1, modo2, modo3, modo4), na.rm = TRUE)
  ) %>%
  ungroup()
  

OD_2017_v1 %>%
  filter(!is.na(modoprin)) %>%
  ggplot() +
  geom_bar(aes(n_modes, weight=fe_via)) +
  scale_y_continuous(name = "Milhares de viagens / dia",
                     labels=function(x) format(x/1000,
                                               big.mark = ",",
                                               decimal.mark = ".",
                                               scientific = FALSE)
                     ) +
  ggtitle("# de modos usados em cada viagem") +
  xlab("") +
  theme_bw()

```

### Recategorisation of modes

I tried two different recategorisations (see code for the choices).The first is based on the "vehicle" itself, and the other relies on the "role" of the individual while in the trip.

```{r, echo=FALSE}

OD_2017_v1 = OD_2017_v1 %>%
  dplyr::mutate(
    modo_simples = case_when(
      modoprin %in% 1:6 ~ "Público",
      modoprin %in% 9:14 ~ "Individual motorizado",
      modoprin == 15 | modoprin == 16 ~ "Ativo",
      modoprin %in% c(7, 8, 17) ~ "Outro"
    )
  )

OD_2017_v1 %>%
  filter(!is.na(modoprin)) %>%
  ggplot() +
  geom_bar(aes(modo_simples, weight=fe_via)) +
  scale_y_continuous(name = "Milhares de viagens / dia",
                     labels=function(x) format(x/1000,
                                               big.mark = ",",
                                               decimal.mark = ".",
                                               scientific = FALSE)
                     ) +
  ggtitle("Distribuição das viagens por modo") +
  xlab("") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 60, vjust = 1, hjust = 1))

```

Now the A/B Streets classification (foot, bike, car, public transport and other):

```{r, echo=FALSE}

OD_2017_v1 = OD_2017_v1 %>%
  dplyr::mutate(
    mode_ab_streets = case_when(
      modoprin %in% 1:6 ~ "public",
      modoprin == 16 ~ "foot",
      modoprin==9 | modoprin==10 ~ "car",
      modoprin %in% c(7, 8, 11, 12, 13, 14, 17) ~ "other",
      modoprin==15 ~ "bike"
    )
  )

OD_2017_v1 %>%
  filter(!is.na(modoprin)) %>%
  ggplot() +
  geom_bar(aes(mode_ab_streets, weight=fe_via)) +
  scale_y_continuous(name = "Milhares de viagens / dia",
                     labels=function(x) format(x/1000,
                                               big.mark = ",",
                                               decimal.mark = ".",
                                               scientific = FALSE)
                     ) +
  ggtitle("Classificação A/B Streets") +
  xlab("") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 60, vjust = 1, hjust = 1))

```

Now the mode share by distance (using the first recategorisation). The average trips' distance is `r weighted.mean(OD_2017_v1$distancia, OD_2017_v1$fe_via, na.rm = TRUE)`.

```{r, echo=FALSE}

col_modes = c("#fe5f55", "grey", "#ffd166", "#457b9d")

OD_2017_v1 %>%
  filter(!is.na(modoprin)) %>%
  mutate(
    dist_km = distancia/1000,
    dist_bands = cut(dist_km,
                     breaks = c(0, 1, 3, 6, 10, 15, 20, 30, 83),     # max is 82.7
                     include.lowest = TRUE)
  ) %>%
  select(dist_bands, modo_simples, fe_via) %>%
  group_by(dist_bands, modo_simples) %>%
  summarise(Trips = sum(fe_via)) %>%
  mutate(modo_simples = factor(modo_simples, levels = c("Individual motorizado",
                                                        "Outro",
                                                        "Público",
                                                        "Ativo")
                               )
         ) %>%
  ggplot(aes(x = dist_bands, y = Trips, fill = modo_simples)) +
  geom_col() +
  xlab("Distância (km)") +
  scale_y_continuous(name = "Milhares de viagens / dia",
                     labels=function(x) format(x/1000,
                                               big.mark = ",",
                                               decimal.mark = ".",
                                               scientific = FALSE)) +
  scale_fill_manual(values = col_modes, name = "Modo") +
  theme_bw()

```


```{r, echo=FALSE}

col_modes = c("#fe5f55", "grey", "#ffd166", "#90be6d", "#457b9d")

OD_2017_v1 %>%
  filter(!is.na(modoprin)) %>%
  mutate(
    dist_km = distancia/1000,
    dist_bands = cut(dist_km,
                     breaks = c(0, 1, 3, 6, 10, 15, 20, 30, 83),     # max is 82.7
                     include.lowest = TRUE)
    ) %>%
  select(dist_bands, mode_ab_streets, fe_via) %>%
  group_by(dist_bands, mode_ab_streets) %>%
  summarise(Trips = sum(fe_via)) %>%
  mutate(mode_ab_streets = factor(mode_ab_streets, levels = c("car",
                                                              "other",
                                                              "public",
                                                              "bike",
                                                              "foot")
                                  )
         ) %>%
  ggplot(aes(x = dist_bands, y = Trips, fill = mode_ab_streets)) +
  geom_col() +
  xlab("Distance (km)") +
  scale_y_continuous(name = "Thousands of travels / day",
                     labels=function(x) format(x/1000,
                                               big.mark = ",",
                                               decimal.mark = ".",
                                               scientific = FALSE)) +
  scale_fill_manual(values = col_modes, name = "Mode") +
  theme_bw()

```


```{r, echo=FALSE}

OD_2017_v1 %>%
  filter(!is.na(modoprin)) %>%
  mutate(
    dist_km = distancia/1000,
    dist_bands = cut(dist_km,
                     breaks = c(0, 1, 3, 6, 10, 15, 20, 30, 83),     # max is 82.7
                     include.lowest = TRUE),
    motivo = case_when(
      motivo_o == 8 ~ motivo_d,
      motivo_d == 8 ~ motivo_o,
      (motivo_o != 8 & motivo_d != 8) ~ motivo_d
    ),
    motivo_simples = case_when(
      motivo %in% 1:3 ~ "Trabalho",
      motivo == 4 ~ "Escola",
      TRUE ~ "Outro")
    ) %>%
  select(dist_bands, motivo_simples, fe_via) %>%
  group_by(dist_bands, motivo_simples) %>%
  summarise(Trips = sum(fe_via)) %>%
  ggplot(aes(x = dist_bands, y = Trips, fill = motivo_simples)) +
  geom_col() +
  ggtitle("") +
  xlab("Distância (km)") +
  scale_y_continuous(name = "Milhares de viagens / dia",
                     labels=function(x) format(x/1000,
                                               big.mark = ",",
                                               decimal.mark = ".",
                                               scientific = FALSE)) +
  scale_fill_brewer(palette = "Set1", name = "Motivo") +
  theme_bw()

```


```{r, echo=FALSE}

OD_2017_v1 = OD_2017_v1 %>%
  mutate(
    saida = h_saida + min_saida/60,
    chegada = h_cheg + min_cheg/60,
    motivo = case_when(
      motivo_o == 8 ~ motivo_d,
      motivo_d == 8 ~ motivo_o,
      (motivo_o != 8 & motivo_d != 8) ~ motivo_d
    ),
    motivo_simples = case_when(
      motivo %in% 1:3 ~ "Trabalho",
      motivo == 4 ~ "Escola",
      TRUE ~ "Outro"
    ),
    dist_km = distancia/1000,
    distancia_viagens = cut(dist_km,
                            breaks = c(0, 1, 3, 6, 83),
                            include.lowest = TRUE)
  )

```

Number of trips / reason

```{r, echo=FALSE}

count(x = OD_2017_v1, motivo, wt = fe_via)

```

```{r, echo=FALSE}

count(x = OD_2017_v1, motivo_simples, wt = fe_via)

```


### Basic temporal analysis

Distributions of trips in the day

```{r, echo=FALSE}

OD_2017_v1 %>%
  filter(!is.na(modoprin)) %>%
  ggplot() +
  geom_density(aes(x=saida,
                   y=after_stat(count),
                   weight=fe_via,
                   color=motivo_simples,
                   linetype="Saída")
               ) +
  geom_density(aes(x=chegada,
                   y=after_stat(count),
                   weight=fe_via,
                   color=motivo_simples,
                   linetype="Chegada")
               ) +
  scale_colour_brewer(name="Motivo",
                      palette="Set1") +
  scale_linetype_manual(name="",
                        values = c("Saída"="solid", "Chegada"="dashed")
                        ) +
  ylab("# de viagens") +
  xlab("Hora") +
  facet_wrap(~modo_simples) +
  theme_bw()

```

```{r, echo=FALSE}

OD_2017_v1 %>%
  filter(!is.na(modoprin)) %>%
  ggplot() +
    geom_density(aes(x=saida,
                   y=after_stat(count),
                   weight=fe_via,
                   color=motivo_simples,
                   linetype="Saída")
               ) +
  geom_density(aes(x=chegada,
                   y=after_stat(count),
                   weight=fe_via,
                   color=motivo_simples,
                   linetype="Chegada")
               ) +
  scale_colour_brewer(name="Motivo",
                      palette="Set1") +
  scale_linetype_manual(name="",
                        values = c("Saída"="solid", "Chegada"="dashed")
                        ) +
  ggtitle("Classificação A/B Streets") +
  ylab("# de viagens") +
  xlab("Hora") +
  facet_wrap(~mode_ab_streets, ncol = 2) +
  theme_bw()

```


### Basic geographic analysis (now with times  )

```{r, include=FALSE}
# # list.files(folder)
# zonas_od = st_read("~/hd/data/project-data/spshared//shapefiles/Zonas_2017_region.shp")
# zonas_od = st_read(file.path(folder, "shapefiles/Zonas_OD2017/Zonas_2017_region.shp")
#                    , options = "ENCODING=LATIN1"
#                    )
# sf::st_write(zonas_od, "zonas_od.gpkg")
# piggyback::pb_upload("zonas_od.gpkg")
# zonas_od = sf::read_sf("zonas_od.gpkg")
# piggyback::pb_download("zonas_od.Rds")
zonas_od = readRDS("zonas_od.Rds")
tmap_mode("plot")                     # mode = view breaks some plots (why? ask Robin)



```

Here I plot the main zones of origin and destination using the total number of trips -- always using the survey's expansion factor. It is interesting to note that there is no centrality pattern when we look at the main Zones of Destination. Maybe we should restrict the analysis to specific types of trips (e.g., commute strictly defined). I checked these results with the official Tables and they are the same in terms of the main zones of origin and destination.

```{r, echo=FALSE}
OD = OD_2017_v1 %>%
  filter(!is.na(modoprin)) %>%
  select(zona_o, zona_d, modoprin, h_saida, fe_via) %>%
  mutate(
    public = ifelse(modoprin %in% 1:6, 1, 0)*fe_via,
    car = ifelse(modoprin %in% 9:10, 1, 0)*fe_via,
    motorcycle = ifelse(modoprin %in% 13:14, 1, 0)*fe_via,
    bike = ifelse(modoprin == 15, 1, 0)*fe_via,
    foot = ifelse(modoprin == 16, 1, 0)*fe_via,
    other = ifelse(modoprin %in% c(7, 8, 11, 12, 17), 1, 0)*fe_via,
    time = case_when(
      h_saida >= 0 & h_saida < 10 ~ "Manhã",
      h_saida >= 10 & h_saida < 15 ~ "Tarde",
      h_saida >=15 & h_saida <= 23 ~ "Noite"
    ),
    across(time, factor, levels=c("Manhã","Tarde","Noite"))
         ) %>%
  select(zona_o, zona_d, time, public, car, motorcycle, bike, foot, other) %>%
  group_by(zona_o, zona_d, time) %>%
  summarise_all(sum) %>%
  ungroup() %>%
  mutate(trips = public + car + motorcycle + bike + foot + other)

# check:
# sum(OD$public) + sum(OD$car) + sum(OD$motorcycle) + sum(OD$bike) + sum(OD$foot) + sum(OD$other)

origin = OD %>%
  group_by(zona_o, time) %>%
  summarise_all(sum) %>%
  select(zona_o, time, trips_o = trips) %>%
  ungroup() %>%
  left_join(zonas_od, by = c("zona_o" = "NumeroZona")) %>%
  st_as_sf(sf_column_name = "geometry") %>%
  mutate(area = units::set_units(st_area(.), km^2),
         trips_o_km2 = as.numeric(trips_o/area)
         )

destination = OD %>%
  group_by(zona_d, time) %>%
  summarise_all(sum) %>%
  select(zona_d, time, trips_d = trips) %>%
  ungroup() %>%
  left_join(zonas_od, by = c("zona_d" = "NumeroZona")) %>%
  st_as_sf(sf_column_name = "geometry") %>%
  mutate(area = units::set_units(st_area(.), km^2),
         trips_d_km2 = as.numeric(trips_d/area)
         )

sp_bounds = read_municipality(code_muni = 3550308) %>%
  st_transform(crs = 4326)
  
```

Plotting the origin zones by time.

```{r, echo=FALSE}

tm_shape(origin) + tm_fill("trips_o_km2", 
                           style = "jenks",
                           title = "Viagens/km²")  +
  tm_facets(by="time", nrow=3) +
  tm_layout(title = "RMSP: Origens",
            inner.margins = c(0, .02, .02, .02),
            legend.outside = TRUE) +
  tm_shape(sp_bounds) +
  tm_borders(col = "red")

```

Plotting the destination zones by time.

```{r, echo=FALSE}

tm_shape(destination) + tm_fill("trips_d_km2", 
                                style = "jenks",
                                title = "Viagens/km²") +
  tm_facets(by="time", nrow=3) +
  tm_layout(title = "RMSP: Destinos",
            inner.margins = c(0, .02, .02, .02),
            legend.outside = TRUE) +
  tm_shape(sp_bounds) +
  tm_borders(col = "red")

```

Looking at the São Paulo City. First, the origins.

```{r, echo=FALSE}

origin %>%
  filter(NomeMunici == "São Paulo") %>%
  tm_shape() +
  tm_fill("trips_o_km2",
          style = "jenks",
          title = "Viagens/km²") +
  tm_facets(by="time", ncol=3) +
  tm_layout(title = "São Paulo - Origens",
            inner.margins = c(0, .02, .02, .02),
            legend.outside = TRUE)

```

Now, the destinations for São Paulo City.

```{r, echo=FALSE}

destination %>%
  filter(NomeMunici == "São Paulo") %>%
  tm_shape() +
  tm_fill("trips_d_km2",
          style = "jenks",
          title = "Viagens/km²") +
  tm_facets(by="time", ncol=3) +
  tm_layout(title = "São Paulo - Destinos",
            inner.margins = c(0, .02, .02, .02),
            legend.outside = TRUE)

```

#### Distances by zone of residence

```{r, echo=FALSE}

pop_zona = OD_2017_v1 %>%
  select(id_pess, zona, fe_pess) %>%
  unique() %>%
  group_by(zona) %>%
  summarise(pop = sum(fe_pess)) %>%
  left_join(zonas_od, by = c("zona"="NumeroZona")) %>%
  st_as_sf(sf_column_name = "geometry") %>%
  mutate(densidade = pop/Area_ha_2)

distancias_zona = OD_2017_v1 %>%
  select(zona, id_pess, distancia, motivo_simples, fe_via, fe_pess) %>%
  group_by(id_pess, zona, motivo_simples) %>%
  mutate(fe_via = replace_na(fe_via, 0),
         distancia = replace_na(distancia, 0)) %>%
  summarise(dist_km = sum((distancia/1000)*fe_via)) %>%
  ungroup() %>%
  group_by(zona, motivo_simples) %>%
  summarise(dist_km = sum(dist_km)) %>%
  ungroup() %>%
  left_join(pop_zona, by = c("zona"="zona")) %>%
  mutate(distancia_media = dist_km/pop) %>%
  st_as_sf(sf_column_name = "geometry")

distancias_commute_zona = OD_2017_v1 %>%
  filter(motivo_simples == "Trabalho") %>%
  select(zona, id_pess, distancia, modo_simples, fe_via, fe_pess) %>%
  group_by(id_pess, zona, modo_simples) %>%
  mutate(fe_via = replace_na(fe_via, 0),
         distancia = replace_na(distancia, 0)) %>%
  summarise(dist_km = sum((distancia/1000)*fe_via)) %>%
  ungroup() %>%
  group_by(zona, modo_simples) %>%
  summarise(dist_km = sum(dist_km)) %>%
  ungroup() %>%
  left_join(pop_zona, by = c("zona"="zona")) %>%
  mutate(distancia_media = dist_km/pop)%>%
  st_as_sf(sf_column_name = "geometry")

```

Looking at the population data here.

```{r, echo=FALSE, eval=FALSE}

populacao = tm_shape(pop_zona) +
  tm_polygons("pop",
              style = "jenks",
              title = "População",
              lwd = .25)

densidade = tm_shape(pop_zona) +
  tm_polygons("densidade",
              style = "jenks",
              title = "População/ha",
              lwd = .25)

tmap_arrange(populacao, densidade, ncol=2)

```

Comparing the population with the Census (2010) using census tracts.

```{r, echo=FALSE, eval=FALSE}

mapa_setores = read_census_tract(code_tract = 35,
                                 year = 2010)

pop_setores = readRDS("./pop_census_tract_2010.Rds")

sf_use_s2(FALSE)       # to avoid that "invalid geometry" error when calculating the centroids 

pop_setores_sp = pop_setores %>%
  filter(id_setor_censitario %in% unique(mapa_setores$code_tract)) %>%
  left_join(mapa_setores, by = c("id_setor_censitario"="code_tract")) %>%
  select(id_setor_censitario, pop = v001, geom) %>%
  st_as_sf() %>%
  st_centroid()     # st_is_longlat(pop_setores_sp) is TRUE

zonas = zonas_od %>%
  select(NumeroZona, NomeZona) %>%
  st_transform(crs = st_crs(pop_setores_sp))

setor_zona = st_join(pop_setores_sp,
                     zonas,
                     st_covered_by)

st_geometry(setor_zona) = NULL

pop_zona_censo = setor_zona %>%
  group_by(NumeroZona) %>%
  summarise(pop_censo = sum(pop)) %>%
  filter(!is.na(NumeroZona))       # dropping the rest of SP State 

pop_zona_od = OD_2017_v1 %>%
  select(id_pess, zona, fe_pess) %>%        # OD is at the trip level, I need one obs/person
  unique() %>%
  group_by(zona) %>%
  summarise(pop_zona = sum(fe_pess))

zonas_pop = zonas %>%
  left_join(pop_zona_od, by = c("NumeroZona"="zona")) %>%
  left_join(pop_zona_censo, by = c("NumeroZona"="NumeroZona")) %>%
  mutate(
    diff_zona_censo = pop_zona - pop_censo,
    diff_perc_od_censo = ((diff_zona_censo)/pop_censo)*100
  )

```

```{r, echo=FALSE, eval=FALSE}

tm_shape(zonas_pop) +
  tm_fill(c("pop_zona", "pop_censo"),
          style = 'jenks',
          title = "População") +
  tm_facets(sync = TRUE, ncol = 2) +
  tm_layout(panel.labels = c("OD2017", "Censo 2010"),
            legend.outside = TRUE)

```

Some people left São Paulo City, but is not significant relative to the total.

```{r, echo=FALSE, eval=FALSE}

diff_map = tm_shape(zonas_pop) +
  tm_fill(c("diff_zona_censo"),
          title = "Diferença (OD - Censo)",
          palette = c("red","white","blue"),
          midpoint = 0) +
  tm_layout(legend.outside = TRUE,
            legend.text.size = 0.3)

perc_change = tm_shape(zonas_pop) +
  tm_fill(c("diff_perc_od_censo"),
          title = "Diferença %",
          breaks = c(seq(-150, 150, by=25)),
          palette = c("red","white","blue"),
          midpoint = 0) +
  tm_layout(legend.outside = TRUE,
            legend.text.size = 0.3)

tmap_arrange(diff_map, perc_change, nrow = 1, ncol = 2)

```


```{r, echo=FALSE, eval=FALSE}

tm_shape(distancias_zona) +
  tm_fill("distancia_media",
          style = "cont",
          title = "Distância média das viagens (km)") +
  tm_facets(by="motivo_simples", ncol=2)

```

```{r, echo=FALSE, eval=FALSE}

tm_shape(distancias_commute_zona) +
  tm_fill("distancia_media",
          style = "jenks",
          title = "Distância média das viagens (km)") +
  tm_facets(by="modo_simples", ncol=2) +
  tm_layout(main.title = "Viagens com motivo trabalho",
            main.title.position = "center")

```


#### Active travels and trip patterns (Desire lines)

(I removed the previous visualizations since they were not saying anything about the data.)

Walking and cycling are highly concentrated in few OD pairs (>90% in 1,000 pairs for walking and >90% in 500 pairs for cycling). Car trips are more dispersed, we have >90% of trips in 10,000 pairs. 

```{r, include=FALSE, eval=FALSE}

bike_pairs = OD %>%
  select(-time) %>%
  group_by(zona_o, zona_d) %>%
  summarise_all(sum) %>%
  ungroup() %>%
  mutate(bike_share = (bike/trips)*100) %>%
  slice_max(bike,
            n = 500,
            with_ties = FALSE)

foot_pairs = OD %>%
  select(-time) %>%
  group_by(zona_o, zona_d) %>%
  summarise_all(sum) %>%
  ungroup() %>%
  mutate(foot_share = (foot/trips)*100) %>%
  slice_max(foot,
            n = 1000,
            with_ties = FALSE)

car_pairs = OD %>%
  select(-time) %>%
  group_by(zona_o, zona_d) %>%
  summarise_all(sum) %>%
  ungroup() %>%
  mutate(car_share = (car/trips)*100) %>%
  slice_max(car,
            n = 10000,
            with_ties = FALSE)

```


Bike trips:

```{r, echo=FALSE, eval=FALSE}

desire_lines_bike = od_jitter(bike_pairs, zonas_od)

tm_shape(desire_lines_bike) +
  tm_lines(palette = "plasma",
           breaks = c(0, 10, 25, 50, 75, 100),
           lwd = "trips",
           scale = 15,
           title.lwd = "# de viagens",
           alpha = 0.8,
           col = "bike_share",
           title = "Viagens de bicicleta (%)") +
  tm_shape(zonas_od) +
  tm_borders(col = "lightgrey", lwd = 0.2) +
  tm_scale_bar()

```

Foot trips:

```{r, echo=FALSE, eval=FALSE}

desire_lines_foot = od_jitter(foot_pairs, zonas_od)

tm_shape(desire_lines_foot) +
  tm_lines(palette = "plasma",
           breaks = c(0, 10, 25, 50, 75, 100),
           lwd = "trips",
           scale = 15,
           title.lwd = "# de viagens",
           alpha = 0.8,
           col = "foot_share",
           title = "Viagens a pé (%)") +
  tm_shape(zonas_od) +
  tm_borders(col = "lightgrey", lwd = 0.2) +
  tm_scale_bar()

```

Car trips, still too many lines (10,000) to account for >92% of total car trips.

```{r, echo=FALSE, eval=FALSE}

desire_lines_car = od_jitter(car_pairs, zonas_od)

tm_shape(desire_lines_car) +
  tm_lines(palette = "plasma",
           breaks = c(0, 10, 25, 50, 75, 100),
           lwd = "trips",
           scale = 15,
           title.lwd = "# de viagens",
           alpha = 0.8,
           col = "car_share",
           title = "Viagens de carro (%)") +
  tm_shape(zonas_od) +
  tm_borders(col = "lightgrey", lwd = 0.2) +
  tm_scale_bar()

```

Using the first 5,000 OD pairs, we have >78% of all car trips.

```{r, echo=FALSE, eval=FALSE}

car_pairs_5 = desire_lines_car %>%
    slice_max(car,
              n = 5000,
              with_ties = FALSE)


desire_lines_car = od2line(car_pairs_5, zonas_od)

tm_shape(desire_lines_car) +
  tm_lines(palette = "plasma",
           breaks = c(0, 10, 25, 50, 75, 100),
           lwd = "trips",
           scale = 15,
           title.lwd = "Number of trips",
           alpha = 0.8,
           col = "car_share",
           title = "Car trips (%)") +
  tm_shape(zonas_od) +
  tm_borders(col = "lightgrey", lwd = 0.2) +
  tm_scale_bar()


```

Using the biggest 1,000 pairs, we have >45% of all car trips.

```{r, echo=FALSE, eval=FALSE}

car_pairs_1 = desire_lines_car %>%
    slice_max(car,
              n = 1000,
              with_ties = FALSE)


desire_lines_car = od2line(car_pairs_1, zonas_od)

tm_shape(desire_lines_car) +
  tm_lines(palette = "plasma",
           breaks = c(0, 10, 25, 50, 75, 100),
           lwd = "trips",
           scale = 15,
           title.lwd = "Number of trips",
           alpha = 0.8,
           col = "car_share",
           title = "Car trips (%)") +
  tm_shape(zonas_od) +
  tm_borders(col = "lightgrey", lwd = 0.2) +
  tm_scale_bar()


```


### Main routes by mode (takes a long time here) -- rewriting this part.

```{r, eval=FALSE, eval=FALSE}

bike_routes = od_jitter(bike_pairs, zonas_od) %>%
  st_transform(crs = 4326) %>%
  mutate(route = st_geometry(line2route(.,
                                        route_fun = route_osrm,
                                        osrm.profile = "bike",
                                        time_delay = 1)
                             )
         )

st_geometry(bike_routes) = bike_routes$route

saveRDS(bike_routes, file = "bike_routes_2021-09-16.Rds")

```

```{r, eval=FALSE, eval=FALSE}

foot_routes = od_jitter(foot_pairs, zonas_od) %>%
  st_transform(crs = 4326) %>%
  mutate(route = st_geometry(line2route(.,
                                        route_fun = route_osrm,
                                        osrm.profile = "foot",
                                        time_delay = 0.5)
                             )
         )

st_geometry(foot_routes) = foot_routes$route

saveRDS(foot_routes, file = "foot_routes_2021-09-16.Rds")

```

```{r, eval=FALSE, eval=FALSE}

car_routes = od_jitter(car_pairs, zonas_od) %>%
  st_transform(crs = 4326) %>%
  mutate(route = st_geometry(line2route(.,
                                        route_fun = route_osrm,
                                        osrm.profile = "car",
                                        time_delay = 0.02)
                             )
         )

st_geometry(car_routes) = car_routes$route

saveRDS(car_routes, file = "car_routes_2021-09-16.Rds")

```



### Plotting the routes

Load the routes I calculated.

```{r, echo=FALSE, eval=FALSE}
bike_routes = readRDS("bike_routes_2021-09-16.Rds")
foot_routes = readRDS("foot_routes_2021-09-16.Rds")
car_routes = readRDS("car_routes_2021-09-16.Rds")           
```

something wrong here, check later...

```{r, echo=FALSE, eval=FALSE}
# remotes::install_cran("archive")
# archive::archive("SAD69-96_SHP_redecicloviaria.rar") # fails
# unrar SAD69-96_SHP_redecicloviaria.rar = bash command
# bike_paths = st_read("SAD69-96_SHP_redecicloviaria.shp")
# st_crs(bike_paths) = 5533     
# bike_paths = st_transform(bike_paths, crs = 4326)
# sf::write_sf(bike_paths, "bike_paths.gpkg")
# piggyback::pb_upload("bike_paths.gpkg")
bike_paths = read_sf("bike_paths.gpkg")
```


```{r, eval=FALSE}
route_net_bike = overline(sl = bike_routes, attrib = "bike")

leaflet(route_net_bike) %>% 
  addProviderTiles(providers$CartoDB.DarkMatter) %>%
  addPolylines(data = bike_paths, color = "green", weight = 1) %>%
  addPolylines(weight = route_net_bike$bike / 1000, color = "white")

```

```{r, eval=FALSE}

route_net_car = overline(sl = car_routes, attrib = "car")

leaflet(route_net_car) %>% 
  addProviderTiles(providers$CartoDB.DarkMatter) %>%
  addPolylines(weight = route_net_car$car / 10000, color = "white")

```

```{r, eval=FALSE}

route_net_foot = overline(sl = foot_routes, attrib = "foot")

leaflet(route_net_foot) %>% 
  addProviderTiles(providers$CartoDB.DarkMatter) %>%
  addPolylines(weight = route_net_foot$foot / 10000, color = "white")

```


### Infrastructure data (São Paulo municipality only)

```{r, echo=FALSE, eval=FALSE}

OD_sp = OD_2017_v1 %>%
  filter(!is.na(modoprin) &
           (muni_d == 36 & muni_o == 36)
         ) %>%
  select(zona_o, zona_d, modoprin, fe_via) %>%
  mutate(
    public = ifelse(modoprin %in% 1:6, 1, 0)*fe_via,
    car = ifelse(modoprin %in% 9:10, 1, 0)*fe_via,
    motorcycle = ifelse(modoprin %in% 13:14, 1, 0)*fe_via,
    bike = ifelse(modoprin == 15, 1, 0)*fe_via,
    foot = ifelse(modoprin == 16, 1, 0)*fe_via,
    other = ifelse(modoprin %in% c(7, 8, 11, 12, 17), 1, 0)*fe_via
    ) %>%
  select(zona_o, zona_d, public, car, motorcycle, bike, foot, other) %>%
  group_by(zona_o, zona_d) %>%
  summarise_all(sum) %>%
  ungroup() %>%
  mutate(trips = public + car + motorcycle + bike + foot + other)

```

```{r, echo=FALSE, eval=FALSE}
# I have two identical geometries in each route dataframe
car_routes$route = NULL
bike_routes$route = NULL
foot_routes$route = NULL

bike_routes$tot_length = (units::set_units(st_length(bike_routes), km))*bike_routes$bike
car_routes$tot_length = (units::set_units(st_length(car_routes), km))*car_routes$car
foot_routes$tot_length = (units::set_units(st_length(foot_routes), km))*foot_routes$foot

```

obs.: very strange behavior here: st_buffer + st_make_valid makes one single polygon very small. without st_buffer, R does not even plot -- "C stack usage too close to ...". These operations seems to be necessary in the other geometries. 
```{r, eval=FALSE}

 od_zones_sp = zonas_od %>%
   filter(NumeroMuni == 36) %>%
#  st_buffer(dist = 0) %>%
   st_transform(4326) %>%
   st_make_valid()

# Calcadas_SP = st_read(file.path(folder,"Poligonos_calcadas/Calcadas_SP.shp"))
# sw = Calcadas_SP %>%  # sidewalks
#   st_buffer(dist = 0) %>%
#   st_transform(4326) %>%
#   st_make_valid()
# saveRDS(sw, "sw.Rds")
# sf::write_sf(sw, "sw.gpkg")
sw = readRDS("sw.Rds")

# streets = st_read(file.path(folder,"Poligonos_vias/Poligonos_vias_total.gpkg")) %>%  # streets
#   st_buffer(dist = 0) %>%
#   st_transform(4326) %>%
#   st_make_valid()
# sf::write_sf(streets, "streets.gpkg")
# saveRDS(streets, "streets.Rds")
# piggyback::pb_upload("streets.gpkg",
#                      repo = "spstreets/OD2017")
# piggyback::pb_upload("streets.Rds",
#                      repo = "spstreets/OD2017")
streets = readRDS("streets.Rds")
 
 area_sw = sw %>%
   st_intersection(od_zones_sp) %>%
   mutate(area_sw = units::set_units(st_area(.), km^2)
          ) %>%
   select(NumeroZona, area_sw) %>%
   st_drop_geometry() %>%
   group_by(NumeroZona) %>%
   summarise_all(sum)
 
 area_streets = streets %>%
   st_intersection(od_zones_sp) %>%
   mutate(area_streets = units::set_units(st_area(.), km^2)
          ) %>%
   select(NumeroZona, area_streets) %>%
   st_drop_geometry() %>%
   group_by(NumeroZona) %>%
   summarise_all(sum)

```

```{r, eval=FALSE}

od_zones = od_zones_sp %>%
  left_join(area_streets, by = c("NumeroZona")) %>%
  left_join(area_sw, by = c("NumeroZona")) %>%
  mutate(area = st_area(.),
         share_sw = as.numeric(area_sw/area)*100,
         share_streets = as.numeric(area_streets/area)*100
         )

```

```{r, eval=FALSE}

street = tm_shape(od_zones) +
  tm_fill(col = "area_streets",
          title = "Streets area (km^2)")

street_share = tm_shape(od_zones) +
  tm_fill(col = "share_streets",
          title = "Street share (%)")

sw = tm_shape(od_zones) +
  tm_fill(col = "area_sw",
          title = "Sidewalk area (km^2)")

sw_share = tm_shape(od_zones) +
  tm_fill(col = "share_sw",
          title = "Sidewalk share (%)")

tmap_arrange(street, sw, street_share, sw_share, ncol = 2)

```

OD pairs that generate the most intensive usage of space (kms traveled in the routes) have more area allocated to the transport mode (both foot and car) 

Very strong pattern here.

```{r, eval=FALSE}

foot_routes %>%
  filter(
    (zona_o %in% unique(OD_sp$zona_o)) & (zona_d %in% unique(OD_sp$zona_d))
    ) %>%
  left_join(area_sw, by = c("zona_o" = "NumeroZona")) %>%
  rename(sw_o = area_sw) %>%
  left_join(area_sw, by = c("zona_d" = "NumeroZona")) %>%
  rename(sw_d = area_sw) %>%
  ggplot(aes(x = log(as.numeric(0.5*sw_o+0.5*sw_d)),
             y = log(as.numeric(tot_length))
             )
         ) +
  geom_point() +
  xlab("log(OD mean sidewalk area)") +
  ylab("log(total km walked OD pair)") +
  theme_bw()

```

Some outliers in this graph, may be worth to check later...

```{r, eval=FALSE}

car_routes %>%
  filter(
    (zona_o %in% unique(OD_sp$zona_o)) & (zona_d %in% unique(OD_sp$zona_d))
    ) %>%
  left_join(area_streets, by = c("zona_o" = "NumeroZona")) %>%
  rename(streets_o = area_streets) %>%
  left_join(area_streets, by = c("zona_d" = "NumeroZona")) %>%
  rename(streets_d = area_streets) %>%
  ggplot(aes(x = log(as.numeric(0.5*streets_o+0.5*streets_d)),
             y = log(as.numeric(tot_length))
             )
         ) +
  geom_point() +
  xlab("log(OD mean streets area)") +
  ylab("log(total km drived in OD pair)") +
  theme_bw()

```