Video con la explicación del código https://youtu.be/5gJYygomOIo?si=Ru7GOMKJv3HGn0i2

Datos del Directorio Estadístico Nacional de Unidades Económicas (DENUE) en la plataforma del INEGI

https://www.inegi.org.mx/app/descarga/

Para el estado de Coahuila de Zaragoza

https://www.inegi.org.mx/contenidos/masiva/denue/denue_05_shp.zip

Contexto geográfico

ruta<-"C:/Users/cguer/Documents/Claudia/Midropbox/Investigacion y escritos/karamanis/data/conjunto_de_datos/"
infile <- "denue_inegi_05_.shp"
nomarchi<-as.filename(paste0(ruta,infile))


a<-"SELECT * FROM denue_inegi_05 WHERE cve_mun ='035'AND codigo_act in ('622111','622112','622211','622212','622311','622312')"


a<-"SELECT * FROM denue_inegi_05 WHERE cve_mun ='035'AND codigo_act in ('622111','622112','622211','622212')"
datadenuet<- read_sf(nomarchi,options = "ENCODING=latin1",query=a)


#datadenuet<- read_sf(nomarchi,options = "ENCODING=latin1")



datadenuet= st_transform(datadenuet, crs = 4326)
datadenuet_geom <- st_geometry(datadenuet)

Incorporación de la zona municipal de Torreón en el polígono urbano Uso de consultas para la lectura filtrada del Shapefile del contexto estatal filtrado municipal

Contexto geográfico

https://inegi.org.mx/contenidos/productos/prod_serv/contenidos/espanol/bvinegi/productos/geografia/marcogeo/794551132173_s.zip

ruta<-"C:/Users/cguer/Documents/Claudia/Midropbox/Investigacion y escritos/karamanis/data/conjunto_de_datos/"
infile <- "05mun.shp"
nomarchi<-as.filename(paste0(ruta,infile))
nomfil<-"05mun"
a<-"SELECT * FROM \"05mun\" WHERE CVE_ENT = '05' AND CVE_MUN='035'"


datageot<- read_sf(nomarchi,options = "ENCODING=latin1",query=a)
datageot= st_transform(datageot, crs = 4326)
datageo_geom <- st_geometry(datageot)
polygon = st_cast(datageo_geom, "POLYGON") 


tor<-polygon[2]

Incorporación de contexto vial como referencia local de ubicación uso de librería OSMDATA Elección territorial del polígono de Torreón urbano (polygon[2])

Delimitación en el área urbana del municipio de Torreón

Código de actividad SCIAN https://www.inegi.org.mx/scian/#

city = "Torreón Coahuila"

viasosm <- opq(city) %>% 
  add_osm_feature(key = "highway") %>% 
  osmdata_sf()

viasosm_sf <- viasosm$osm_lines 
st_crs(viasosm_sf) <-4326

recorteviasosm_sf <-viasosm_sf[tor, ] 


recorteviasosm_sf <-viasosm_sf|> st_intersection(tor)
recorteviasosm_sf <-subset(recorteviasosm_sf,recorteviasosm_sf$highway %in% c("primary","secondary","tertiary"))

recortedenue_sf <-datadenuet|> st_intersection(tor)
recortedenue_sf <-subset(recortedenue_sf,recortedenue_sf$codigo_act %in% c("622111","622112","622211","622212","622311","622312"))
#recortedenue_sf <-subset(recortedenue_sf,recortedenue_sf$codigo_act %in% c("541120"))

recortedenue_sf <-subset(recortedenue_sf,recortedenue_sf$per_ocu %in% c("101 a 250 personas","251 y más personas","51 a 100 personas"))
recortedenue_sf$id2<-paste0(recortedenue_sf$latitud,recortedenue_sf$longitud)
recortedenue_sf$id3<- as.integer(row.names(recortedenue_sf)) 

# Plotting the rough maps
g <- ggplot() +
  geom_sf(data = tor, fill = NA) +
  geom_sf(data = recorteviasosm_sf, fill = NA, linewidth = 0.1) +
  geom_sf(
    data = recortedenue_sf,
   # alpha = 0.3,
    size = 1,
    color="red"
  ) +
  labs(
    title = "Puntos DENUE"
  ) +
  theme_minimal(
    base_family = "body_font",
    base_line_size = unit(0.1, "mm")
  ) +
  theme(
    plot.title = element_text(size = 30)
  )
g

Preparación con la librería de OSRM que proporciona la información de las rutas más cortas Con el paquete osrm que dirige las consultas al servidor OSRM

Generación de la simulación de puntos dentro del polígono municipal de Torreón para la deerminación de la ruta más corta al hospital más cercano

torpuntos <- st_cast(tor, "POINT")

#Get unique coords
coords = as.data.frame(unique(st_coordinates(torpuntos)))
lonmin<-min(coords$X)
lonmax<-max(coords$X)
latmin<-min(coords$Y)
latmax<-max(coords$Y)
#Asses minimimun and unique

#Generación de 1000 puntos distribuidos dentro de esos márgenes
# Generación de dos corridas de 1000 con seed 123 y seed 321  para crear dftt y no perder las coordenadas de fuentes y destinos
n <- 1000

set.seed(321)
puntosaleaotrios <- tibble(id=1:n,
                    lon = runif(n, min = lonmin, max = lonmax),
                    lat = runif(n, min = latmin, max = latmax))

#class(puntosaleaotrios)

puntosaleatorios_sf<-puntosaleaotrios %>% st_as_sf(coords = c("lon", "lat"), crs = st_crs(4326))
puntosaleatorios_sf <-puntosaleatorios_sf|> st_intersection(tor)


#Acotar los puntos a las zonas habitadas por las AGEB


ruta<-"C:/Users/cguer/Documents/Claudia/Midropbox/Investigacion y escritos/karamanis/data/conjunto_de_datos/"
infile <- "05a.shp"
nomarchi<-as.filename(paste0(ruta,infile))
nomfil<-"05a"
a<-"SELECT * FROM \"05a\" WHERE CVE_ENT = '05' AND CVE_MUN='035'"


datageotab<- read_sf(nomarchi,options = "ENCODING=latin1",query=a)
datageotab= st_transform(datageotab, crs = 4326)

puntosaleatorios_sf <-puntosaleatorios_sf|> st_intersection(datageotab)
puntosaleatorios_sf$id<- as.integer(row.names(puntosaleatorios_sf)) 


# Plotting the rough maps
g <- ggplot() +
  geom_sf(data = tor, fill = NA) +
  geom_sf(data = recorteviasosm_sf, fill = NA, linewidth = 0.1) +
  geom_sf(
    data = recortedenue_sf,
   # alpha = 0.3,
    size = 1,
    color="red"
  ) +
  
  geom_sf(
    data = puntosaleatorios_sf,
   # alpha = 0.3,
    size = 1,
    color="navyblue"
  ) +
  labs(
    title = "Hospitales y puntos aleatorios"
  ) +
  theme_minimal(
    base_family = "body_font",
    base_line_size = unit(0.1, "mm")
  ) +
  theme(
    plot.title = element_text(size = 30)
  )
g

La identificación de la ruta a través del servido OSRM Se utiliza la función osrmTable(origen,destino) El origen será simulado por los puntos aleatorios y el destino serán las unidades hospitalarias identificadas en el DENUE. Se alcula el tiempo de traslado

#Tiempo de cálculo
tinicial <- Sys.time()

tabladedistancias<- osrmTable(src = puntosaleatorios_sf, dst = recortedenue_sf)

#diferencia del tiempo inical y de cuando concluye

Sys.time() - tinicial

Time difference of 0.884567 secs

## Time difference of 0.7395082 secs
### esto se genero y guardó en dos corridas 
#tabladedistancias %>% summary()
#tabladedistancias$sources
#tabladedistancias$durations
#tabladedistancias$destinations

tabdes<-bind_cols(tabladedistancias$destinations[1],tabladedistancias$destinations[2])
tabdes$id2<-paste0(tabdes$lat,tabdes$lon)
tabdur<-as.tibble(tabladedistancias$durations)


nh<-as.data.frame(tabdur[1])
nh$dest<-1
colnames(nh)[1]<-"tiempo"
nh$source<- as.integer(row.names(nh)) 
nh$minumin<-min(nh$tiempo)
#nht<-nh

for (i in 2:length(tabdur)) {

nh2<-as.data.frame(tabdur[i])
nh2$dest<-i
colnames(nh2)[1]<-"tiempo"
nh2$source<- as.integer(row.names(nh2)) 
nh2$minumin<-min(nh2$tiempo)

nh<-bind_rows(nh,nh2)
}
nht<-nh

## df3 y #df4 oara juntar las dos corridas y guardadas en dftt
df3 <- left_join(nht, puntosaleatorios_sf, by=c('source'='id'))
df4 <- left_join(nh, puntosaleatorios_sf, by=c('source'='id'))
dft<-bind_rows(df3,df4)
dft<-dplyr::select(dft, tiempo, dest, geometry)
dft<-dft %>% distinct(across(everything()))
dftt <- left_join(dft, recortedenue_sf, by=c('dest'='id3'))
#st_write(dftt, nomfile)

De las distancias obtenidas de cada punto a cada hospital, se obtienen los tiempos más cortos

distanciasminimas <- bind_cols(tabladedistancias$sources, mintime = apply(tabladedistancias$durations, 1, min))%>% 
  as_tibble() %>% 
  mutate(mintime_str = as.character(mintime)) %>% 
  distinct()


#summary(distanciasminimas$mintime)

distanciasminimas_sf<-distanciasminimas %>% st_as_sf(coords = c("lon", "lat"), crs = st_crs(4326))
#install.packages("devtools")
#devtools::install_github("yutannihilation/ggsflabel")
library(ggsflabel)
g <- ggplot() +
  geom_sf(data = tor, fill = NA) +
  geom_sf(data = recorteviasosm_sf, fill = NA, linewidth = 0.1) +
  geom_sf(
    data = recortedenue_sf,
   # alpha = 0.3,
    size = 1,
    color="red"
  ) +
  
  geom_sf(
    data = distanciasminimas_sf,
   # alpha = 0.3,
    size = 1,
    color="navyblue"
  ) +
    #geom_sf_text_repel(data=distanciasminimas_sf,aes(label=mintime_str), size=2, color="navyblue")+

  labs(
    title = "Hospitales y puntos aleatorios"
  ) +
  theme_minimal(
    base_family = "body_font",
    base_line_size = unit(0.1, "mm")
  ) +
  theme(
    plot.title = element_text(size = 30)
  )
g

Simulación de 10000 puntos de 100 en 100

nruns <- 100
n <- 100

set.seed(13)

for (r in 1:nruns) {
  
 puntosaleaotrios <- tibble(id=1:n,
                    lon = runif(n, min = lonmin, max = lonmax),
                    lat = runif(n, min = latmin, max = latmax))
puntosaleatorios_sf<-puntosaleaotrios %>% st_as_sf(coords = c("lon", "lat"), crs = st_crs(4326))
puntosaleatorios_sf <-puntosaleatorios_sf|> st_intersection(tor)
puntosaleatorios_sf <-puntosaleatorios_sf|> st_intersection(datageotab)
  
  tabladedistancias<- osrmTable(src = puntosaleatorios_sf, dst = recortedenue_sf) 
  distanciasminimasi <- bind_cols(tabladedistancias$sources, mintime = apply(tabladedistancias$durations, 1, min))%>% 
  as_tibble() %>% 
  mutate(mintime_str = as.character(mintime)) %>% 
  distinct()

  distanciasminimas <- bind_rows(distanciasminimas, distanciasminimasi )
  #tabladedistancias$sources
  #tabladedistancias$durations
  #tabladedistancias$destinations

  tabdes<-bind_cols(tabladedistancias$destinations[1],tabladedistancias$destinations[2])
  tabdes$id2<-paste0(tabdes$lat,tabdes$lon)
  tabdur<-as.tibble(tabladedistancias$durations)

  length(tabdur)
  nh<-as.data.frame(tabdur[1])
  nh$dest<-1
  colnames(nh)[1]<-"tiempo"
  nh$source<- as.integer(row.names(nh)) 
  nh$minumin<-min(nh$tiempo)


for (i in 2:length(tabdur)) {

nh2<-as.data.frame(tabdur[i])
nh2$dest<-i
colnames(nh2)[1]<-"tiempo"
nh2$source<- as.integer(row.names(nh2)) 
nh2$minumin<-min(nh2$tiempo)
nh<-bind_rows(nh,nh2)
  
}
  nht<-bind_rows(nht,nh)
}
nht %>% distinct(across(everything()))
distanciasminimas_sf<-distanciasminimas %>% st_as_sf(coords = c("lon", "lat"), crs = st_crs(4326))


g <- ggplot() +
  geom_sf(data = tor, fill = NA) +
  geom_sf(data = recorteviasosm_sf, fill = NA, linewidth = 0.1) +
  geom_sf(
    data = recortedenue_sf,
   # alpha = 0.3,
    size = 1,
    color="red"
  ) +
  
  geom_sf(
    data = distanciasminimas_sf,
   # alpha = 0.3,
    size = 1,
    color="navyblue"
  ) +
    #geom_sf_text_repel(data=distanciasminimas_sf,aes(label=mintime_str), size=2, color="navyblue")+

  labs(
    title = "Hospitales y puntos aleatorios"
  ) +
  theme_minimal(
    base_family = "body_font",
    base_line_size = unit(0.1, "mm")
  ) +
  theme(
    plot.title = element_text(size = 30)
  )
g

Agrupamiento de distanciasminimas y escala de color

#rtistry280
library(cptcity)
#http://soliton.vm.bytemark.co.uk/pub/cpt-city/notes/formats.html
#http://soliton.vm.bytemark.co.uk/pub/cpt-city/index.html
pals <- find_cpt("ncl")
pals <- find_cpt("seq")
#scale_fill_gradientn(name = "medallas",colours = cpt(n = 10, "cb_seq_PuBuGn_09"),breaks = round(breaks, 0))+

breaks <- classInt::classIntervals(
  distanciasminimas_sf$mintime,
  n = 10,
  #style = "pretty"
  #style = "quantile"
  #style = "kmeans"
  #style = "fisher"
  style = "jenks"
)$brks

breaks <- classInt::classIntervals(distanciasminimas_sf$mintime, n=10, style='jenks')
names(breaks)

[1] “var” “brks”

distanciasminimas_sf$clasif<-cut(distanciasminimas_sf$mintime, breaks = breaks$brks, labels=as.character(1:10))
distanciasminimas_sf$clasif2<-cut(distanciasminimas_sf$mintime, breaks = breaks$brks, include.lowest = TRUE)



distanciasminimast<-distanciasminimas_sf%>% 
  group_by(clasif, clasif2) %>%
  summarize(Cuenta = n()) 

plot(distanciasminimast)

g <- ggplot() +
  geom_sf(data = tor, fill = NA) +
  geom_sf(data = recorteviasosm_sf, fill = NA, linewidth = 0.1) +
  geom_sf(
    data = recortedenue_sf,
   # alpha = 0.3,
    size = 1,
    color="red"
  ) +
  

  geom_sf(data = distanciasminimast,aes(color=Cuenta))+
  scale_color_gradientn(colours = cpt(n = 10, "jjg_serrate_seq_srtYlGnBu05"))+
  
    #geom_sf_text_repel(data=distanciasminimas_sf,aes(label=mintime_str), size=2, color="navyblue")+

 labs(
    title = "Hospitales y puntos aleatorios"
  ) +
  theme_minimal(
    base_family = "body_font",
    base_line_size = unit(0.1, "mm")
  ) +
  theme(
    plot.title = element_text(size = 30)
  )
g

Se interpolarán los puntos y se creará una capa raster

Interpolation For this we prepare the following function. We use a nearest neighbor (k = 100) algorithm is used.

library(gstat)

interpolateSurface <- function(data, gridRes = 100){
  
  # Make new dataset, that will be spatial - sp:: class
  data_sp <- data %>% as.data.frame()
  # Make data as spatial object - spatialPointsDataFrame
  coordinates(data_sp) <- ~lon+lat
  # Define CRS - coordinate reference system
  proj4string(data_sp) <- crs("+init=epsg:4326")
  # Make it as of object of sf:: class
  crs <- "+init=epsg:4326"
  coords <- c("lon", "lat")
  # (sf: Simple feature collection)
  data_sf <- st_as_sf(data_sp, coords = coords, crs = crs)
  # WEB Mercator projection
  data_sp_mp <- spTransform(data_sp, CRSobj = "+init=epsg:3857")

  # Bounding box, resolution and grid for interpolation
  boxx = data_sp_mp@bbox
  deltaX = as.integer((boxx[1,2] - boxx[1,1]) + 1.5)
  deltaY = as.integer((boxx[2,2] - boxx[2,1]) + 1.5)
  gridSizeX = deltaX / gridRes
  gridSizeY = deltaY / gridRes
  grd = GridTopology(boxx[,1], c(gridRes,gridRes), c(gridSizeX,gridSizeY))
  pts = SpatialPoints(coordinates(grd))
  proj4string(pts) <- crs("+init=epsg:3857")
  
  # Interpolate the grid cells Nearest neighbour
  r.raster <- raster::raster()
  extent(r.raster) <- extent(pts) # set extent
  res(r.raster) <- gridRes # 500 # set cell size
  crs(r.raster) <- crs("+init=epsg:3857") # set CRS
  gs <- gstat(formula = mintime~1, 
              locations = data_sp_mp, 
              nmax = 100, 
              set = list(idp = 0))
  nn <- interpolate(r.raster, gs)
  
  # Result rasters - surfaces
  data_list <- list(data_sf, nn)
  names(data_list) <- c("data", "nn")
  return(data_list)
}


#Now that we have a function and know about the computation times, let’s run it one last and definitive time with #gridRes = 500 and save it for later use.

#mindistances <- distanciasminimas %>% 
#  distinct()
mindistances <- distanciasminimas 

minraster <- interpolateSurface(mindistances, gridRes = 500)

[inverse distance weighted interpolation]

rutrast<-"C:/Users/cguer/Documents/Claudia/Midropbox/Investigacion y escritos/mapasnacionales//CEM_V3_20170619_R15_E05_TIF/"
rst_mosaic <-raster(paste0(rutrast,'Coahuila_r15m.tif'))

sp::proj4string(rst_mosaic) <- CRS("+init=epsg:3857")
tor <- tor |> st_transform("EPSG:3857")  
coords = as.data.frame(unique(st_coordinates(tor)))
lonmin<-min(coords$X)
lonmax<-max(coords$X)
latmin<-min(coords$Y)
latmax<-max(coords$Y)



full_extent <- extent(lonmin, latmin,
                      lonmax, latmax)

recorterastertor <- crop(rst_mosaic, full_extent, snap="out")      

minraster<- crop(minraster$nn, full_extent, snap="out")      

minraster_shaped_aggr <- aggregate(minraster, fact = 1)

minraster_shaped_rtp <- rasterToPolygons(minraster_shaped_aggr)

minraster_shaped_rtp <- spTransform(minraster_shaped_rtp, 
                                    CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))

minraster_shaped_rtp@data$id <- 1:nrow(minraster_shaped_rtp@data)



minraster_shaped_fort <- fortify(minraster_shaped_rtp, 
                                 data = minraster_shaped_rtp@data)

minraster_shaped_fort <- merge(minraster_shaped_fort, 
                               minraster_shaped_rtp@data, 
                               by.x = 'id', 
                               by.y = 'id')



tor= st_transform(tor, crs = 4326)


recorteviasosm_sf= st_transform(recorteviasosm_sf, crs = 4326)
recortedenue_sf= st_transform(recortedenue_sf, crs = 4326)

g <- ggplot() +
  geom_sf(data = tor, fill = NA) +
  geom_sf(data = recorteviasosm_sf, fill = NA, linewidth = 0.1) +
  geom_sf(
    data = recortedenue_sf,
   # alpha = 0.3,
    size = 1,
    color="red"
  ) +
  
geom_polygon(data = minraster_shaped_fort, 
               aes(x = long, y = lat, group = group, fill = var1.pred, alpha = var1.pred), 
               # alpha = 0.8, 
               size = 0)+
 scale_fill_gradient(low = "white", high = "red") +
  scale_alpha_continuous(guide = "none", range = c(0.1,1))+

 labs(
    title = "Hospitales y puntos aleatorios"
  ) +
  theme_minimal(
    base_family = "body_font",
    base_line_size = unit(0.1, "mm")
  ) +
  theme(
    plot.title = element_text(size = 30)
  )
g

library(ggmap)

coords = as.data.frame(unique(st_coordinates(tor)))
lonmin<-min(coords$X)
lonmax<-max(coords$X)
latmin<-min(coords$Y)
latmax<-max(coords$Y)
library(ggmap)
zoom=10
#maptype = "alidade_smooth"
register_stadiamaps("fff6371d-02f3-48d7-9c96-b75506e41250")
bbox<-c(lonmin,latmin,lonmax,latmax)

  mapfon<-ggmap(get_stadiamap(bbox,maptype = "alidade_smooth", 14))

  mapfon+
     geom_sf(data =tor, inherit.aes = FALSE,fill=NA)+
     geom_sf(data = recorteviasosm_sf, inherit.aes = FALSE,fill = NA, linewidth = 0.1) +
  geom_sf(
    data = recortedenue_sf,
   # alpha = 0.3,
   inherit.aes = FALSE,
    size = 1,
    color="red"
  )+
  
geom_polygon(data = minraster_shaped_fort, 
               aes(x = long, y = lat, group = group, fill = var1.pred, alpha = var1.pred), 
               # alpha = 0.8, 
               size = 0)+
 scale_fill_gradient(low = "white", high = "red") +
  scale_alpha_continuous(guide = "none", range = c(0.1,1))+
  labs(
    title = "Qué tan lejos queda un hospital",
    subtitle=" Tiempo en minutos en puntos interpolados utilizando OSRM",
    fill="minutos para llegar al hospital más cercano"
  ) +
  theme_minimal(
    base_family = "body_font",
    base_line_size = unit(0.1, "mm")
  ) +
  theme(
    plot.title = element_text(size = 30)
  )