En este post voy a mostrar el uso de algunas librerias de R para el consumo de geoservicios. Además, veremos como visualizar las capas de geoservicios con leaflet y haremos algunas operaciones para apoyar la gestión local. Espero que este mundo los atrape tanto como a mí.
En muchas áreas de la vida social, económica, cultural, etc, la comprensión del espacio o “lo espacial” se vuelve fundamental para la comprensión de ciertas dinámicas o tendencias. En el área de salud esto se hace más evidente ya que en gran parte la salud de un individuo o de una comunidad está determinada por como esta configurado este espacio que lo rodea.
Por eso aquí vamos a trabajar un ejemplo del uso de los geoservicios para la gestión local en salud, pero ¿que son los geoservicios?.
Vamos a decir que son un Servicio Web específico que permite intercambiar información únicamente de componente geográfica. Para la generación y utilización de los estos se utilizan lenguajes específicos y protocolos estándares definidos.
Es la sigla para representar el término Infraestructura de Datos Espaciales, que cumple la función de permitir:
…acceder a datos, productos y servicios geoespaciales, publicados en internet bajo estándares y normas definidos, asegurando su interoperabilidad y uso, como así también la propiedad sobre la información por parte de los organismos que la publican y su responsabilidad en la actualización (IDERA)
En Argentina, este rol lo cumple IDERA, una comunidad geoespacial que involucra actores estatales, de la sociedad civil y privados.
Luego de la pequeña intro vamos a los datos. El ejemplo que utilizaremos para revisar el uso de las librerias corresponde a los geoservicios de la ciudad de Gualeguaychú, que además de ser mi ciudad de nacimiento, recientemente ha disponibilizado esta información para la comunidad en el siguiente link.
Realizaremos la lectura del servicio con la libreria sf, la la usaremos la lo largo de este post ya que nos permite una amplia gama de operaciones. Es importante saber que existen dos tipos de servicios WMS y WFS 1.
## [1] "gis:AG_Control"
## [2] "gis:AG_Nodos"
## [3] "gis:AG_Symbols"
## [4] "gis:AG_Tramos"
## [5] "gis:aerodromo"
## [6] "gis:aerodromo_poligono"
## [7] "gis:turismo_alojamientos"
## [8] "gis:antenas"
## [9] "gis:antiguos_nombres_de_calles"
## [10] "gis:area_de_fabricacion_y_procesamiento"La función st_layers de sf nos permite explorar todas las capas disponibles en el geoservicio. En nuestro caso tomaremos 3:
Establecimientos de Salud (gis:est_salud)
Areas programáticas (gis:areas_programaticas)
Radios Censales (gis:radios_censales)
Miremos ahora 👀 las capas en un mapa, y usaremos en esta oportunidad la librería leaflet que nos permite acceder a mapas interactivos muy poderosos. Además, utilizaremos la libreria leaflet.extras2 que nos da la posibilidad de ver características extras a nuestras capas.
Usaremos Argenmapcomo tesela 2 base y luego iremos agregando las respectivas capas. Además, generaremos un control de capas que nos va a permitir seleccionar cuál o cuáles queremos ver.
leaflet() %>%
addTiles(urlTemplate = "https://wms.ign.gob.ar/geoserver/gwc/service/tms/1.0.0/mapabase_gris@EPSG%3A3857@png/{z}/{x}/{-y}.png",
tileOptions(tms = TRUE,maxZoom = 14), attribution = '<a target="_blank" href="https://www.ign.gob.ar/argenmap/argenmap.jquery/docs/#datosvectoriales" style="color: black; text-decoration: underline; font-weight: normal;">Datos IGN Argentina // OpenStreetMap</a>',
group = "Argenmap") %>% #Aquí agregamos Argenmap
addWMS(
geo_gchu,
layers = "radios_censales",
options = WMSTileOptions(format = "image/png", transparent = TRUE,info_format = "text/html", tiled=FALSE),
group = "Radios Censales"
) %>% #Agregamos layer de Radios Censales
addWMS(
geo_gchu,
layers = "areas_programaticas",
options = WMSTileOptions(format = "image/png", transparent = TRUE,info_format = "text/html", tiled=FALSE),
group = "Areas Programáticas"
) %>% #Agregamos layer de Areas Programáticas
addWMS(
geo_gchu,
layers = "est_salud",
options = WMSTileOptions(format = "image/png", transparent = TRUE,info_format = "text/html", tiled=FALSE),
group = "Est de salud"
) %>% #Agregamos layer de Establecimientos de Salud
addLayersControl(
baseGroups = "Argenmap",
overlayGroups = c("Areas Programáticas","Est de salud","Radios Censales"),
options = layersControlOptions(collapsed = FALSE)) %>%
setView(lng = -58.52597, lat = -33.00606,zoom = 11)
Si cliquean sobre el área programática pueden obtener información del responsable del área, ubicación, teléfono. Si en cambio cliquean el radio censal les devuelve información de la cantidad de hogares, personas e incluso el área 🙌🏼.
El acceso al sistema de salud puede ser abordado desde distintas perspectivas, barreras culturales, geográficas, presencia de medios de transporte, etc.
En este ejercicio intentaré abordar la accesibilidad que tiene cada área programática a un centro de salud, para poder determinar en términos de tiempo y distancia cuál tiene una mayor accesibilidad y cuál menor.
Antes de avanzar, cuándo hablamos de área programática hacemos referencia al área de influencia de un determinado centro asistencial que se expresa, aunque no únicamente, en un polígono geográfico. La idea central es que las personas que viven en esa área se referencian con el centro de salud o los centros de salud contenidos en ésta 3.
La metodología que voy a usar para determinar la distancia por área programática al centro de salud, es asignar los radios censales a cada uno de los poligonos que representan el área programática.
Luego, voy a calcular los centroides de cada radio censal y voy a medir la distancia en tiempo (a pie) y en kilómetros hasta el centro o los centros de salud que contiene el area.
Resulta un poco confuso 🙄, veamos el código.
#Extraemos las areas programaticas
areas_programaticas <- st_read("WFS:https://geo.gualeguaychu.gov.ar/geoserver/gis/wms?Request=GetCapabilities","gis:areas_programaticas")
## Reading layer `gis:areas_programaticas' from data source
## `WFS:https://geo.gualeguaychu.gov.ar/geoserver/gis/wms?Request=GetCapabilities'
## using driver `WFS'
## Simple feature collection with 9 features and 5 fields
## Geometry type: MULTISURFACE
## Dimension: XY
## Bounding box: xmin: 5628553 ymin: 6342554 xmax: 5641830 ymax: 6350683
## Projected CRS: POSGAR 98 / Argentina 5#Extraemos los establecimientos de salud
est_salud <- st_read("WFS:https://geo.gualeguaychu.gov.ar/geoserver/gis/wms?Request=GetCapabilities", "gis:est_salud")
## Reading layer `gis:est_salud' from data source
## `WFS:https://geo.gualeguaychu.gov.ar/geoserver/gis/wms?Request=GetCapabilities'
## using driver `WFS'
## Simple feature collection with 11 features and 17 fields
## Geometry type: MULTIPOINT
## Dimension: XY
## Bounding box: xmin: -58.54233 ymin: -33.03131 xmax: -58.50246 ymax: -32.99495
## Geodetic CRS: WGS 84#Extraemos los radios censales
radios_censales <- st_read("WFS:https://geo.gualeguaychu.gov.ar/geoserver/gis/wms?Request=GetCapabilities","gis:radios_censales")
## Reading layer `gis:radios_censales' from data source
## `WFS:https://geo.gualeguaychu.gov.ar/geoserver/gis/wms?Request=GetCapabilities'
## using driver `WFS'
## Simple feature collection with 95 features and 10 fields
## Geometry type: MULTISURFACE
## Dimension: XY
## Bounding box: xmin: 5625130 ymin: 6334522 xmax: 5675083 ymax: 6353581
## Projected CRS: POSGAR 98 / Argentina 5#homogeneizamos las proyecciones
# Centros de salud
csalud <- est_salud %>% st_drop_geometry()
csalud <- st_as_sf(csalud[c("fid","nombre","longitud","latitud")], coords = c("longitud", "latitud"),
crs = 4326, remove= FALSE) %>% st_transform(crs = 4326)
est_salud <- st_transform(est_salud, 4326)
# Areas programaticas
areas_programaticas <- st_transform(areas_programaticas, 4326)
areas_programaticas <- st_cast(areas_programaticas, "GEOMETRYCOLLECTION") %>% st_collection_extract("POLYGON")
# Radios censales
radios_censales <- st_transform(radios_censales, 4326)
radios_censales <- st_cast(radios_censales, "GEOMETRYCOLLECTION") %>% st_collection_extract("POLYGON")
Primero utilizamos la función st_centroid de sf y lo agregamos como layers a nuestro mapa.
radios_censales$centroides <- radios_censales %>% st_centroid() %>%
st_geometry()
centroides <- radios_censales %>%
select(gml_id)%>%
st_centroid() %>%
st_as_sf()
leaflet() %>%
addTiles(urlTemplate = "https://wms.ign.gob.ar/geoserver/gwc/service/tms/1.0.0/mapabase_gris@EPSG%3A3857@png/{z}/{x}/{-y}.png",
tileOptions(tms = TRUE,maxZoom = 14), attribution = '<a target="_blank" href="https://www.ign.gob.ar/argenmap/argenmap.jquery/docs/#datosvectoriales" style="color: black; text-decoration: underline; font-weight: normal;">Datos IGN Argentina // OpenStreetMap</a>',
group = "Argenmap") %>% #Aquí agregamos Argenmap
addWMS(
geo_gchu,
layers = "radios_censales",
options = WMSTileOptions(format = "image/png", transparent = TRUE,info_format = "text/html", tiled=FALSE),
group = "Radios Censales"
) %>% #Agregamos layer de Radios Censales
addWMS(
geo_gchu,
layers = "areas_programaticas",
options = WMSTileOptions(format = "image/png", transparent = TRUE,info_format = "text/html", tiled=FALSE),
group = "Areas Programáticas"
) %>% #Agregamos layer de Areas Programáticas
addWMS(
geo_gchu,
layers = "est_salud",
options = WMSTileOptions(format = "image/png", transparent = TRUE,info_format = "text/html", tiled=FALSE),
group = "Est de salud"
) %>% #Agregamos layer de Establecimientos de Salud
addCircles(data= centroides,group = "centroides",color = "black")%>% #Agregamos el circulo con los centroides
addLayersControl(
baseGroups = "Argenmap",
overlayGroups = c("Areas Programáticas","Est de salud","Radios Censales","centroides"),
options = layersControlOptions(collapsed = FALSE)) %>%
setView(lng = -58.52597, lat = -33.00606,zoom = 11)
En esta etapa vamos a asignar los centros de salud y los radios censales a cada área programática, para luego poder hacer los cálculos. Utilizamos la función st_join de sf.
join <- st_join(areas_programaticas,csalud,join= st_intersects)
Veamos que los centros de salud han sido asignado a areas programatica
Hacemos la asignación de los radios censales a las áreas programáticas.
join2 <- st_join(join,radios_censales,join= st_intersects)
Revisamos el resultado.
Llegamos a una etapa por demás interesante, en la que obtendremos las distancias entre los centroides y los centros de salud. Para ello, utilizo el paquete OSRM que se basa en el proyecto que lleva el mismo nombre y que es un Servicio de enrutamiento basado en datos de OpenStreetMap. Vale decir, que el resultado de las mediciones es la distancia en tiempo a pie 🚶y en kilómetros.
Vale aclarar que en el código de lectura de los radios censales se excluyen los radios 11 (Pueblo Belgrano) y 2, ya que abarcan territorio muy por fuera del ejido.
Para tener las distancias, voy a:
crear un data frame con los puntos de los establecimientos de salud y los centroides por area programatica.
separar el data frame anterior en uno de origen (from) y uno de destino (to).
obtener las distancias
unificar los dataset.
distancias_x_cs <- cbind(distancias[,c(1,2)] %>% st_drop_geometry(),df)
join %>%
st_drop_geometry()%>%
select(gml_id,fid)%>%
left_join(distancias_x_cs, by= "fid") %>%
group_by(gml_id)%>%
filter(!is.na(duration)) %>%
summarise(min= round(min(duration),1),
prom= round(mean(duration),1),
mediana= round(median(duration),1),
max= round(max(duration),1)) %>%
arrange(prom)%>%
kableExtra::kbl(format.args = list(big.mark = ".",
decimal.mark= ",")) %>%
kableExtra::kable_classic_2()
| gml_id | min | prom | mediana | max |
|---|---|---|---|---|
| areas_programaticas.4 | 1,0 | 2,5 | 2,1 | 7,2 |
| areas_programaticas.12 | 0,8 | 2,6 | 2,5 | 5,4 |
| areas_programaticas.11 | 0,6 | 2,8 | 3,0 | 4,3 |
| areas_programaticas.3 | 0,1 | 3,0 | 3,0 | 5,2 |
| areas_programaticas.9 | 0,4 | 3,7 | 3,8 | 10,1 |
| areas_programaticas.8 | 1,1 | 3,9 | 3,9 | 6,5 |
| areas_programaticas.2 | 2,1 | 4,6 | 3,6 | 11,9 |
| areas_programaticas.5 | 1,0 | 5,1 | 2,8 | 11,9 |
join %>%
st_drop_geometry()%>%
select(gml_id,fid)%>%
left_join(distancias_x_cs, by= "fid") %>%
group_by(gml_id)%>%
filter(!is.na(duration)) %>%
summarise(min= round(min(distance),1),
prom= round(mean(distance),1),
mediana= round(median(distance),1),
max= round(max(distance),1)) %>%
arrange(prom)%>%
kableExtra::kbl(format.args = list(big.mark = ".",
decimal.mark= ",")) %>%
kableExtra::kable_classic_2()
| gml_id | min | prom | mediana | max |
|---|---|---|---|---|
| areas_programaticas.4 | 0,4 | 1,1 | 0,8 | 3,1 |
| areas_programaticas.11 | 0,1 | 1,3 | 1,4 | 2,5 |
| areas_programaticas.12 | 0,3 | 1,3 | 1,2 | 3,0 |
| areas_programaticas.3 | 0,1 | 1,4 | 1,4 | 2,5 |
| areas_programaticas.8 | 0,4 | 1,7 | 1,6 | 2,6 |
| areas_programaticas.9 | 0,2 | 1,9 | 1,9 | 5,0 |
| areas_programaticas.2 | 0,8 | 2,0 | 1,4 | 6,1 |
| areas_programaticas.5 | 0,5 | 3,4 | 1,2 | 9,0 |
✔️ En relación al objetivo principal del ejercicio, puede observase que el área programática 4, es la que tiene mayor accesibilidad con un promedio de 2,5 minutos y 1,1 Km al centro de salud San Francisco. En el opuesto el área programática 5 mostró la mayor distancia con un promedio de 5,1 minutos y 3,4 Km al centro de salud Suburbio Sur. Sin dudas un promedio de 5 minutos al centro de salud en el área con menor accesibilidad da cuentas de una amplia cobertura del sistema de salud.
✔️En lo que respecta a la metodología sin dudas muestra limitaciones. Una de ellas es el hecho de que algunos radios censales caen en mas de un área programática. Por ej: el radio censal 13 cae dentro de Villa Maria, Suburbio Sur, Médanos. Esto implica que un habitante de esta zona cuenta con mayor oferta. Otra limitación es haber tomado los centroides, se podría mejorar esto tomando puntos al azar dentro de los radios censales.
✔️Por ultimo espero haber podido mostrar la potencialidad que tiene estos datos y la necesidad de bregar porque los mismos sigan siendo de acceso publico y se mantengan los estándares de calidad e interoperabilidad.
WMS: permite la visualización de información geográfica a partir de una representación de ésta, de una imagen del mundo real para un área solicitada por el usuario. Puede organizarse en una o más capas de datos que pueden visualizarse u ocultarse una a una. WFS : permite el acceso y consulta de los atributos de un vector (feature) que representa información geográfica como un río, una ciudad o un lago, con una geometría descrita por un conjunto de coordenadas. El servicio WFS permite no solo visualizar la información tal y como permite un WMS, sino también consultarla y editarla libremente↩︎
Las teselas vectoriales son paquetes de datos geográficos, empaquetados en «mosaicos» predefinidos de forma aproximadamente cuadrada para su transferencia a través de la web.↩︎
No necesariamente las personas buscan atención siguiendo esta racionalidad, ya que existen múltiples razones para orientarse por otro establecimiento, como puede ser la complejidad del tratamiento requerido.↩︎