Prueba 3
Lavezzolo
30/6/2020
library(tidyverse) # Paquete multiuso## -- Attaching packages ----------------------------------------------------------------------------------------- tidyverse 1.3.0 --## v ggplot2 3.3.2 v purrr 0.3.4
## v tibble 3.0.3 v dplyr 1.0.0
## v tidyr 1.1.0 v stringr 1.4.0
## v readr 1.3.1 v forcats 0.5.0## -- Conflicts -------------------------------------------------------------------------------------------- tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()library(sf) # Paquete clave para manipular datos espaciales## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1library(leaflet) # Uno de los paquetes para # Creamos un Data Frame con los datos necesarios
datos <- data.frame(lat = c(-34.714656, 51.532068), long = c(-58.785999,
-0.177305), ubicacion = c("UTDT", "Abbey Road"))
# Lo convertimos a un objeto sf
puntosEspaciales <- st_as_sf(datos, coords = c("long", "lat"),
crs = 4326)
# st_distance() nos permite encontrar la diferencia en la
# unidad que diga el CRS (sistema de coordenadas de
# referencia)
st_distance(puntosEspaciales)/1000 # En kilómetros## Units: [m]
## [,1] [,2]
## [1,] 0.00 11131.51
## [2,] 11131.51 0.00leaflet(puntosEspaciales) %>%
addTiles() %>%
addMarkers()# Cargamos la librería de SF
library(sf)
calles <- read_sf("D:/Documents/01-Ditella/Instrumentos de Analisis Urbano II/callejero-rar/callejero.shp")st_crs(calles)## Coordinate Reference System:
## User input: Argentina_GKBsAs
## wkt:
## PROJCRS["Argentina_GKBsAs",
## BASEGEOGCRS["Campo Inchauspe",
## DATUM["Campo Inchauspe",
## ELLIPSOID["International 1924",6378388,297,
## LENGTHUNIT["metre",1]],
## ID["EPSG",6221]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["Degree",0.0174532925199433]]],
## CONVERSION["unnamed",
## METHOD["Transverse Mercator",
## ID["EPSG",9807]],
## PARAMETER["Latitude of natural origin",-34.6297166,
## ANGLEUNIT["Degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",-58.4627,
## ANGLEUNIT["Degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["Scale factor at natural origin",0.999998,
## SCALEUNIT["unity",1],
## ID["EPSG",8805]],
## PARAMETER["False easting",100000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",100000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]]],
## CS[Cartesian,2],
## AXIS["(E)",east,
## ORDER[1],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]],
## AXIS["(N)",north,
## ORDER[2],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]]]Si hay que descargar archivo, hay que cambiar el txt!! cuando lo guardas(hacerlo como guardar como) poner gjson
radiosCensales <- read_sf("D:/Documents/01-Ditella/Instrumentos de Analisis Urbano II/caba_radios_censales.geojson")
glimpse(radiosCensales)## Rows: 3,554
## Columns: 9
## $ RADIO_ID <chr> "1_1_1", "1_12_1", "1_12_10", "1_12_11", "1_12_2", "1_1...
## $ BARRIO <chr> "RETIRO", "SAN NICOLAS", "SAN NICOLAS", "SAN NICOLAS", ...
## $ COMUNA <chr> "1", "1", "1", "1", "1", "1", "1", "1", "1", "1", "1", ...
## $ POBLACION <dbl> 336, 341, 296, 528, 229, 723, 393, 600, 472, 786, 329, ...
## $ VIVIENDAS <dbl> 82, 365, 629, 375, 445, 744, 341, 505, 504, 546, 275, 8...
## $ HOGARES <dbl> 65, 116, 101, 136, 129, 314, 209, 275, 202, 347, 129, 3...
## $ HOGARES_NBI <dbl> 19, 25, 1, 7, 16, 104, 110, 32, 49, 89, 15, 57, 1, 1, 2...
## $ AREA_KM2 <dbl> 1.79899705, 0.01856469, 0.04438025, 0.36634000, 0.01836...
## $ geometry <MULTIPOLYGON [°]> MULTIPOLYGON (((-58.37189 -..., MULTIPOLYG...hospitales <- read_delim("D:/Documents/01-Ditella/Instrumentos de Analisis Urbano II/Hospitales.csv",delim = ";") ## Parsed with column specification:
## cols(
## .default = col_character(),
## long = col_double(),
## lat = col_double(),
## calle_altura = col_double(),
## calle_altura_2 = col_double(),
## codigo_postal = col_double(),
## camas_medicina = col_double(),
## camas_cirugía = col_double(),
## camas_pediatría = col_double(),
## camas_tocoginecología = col_double(),
## camas_urgencia = col_double(),
## camas_otra = col_double(),
## camas_totales = col_double()
## )## See spec(...) for full column specifications.hospitales <- st_as_sf(hospitales,coords=c("long","lat"), crs=4326)manzanas <- read_sf("D:/Documents/01-Ditella/Instrumentos de Analisis Urbano II/manzanas.geojson")hospitales <- st_transform(hospitales, crs="+proj=tmerc +lat_0=-34.6297166 +lon_0=-58.4627 +k=1 +x_0=100000 +y_0=100000 +ellps=intl +units=m +no_defs")coberturaHospitales <- st_buffer(hospitales,dist = 1000)
coberturaHospitales <- coberturaHospitales %>% summarise(cobertura=TRUE)coberturaHospitales1 <- coberturaHospitales%>%
st_transform(4326)leaflet(coberturaHospitales1) %>%
addTiles() %>%
addPolygons()avenidas <- calles %>% filter(tipo_c == "AVENIDA")plot(avenidas[,c("sentido")])
manzanas <- st_transform(manzanas,crs = st_crs(hospitales))manzanasCentroides <- st_centroid(manzanas)## Warning in st_centroid.sf(manzanas): st_centroid assumes attributes are constant
## over geometries of xggplot() +
geom_sf(data = manzanas) +
geom_sf(data = manzanasCentroides,color="red", size= 0.001)
avenidas <- st_transform(avenidas,crs = st_crs(hospitales))distanciaAvenidas <- st_distance(manzanasCentroides,avenidas)# 12.520 filas (centroides de manzanas) y 6,758 columnas (tramos de avenidas)
dim(distanciaAvenidas)## [1] 12520 6686# 1 significa filas, 2 columnas. functon(x) min(x) significa que para cada fila devuelva el valor mínimo
avenidaMasCercana <- apply(distanciaAvenidas,1,function(x) min(x))# Rendondeamos
avenidaMasCercana <- round(avenidaMasCercana,0)manzanas <- manzanas %>% mutate(distanciaAvenida=avenidaMasCercana)manzanasCentroides <- st_join(manzanasCentroides,coberturaHospitales)manzanasCentroides <- manzanasCentroides %>% mutate(cobertura=ifelse(is.na(cobertura),FALSE,cobertura))manzanasCentroides <- manzanasCentroides %>%
mutate(cobertura= ifelse(is.na(cobertura), FALSE, cobertura))radiosCensales <- radiosCensales %>% mutate(densidadPob=POBLACION/AREA_KM2)
radiosCensales <- radiosCensales %>% st_transform(radiosCensales,crs=st_crs(hospitales))
manzanasCentroides <- st_join(manzanasCentroides,radiosCensales)# Elegimos solo las variables que queremos unir, antes lo convertimos en data frame para poder perder la columna geometry
manzanasCentroides <- manzanasCentroides %>%
as.data.frame() %>%
select(SM,densidadPob,HOGARES_NBI,cobertura)
manzanas <- left_join(manzanas,manzanasCentroides,by="SM")#Paso 5
# Queremos que nos muestre en que porcentaje de estos está cada observación...
quiebres <- c(0,0.2,0.4,0.6,0.8,1)
manzanas <- manzanas %>%
mutate(cat_densidad=cut(densidadPob,breaks = quantile(densidadPob,quiebres,na.rm = TRUE ),include.lowest = TRUE),
cat_NBI=cut(HOGARES_NBI,breaks = quantile(HOGARES_NBI,quiebres,na.rm = TRUE ),include.lowest = TRUE),
cat_distanciaAv=cut(-distanciaAvenida,breaks = quantile(-distanciaAvenida,quiebres,na.rm = TRUE ),include.lowest = TRUE))ggplot() +
geom_sf(data=manzanas %>% filter(!is.na(cat_NBI)) ,aes(fill=cat_NBI), color=NA) +
scale_fill_viridis_d() +
theme_minimal() +
coord_sf(datum=NA)
# En el caso de cobertura convertira 0 cuando era FALSE y 1 cuando era TRUE
manzanas <- manzanas %>%
mutate(cat_densidad=as.numeric(cat_densidad),
cat_NBI=as.numeric(cat_NBI),
cat_distanciaAv=as.numeric(cat_distanciaAv),
cat_cobertura=as.numeric(!cobertura))manzanas <- manzanas %>%
mutate(IDS=cat_densidad*0.1+cat_NBI*0.3+cat_distanciaAv*0.1+cat_cobertura*0.5,
IDS=ifelse(IDS>quantile(IDS,probs = 0.9,na.rm = TRUE),TRUE,FALSE))# Podemos agrupar a los radioscensales por barrio para que nos queden los polígonos de los barrios.
# Tambien es posible bajarlos directamente desde la página del GCBA
barrios <- radiosCensales %>% group_by(BARRIO) %>% summarise(n())## `summarise()` ungrouping output (override with `.groups` argument)ggplot() +
geom_sf(data=barrios) +
geom_sf(data=manzanas %>% filter(!is.na(IDS)) ,aes(fill=IDS), color=NA) +
scale_fill_manual(values = c(NA,"#f03b20")) +
guides(fill=FALSE) +
theme_minimal() +
coord_sf(datum=NA) +
labs(title="Índice de Demanda de Salud",
subtitle="Ciudad de Buenos Aires")