Análisis espacial 2

En este ejercicio se realizará:

Cercanía de venta de frutas y verduras a escuela

Los datos de areas verdes, y otros datos saldrán de Biblioteca digital INEGI

El archivo descargado es un zip que contiene distintos archivos shp que son los siguientes, para hacer más sencilla la lectura (al menos para mí) lo agregué en un objeto que básicamente e suna tabla donde cada renglón contiene la ruta del archivo

archivos = list.files(
  path = "../01_input/702825585341_s/090070001/",
  pattern = "shp",
  full.names = T) %>% 
  as_tibble()

archivos
## # A tibble: 8 × 1
##   value                                                
##   <chr>                                                
## 1 ../01_input/702825585341_s/090070001/090070001A.shp  
## 2 ../01_input/702825585341_s/090070001/090070001E.shp  
## 3 ../01_input/702825585341_s/090070001/090070001L.shp  
## 4 ../01_input/702825585341_s/090070001/090070001M.shp  
## 5 ../01_input/702825585341_s/090070001/090070001SIA.shp
## 6 ../01_input/702825585341_s/090070001/090070001SIL.shp
## 7 ../01_input/702825585341_s/090070001/090070001SIP.shp
## 8 ../01_input/702825585341_s/090070001/090070001T.shp

El archivo no tiene un descriptor de variables y no indica, diccionario o alguna referencia para que sepamos de que va cada archivo, por lo que lo exploraremos, podría usar una función de purr para leerlos todos, pero la ídea es ir viendo uno por uno

p1=sf::st_read(archivos$value[1],quiet = T) %>% sf::st_transform(crs=4326)
p2=sf::st_read(archivos$value[2],quiet = T) %>% sf::st_transform(crs=4326)
p3=sf::st_read(archivos$value[3],quiet = T) %>% sf::st_transform(crs=4326)
p4=sf::st_read(archivos$value[4],quiet = T) %>% sf::st_transform(crs=4326)
p5=sf::st_read(archivos$value[5],quiet = T) %>% sf::st_transform(crs=4326)
p6=sf::st_read(archivos$value[6],quiet = T) %>% sf::st_transform(crs=4326)
p7=sf::st_read(archivos$value[7],quiet = T) %>% sf::st_transform(crs=4326)
p8=sf::st_read(archivos$value[8],quiet = T) %>% sf::st_transform(crs=4326)

son 8 archivos por lo que haré el plot de su primer variable para saber qué es cada uno.

card(
for (i in paste0("p",c(1:8))) {
  get(i) %>%
    select(1) %>% 
    plot()
}
,full_screen = T)

Con las visualizaciones se identifica que son ageb, vialidades, figura del municipio, etc .., los nombres la columna GEOGRÁFICO tienen una codificación, los más comúnes son utf8 o latin 1 por lo que se hará la reclasificación

p5
## Simple feature collection with 991 features and 8 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -99.1397 ymin: 19.28801 xmax: -98.96317 ymax: 19.40069
## Geodetic CRS:  WGS 84
## First 10 features:
##     CODIGO CONDICION  GEOGRAFICO FECHAACT GEOMETRIA INSTITUCIO    NOMBRE
## 1  LA_6161 NO APLICA CAMELL\xd3N  06/2009   \xc1REA      INEGI NO APLICA
## 2  LA_6161 NO APLICA CAMELL\xd3N  06/2009   \xc1REA      INEGI NO APLICA
## 3  LA_6161 NO APLICA CAMELL\xd3N  06/2009   \xc1REA      INEGI NO APLICA
## 4  LA_6161 NO APLICA CAMELL\xd3N  06/2009   \xc1REA      INEGI NO APLICA
## 5  LA_6162 NO APLICA CAMELL\xd3N  06/2009   \xc1REA      INEGI NO APLICA
## 6  LA_6161 NO APLICA CAMELL\xd3N  06/2009   \xc1REA      INEGI NO APLICA
## 7  LA_6162 NO APLICA CAMELL\xd3N  06/2009   \xc1REA      INEGI NO APLICA
## 8  LA_6161 NO APLICA CAMELL\xd3N  06/2009   \xc1REA      INEGI NO APLICA
## 9  LA_6161 NO APLICA CAMELL\xd3N  06/2009   \xc1REA      INEGI NO APLICA
## 10 LA_6161 NO APLICA CAMELL\xd3N  06/2009   \xc1REA      INEGI NO APLICA
##           TIPO                       geometry
## 1  CAMELL\xd3N POLYGON ((-98.99842 19.3592...
## 2  CAMELL\xd3N POLYGON ((-98.99816 19.3598...
## 3  CAMELL\xd3N POLYGON ((-98.99514 19.3585...
## 4  CAMELL\xd3N POLYGON ((-98.9962 19.3591,...
## 5     GLORIETA POLYGON ((-98.99539 19.3585...
## 6  CAMELL\xd3N POLYGON ((-98.99552 19.3583...
## 7     GLORIETA POLYGON ((-98.99821 19.3590...
## 8  CAMELL\xd3N POLYGON ((-98.99649 19.3579...
## 9  CAMELL\xd3N POLYGON ((-98.99541 19.3479...
## 10 CAMELL\xd3N POLYGON ((-98.98652 19.3461...

la función encoding permite transformar distintos tipos de carácteres

Encoding(p5$GEOGRAFICO) ="latin1"
Encoding(p5$NOMBRE) ="latin1"
Encoding(p5$GEOMETRIA) ="latin1"
Encoding(p5$TIPO) ="latin1"
Encoding(p7$GEOGRAFICO) ="latin1"
Encoding(p7$NOMBRE) ="latin1"
Encoding(p7$GEOMETRIA) ="latin1"
Encoding(p7$TIPO) ="latin1"

así se ve ahora la base

p5
## Simple feature collection with 991 features and 8 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -99.1397 ymin: 19.28801 xmax: -98.96317 ymax: 19.40069
## Geodetic CRS:  WGS 84
## First 10 features:
##     CODIGO CONDICION GEOGRAFICO FECHAACT GEOMETRIA INSTITUCIO    NOMBRE
## 1  LA_6161 NO APLICA   CAMELLÓN  06/2009      ÁREA      INEGI NO APLICA
## 2  LA_6161 NO APLICA   CAMELLÓN  06/2009      ÁREA      INEGI NO APLICA
## 3  LA_6161 NO APLICA   CAMELLÓN  06/2009      ÁREA      INEGI NO APLICA
## 4  LA_6161 NO APLICA   CAMELLÓN  06/2009      ÁREA      INEGI NO APLICA
## 5  LA_6162 NO APLICA   CAMELLÓN  06/2009      ÁREA      INEGI NO APLICA
## 6  LA_6161 NO APLICA   CAMELLÓN  06/2009      ÁREA      INEGI NO APLICA
## 7  LA_6162 NO APLICA   CAMELLÓN  06/2009      ÁREA      INEGI NO APLICA
## 8  LA_6161 NO APLICA   CAMELLÓN  06/2009      ÁREA      INEGI NO APLICA
## 9  LA_6161 NO APLICA   CAMELLÓN  06/2009      ÁREA      INEGI NO APLICA
## 10 LA_6161 NO APLICA   CAMELLÓN  06/2009      ÁREA      INEGI NO APLICA
##        TIPO                       geometry
## 1  CAMELLÓN POLYGON ((-98.99842 19.3592...
## 2  CAMELLÓN POLYGON ((-98.99816 19.3598...
## 3  CAMELLÓN POLYGON ((-98.99514 19.3585...
## 4  CAMELLÓN POLYGON ((-98.9962 19.3591,...
## 5  GLORIETA POLYGON ((-98.99539 19.3585...
## 6  CAMELLÓN POLYGON ((-98.99552 19.3583...
## 7  GLORIETA POLYGON ((-98.99821 19.3590...
## 8  CAMELLÓN POLYGON ((-98.99649 19.3579...
## 9  CAMELLÓN POLYGON ((-98.99541 19.3479...
## 10 CAMELLÓN POLYGON ((-98.98652 19.3461...

los añadimos a un objeto llamado: área verde y otro llamado escuelas, aunque podríamos reescribir p5 y p7

area_verde = p5 %>% 
  filter(GEOGRAFICO=="ÁREA VERDE")

escuelas = p7 %>%
  filter(GEOGRAFICO=="ESCUELA")

DENUE

Directorio nacional de unidades económicas, DENEU acceder a los datos d ela denue puede ser de distintas formas, una forma sencilla en su página seleccionar el área geográfica y visualizar los resultados, las clasificaciones de negocios se basan en el SCIAN por lo que es posible obtener sectores,subsectores o negocios muy específicos, para este ejercicio buscaremos los que estén relacionados al comercio al por menor de frutas y verduras frescas

posteriormente debemos darle en descargar, sin son muchos negocios nos mandará a la descarga masiva, pero si son pocos desde esa interfaz.

denue = sf::st_read("../01_input/denue_iztapalapa/denue_iztapalapa.shp")
## Reading layer `denue_iztapalapa' from data source 
##   `C:\Users\DELL\Desktop\Dip_UNAM_2025\01_input\denue_iztapalapa\denue_iztapalapa.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 87905 features and 42 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 479522.1 ymin: 2100827 xmax: 590030.7 ymax: 2148594
## Projected CRS: WGS 84 / UTM zone 14N
denue = denue %>% 
  filter(nombre_act=="Comercio al por menor de frutas y verduras frescas") %>% 
  sf::st_transform(crs=4326)
library(leaflet)
card(
     class = "p-0",
    card_header(
    class = "bg-dark",
    "Escuelas"
  ),
leaflet() %>% 
  addProviderTiles(providers$OpenStreetMap) %>% 
  addPolygons(data=p3,
              weight = 2,
              fill = NA,
              color="black",
              opacity = 1) %>% 
  addCircles(data=denue)
,full_screen = T)
Escuelas

Visualización con otras capas

En el siguiente mapa se muestran:

vialidades = p2

card(
   class = "p-0",
    card_header(
    class = "bg-dark",
    "Escuelas, áreas verdes y avenidas"
  ),
leaflet() %>% 
  addProviderTiles(providers$CartoDB,group = "carto") %>%
  addProviderTiles(providers$OpenStreetMap,group = "OSM") %>% 
  addPolygons(data=area_verde,
              weight = 1,
              color = "green",
              fillColor = "green",
              fillOpacity = 1) %>% 
  addCircles(data=escuelas,opacity = 1) %>% 
  addCircles(data=denue,opacity = 1,color = "red") %>% 
  addPolylines(data=vialidades %>%
                 filter(TIPOVIAL=="AVENIDA"),
               weight = 2,
               color = "black",
               opacity = 1) %>% 
  addPolygons(data=p3,
              fill = NA,
              color = "black",
              opacity = 1,
              weight = 2
              ) %>% 
    addLayersControl(baseGroups = c("carto",
                                    "OSM"),
                     options = layersControlOptions(collapsed=F)
                     )
,full_screen = T)
Escuelas, áreas verdes y avenidas

Reproyectar y cálculo de distancias

se reproyectará a 32614 que corresponde a la zona 14 y se hará el cálculo para a cada escuela asignarle la unidad de negocio más cercano, calcularemos las distancias así como sus rutas.

Reproyectar a 32614

denue=denue %>% sf::st_transform(crs=32614)
area_verde = area_verde %>% sf::st_transform(32614)

Calculo de distancias

con la función sf_nearest feature obtenemos la posición más cercana a través de la geometría, posteriormente se calcula la distancia entre esos dos puntos con la función de distance, el mismo proceso se hace para las áreas verdes, al tener las geometrías en 32614 mostrará la distancia en metros.

distancias = escuelas %>% 
  sf::st_transform(crs=32614) %>% 
  mutate(id_escuela = 1:n()) %>% 
  select(CODIGO, id_escuela) %>% 
  mutate(nearest_id = sf::st_nearest_feature(geometry, denue)) %>%
  mutate(distancia = sf::st_distance(geometry, denue$geometry[nearest_id], by_element = TRUE))

distancias_verdes = escuelas %>% 
  sf::st_transform(crs=32614) %>% 
  mutate(id_escuela = 1:n()) %>% 
  select(CODIGO, id_escuela) %>% 
  mutate(nearest_id = sf::st_nearest_feature(geometry, area_verde)) %>%
  mutate(distancia = sf::st_distance(geometry, area_verde$geometry[nearest_id], by_element = TRUE))

Visualización de distancias

El histograma nos ayuda a visualizar la distribución de la variables

grafico  = distancias %>% 
  ggplot(aes(distancia %>% as.numeric()))+
  geom_histogram(aes(fill=..x..))+
  theme_bw()+
  viridis::scale_fill_viridis(option = "D",labels=scales::comma)+
  labs(title = "Distribución de distancia",
       subtitle = "Escuelas y unidades de verduras más cercanas",
       x="Frecuencia de distancia en metros"
       )
  

grafico_verde = distancias_verdes %>% 
  ggplot(aes(distancia %>% as.numeric()))+
  geom_histogram(aes(fill=..x..))+
  theme_bw()+
  viridis::scale_fill_viridis(option = "D",labels=scales::comma)+
  labs(title = "Distribución de distancia",
       subtitle = "Escuelas y area verde más cercana",
       x="Frecuencia de distancia en metros"
  )


grafico
## Warning: The dot-dot notation (`..x..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(x)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

grafico_verde
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

En escuelas a la unidad de frutas y verduras la distancia más grande es de 600 metros, mientras que una escuela a un área verde es de hasta 2,000 metros. Hay que tomar en cuenta que estas distancias son en línea recta sin considerar las vialidades o barreras físicas que podrían limitar la movilidad.

EN promedio las escuelas tienen una distancia a la unidade de frutas más cercana a 119.52 metros, mientras la mayor distancia es de 627 y la menor de 16.31 metris para las áreas verdes la distancia promedio es de 470 metros y la mayor de 2,034 metros y la menor de 11 metros.

summary(distancias$distancia)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   16.31   71.82  119.52  143.55  183.64  627.07
summary(distancias_verdes$distancia)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   11.57  187.91  470.77  553.33  815.49 2035.75

Distancia máxima

distancias %>% 
  filter(distancia==max(distancia))
## Simple feature collection with 1 feature and 4 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 493661.2 ymin: 2141462 xmax: 493661.2 ymax: 2141462
## Projected CRS: WGS 84 / UTM zone 14N
##   CODIGO id_escuela                 geometry nearest_id    distancia
## 1 DP_322         50 POINT (493661.2 2141462)       1651 627.0726 [m]
escuela_cercana = denue[1651,]

card(
    class = "p-0",
    card_header(
    class = "bg-dark",
    "Escuela con la Distancia máxima en unidad de frutas y verduras"
  ),
leaflet() %>% 
  addProviderTiles(providers$OpenStreetMap) %>% 
  addCircles(data=distancias %>% 
               filter(distancia==max(distancia)) %>% sf::st_transform(4326),opacity = 1) %>% 
  addCircles(data=escuela_cercana %>% sf::st_transform(crs=4326),color = "red",opacity = 1)
,full_screen = T
)
Escuela con la Distancia máxima en unidad de frutas y verduras

Linea de deseo y calculo de ruta

library(sf)
## Linking to GEOS 3.12.2, GDAL 3.9.3, PROJ 9.4.1; sf_use_s2() is TRUE
origen= distancias %>% 
  filter(distancia==max(distancia)) %>% 
  sf::st_transform(4326) %>% 
  select(1)

destino= escuela_cercana %>% 
  sf::st_transform(crs=4326)


linea <- st_linestring(rbind(
  st_coordinates(origen),
  st_coordinates(destino)
))

ruta_caminando = osrm::osrmRoute(
  src = origen,
  dst = destino,
  osrm.profile = "foot"
  )

card(
  class = "p-0",
    card_header(
    class = "bg-dark",
    "Escuela con la Distancia máxima en unidad de frutas y verduras"
  ),
leaflet() %>% 
  addProviderTiles(providers$OpenStreetMap) %>% 
  addCircles(data=distancias %>% 
               filter(distancia==max(distancia)) %>% sf::st_transform(4326),opacity = 1) %>% 
  addCircles(data=escuela_cercana %>% sf::st_transform(crs=4326),color = "red",opacity = 1) %>% 
  addPolylines(data=ruta_caminando,
               color = "green",
               opacity = 1) %>% 
  addPolylines(data=linea,
               color = "black"
               )
,full_screen = T)
Escuela con la Distancia máxima en unidad de frutas y verduras

la distancia calculada con OSRM indica que se hace 11 minutos caminando y son 840 metros

ruta_caminando
## Simple feature collection with 1 feature and 4 fields
## Geometry type: LINESTRING
## Dimension:     XY
## Bounding box:  xmin: -99.06091 ymin: 19.36169 xmax: -99.0591 ymax: 19.36726
## Geodetic CRS:  WGS 84
##        src  dst duration distance                       geometry
## 1_1651   1 1651 11.20833   0.8408 LINESTRING (-99.06083 19.36...

Como ejercicio se calcularan únicamente las líneas de deseo a la unidad más cercana, el ruteo se puede hacer de forma vectorial ya sea con OSRM o con mapbox

contenedor=data.frame()
for (i in distancias$id_escuela) {
  
  origen = distancias %>% 
    filter(id_escuela==i) %>% 
    sf::st_transform(4326) %>% 
    select(1,nearest_id)
  
  destino= denue[origen$nearest_id,] %>% 
  sf::st_transform(crs=4326)
  
  linea = st_linestring(
  rbind(
    st_coordinates(origen),
    st_coordinates(destino)
    )
  ) %>% 
  st_sfc(.) %>% 
  sf::st_sf(id=i)
  
  contenedor= rbind(contenedor,linea)
  
}

Visualización de líneas

library(bslib)

card(
  class = "p-0",
    card_header(
    class = "bg-dark",
    "Unidades de Frutas y verduras asociada a la escuela más cercana"
  ),
  card_body(
    leaflet() %>% 
      addProviderTiles(providers$CartoDB) %>% 
      addPolygons(
        data=p3,
        weight = 2,
        opacity = 1,
        fill = NA,
        color = "black"
        ) %>% 
      addPolylines(data=contenedor,color = "#4daf4a",opacity = 1,weight = 2) %>% 
      addCircles(data=distancias %>% sf::st_transform(crs=4326),opacity = 1,label = "Escuela") %>%
      addCircles(data=denue[distancias$nearest_id,] %>% sf::st_transform(crs=4326),color = "red",opacity = 1,label = ~nom_estab)
    ),full_screen = TRUE
)
Unidades de Frutas y verduras asociada a la escuela más cercana

Versión con rutas

Visualización de rutas y comparación con lineas de deseo.

Las rutas se obtienen de un proceso llamando a la api de MapBox, pero puede ser de distintas formas, con google maps, con osrm, con el análisis de nodos etc.

card(
  class = "p-0",
    card_header(
    class = "bg-dark",
    "Linea de deseo y ruta estimada"
  ),
  card_body(
leaflet() %>% 
  addProviderTiles(providers$OpenStreetMap) %>% 
  addPolylines(data=contenedor %>% 
  filter(id%in%rutas_mb$viaje),color = "black",
  opacity = 1,weight = 2) %>% 
  addPolylines(data=rutas_mb,color="red",opacity = 1,weight = 2)
))
Linea de deseo y ruta estimada

Distribución de tiempo de traslados (muestra de 300 OD)

rutas_mb %>% 
  ggplot(aes(duration))+
  geom_histogram(aes(fill=..x..))+
  theme_bw()
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.