CDMX_bbox <- getbb("Ciudad de Mexico")
CDMX_bbox_poly <- getbb("Ciudad de Mexico", format_out = "sf_polygon")

leaflet(CDMX_bbox_poly) %>%
  addTiles() %>%
  addPolygons()

#Descargamos shapefile de calles, avenidas y autopistas de CDMX

CDMX_calles <- opq(CDMX_bbox) %>%
  add_osm_feature(key = "highway") %>%
  osmdata_sf()
CDMX_calles <- CDMX_calles$osm_lines
CDMX_calles <- st_intersection(CDMX_calles, filter(CDMX_Alcaldias, nomgeo == "Cuauhtémoc" | nomgeo == "Benito Juárez"))
## although coordinates are longitude/latitude, st_intersection assumes that they are planar
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries

#Nos centramos en mapear solo las calles de las alcaldías CUAUHTEMOC y BENITO JUAREZ, sino para toda la ciudad queda una mancha negra imposible de visualizar. No pameamos según ninguno de los atributos (ejemplo: “velocidad máxima”, “carriles”, etc. ya la la mayoría aparecen con todos los campos en NA). En cambio, construimos un layer superpuesto que demarca todas las vías que contienen “Avenida” en su nombre.

ggplot() +
  geom_sf(data= CDMX_calles, color = "gray40") +
  geom_sf(data = CDMX_Alcaldias %>%
            filter(nomgeo == "Cuauhtémoc" | nomgeo == "Benito Juárez"), color = "black", fill = NA) +
  geom_sf(data= filter(CDMX_calles, str_detect(name, "Avenida"),), colour = "red") +
  labs(title= "CDMX: Aveniadas",
       subtitle= "Avenidas de alcaldías Cuauhtemoc y Benito Juarez",
       caption= "Fuente: OpenStreetMap")

#Descargamos shapefile con estaciones de servicio y de concesionarios automotrices

CDMX_fuel_query <- opq(CDMX_bbox) %>%
  add_osm_feature(key = "amenity", value = "fuel")

CDMX_fuel <- osmdata_sf(CDMX_fuel_query)

CDMX_fuel <- CDMX_fuel$osm_points
CDMX_fuel <- st_intersection(CDMX_fuel, CDMX_bbox_poly)
## although coordinates are longitude/latitude, st_intersection assumes that they are planar
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
CDMX_auto <- opq(CDMX_bbox) %>%
  add_osm_feature(key = "shop", value = "car") %>%
  osmdata_sf()
CDMX_auto <- CDMX_auto$osm_points
CDMX_auto <- st_intersection(CDMX_auto, CDMX_bbox_poly)
## although coordinates are longitude/latitude, st_intersection assumes that they are planar
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries

#Mapeamos las gasolineras y concesionarios automotrices

ggplot()+
  geom_sf(data= CDMX_Alcaldias, color= "grey", fill= NA) +
  geom_sf(data= CDMX_fuel, color= "red", size= 0.3) +
  geom_sf(data= CDMX_auto, color= "blue", size = 0.3) +
  labs(title = "CDMX: Gasolineras (rojo) 
y Concesionarios Automotrices (azul)",
       caption= "Fuente: OpenStreetMap")

#Muy bonito, aunque la CDMX es tan tan grande que esta división por ALCALDÍAS no nos permite ver con mayor granularidad la distribución espacial de las estaciones de servicios. Hagamos zoom in a dos ALCALDÍAS céntricas: la Benito Juarez y la Cuauhtémoc. Luego, dividiremos estas dos alcaldías en sus respectivas COLONIAS (suerte de demarcaciones barriales, más pequeñas). No hemos fraccionado la CDMX entera en COLONIAS, ya que muchas de ellas no coinciden exactamente con las demarcaciones por ALCALDÍAS, o están superpuestas en 2 alcaldías distintas o incluso con el Estado de México (provincia). Más aún, la CDMX tiene más de mil COLONIAS, lo cual la vuelve poco “visualizable”. Así mismo, el censo nacional divide la CDMX en AGEBs, que no tinenen relación con las colonias. Es decir, tenemos distintos niveles alternativos de fraccionamiento del territorio según data que nos interese analizar.

CDMX_colonias <- read_sf("/Users/martinsingla/Documents/UTDT - MEU/Ciencia de Datos para Ciudades II/INFO SELECCIONADA/CDMX_Colonias/coloniascdmx.shp")
CDMX_colonias_select <- st_intersection(CDMX_colonias, 
                                        filter(CDMX_Alcaldias, nomgeo == "Cuauhtémoc" | nomgeo == "Benito Juárez"))
## although coordinates are longitude/latitude, st_intersection assumes that they are planar
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
CDMX_fuel_colonias <- st_intersection(CDMX_fuel, CDMX_colonias_select)
## although coordinates are longitude/latitude, st_intersection assumes that they are planar
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
ggplot() +
  geom_sf(data = CDMX_colonias_select) +
  geom_sf(data = CDMX_fuel_colonias, color = "red", size = 1, alpha = 0.5) +
  theme_minimal() +
  labs(title= "Ubicación de gasolineras",
       subtitle= "CDMX, alcaldías Cuauhtémoc y Benito Juárez",
       caption= "Fuente: OpenStreetMap")

CDMX_colonias_select$fuelstation_count <- lengths(st_intersects(CDMX_colonias_select, CDMX_fuel_colonias))
## although coordinates are longitude/latitude, st_intersects assumes that they are planar
CDMX_colonias_select$area_km2 <- as.numeric(st_area(CDMX_colonias_select))

ggplot() +
  geom_sf(data = CDMX_colonias_select, aes(fill= fuelstation_count/ area_km2)) +
  labs(title= "CDMX: Densidad de gasolineras",
       subtitle= "Dens. de gasolineras por colonia 
       (Alcaldías Cuauhtémoc y Benito Juárez)",
       caption= "Fuente: OpenStreetMap") +
  scale_fill_viridis_c()

#Ahora vemamos lo mismo para los concesionarios automotrices.

CDMX_auto_colonias <- st_intersection(CDMX_auto, CDMX_colonias_select)
## although coordinates are longitude/latitude, st_intersection assumes that they are planar
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
ggplot() +
  geom_sf(data= CDMX_colonias_select) +
  geom_sf(data= CDMX_auto_colonias, color= "blue", size= 1) +
  labs(title = "Ubicación de Concesionarios automotrices",
       subtitle = "CDMX, Alcaldías Cuauhtémoc y Benito Juárez",
       caption = "Fuente: OpenStreetMap")

CDMX_colonias_select$auto_count <- lengths(st_intersects(CDMX_colonias_select, CDMX_auto_colonias))
## although coordinates are longitude/latitude, st_intersects assumes that they are planar
ggplot() +
  geom_sf(data= CDMX_colonias_select, aes(fill = auto_count / area_km2)) +
  scale_fill_viridis_c() +
  labs(title = "Densidad de Concesionarios automotrices",
       subtitle = "CDMX, Alcaldías Cuauhtémoc y Benito Juárez",
       caption = "Fuente: OpenStreetMap")

#Observaciones: hemos identificado al menos tres clusters de concesionarios automotrices en estas 2 alcaldías (la Cuauhtemoc y la Benito Juárez) de CDMX.El resto de los concesionarios aparecen desperdigados por distintas colonias. Contrario a eso, las gasolineras aparecen en clusters desperdigadas por muchas colonias de estas dos alcaldías, y especialmente concentradas en la zona centro-oeste. El cluster sur de concesionarios automotrices coincide con uno de gasolineras también. Se atribuye esta zona a un centro de servicios vehiculares.

#Crítica 1: sé por haber vivido en CDMX, que hay puntos faltantes, hay varios concesionarios automotrices que no han sido mapeados en OpenStreetMap, por lo que no figuran en nuestros datos. Un estudio de geo-marketing sobre concesionarios automotrices debería tener en cuenta fuentes más completas. Lo mismo puede llegar a estar aplicando para gasolineras.

#Crítica 2: tal vez podríamos pensar en otra metodología para mapear clusters de locales, ya que por la demarcación por colonias tal vez no se está visualizando correctamente esto. Claramente las colonias contiguas a aquellas donde se concentran concesionarios, también tienen 1 o 2 agencias, pero no están siendo resaltadas debidamente.