Mapas de objetos espaciales

Las visualizaciones basadas en mapas que describen los resultados de la información espacial a menudo utilizan archivos con formato shape, los cuales son comúnmente empleados en sistemas de información geográfica (GIS) para representar datos geoespaciales. Estos archivos constan de un conjunto de archivos que contienen información sobre elementos geográficos, tales como puntos, líneas o polígonos, junto con sus atributos asociados. Para este caso se usará un archivo de la librería rnaturalearth.

la librería rnaturalearth función ne_countries() con los argumentos type = “countries”, scale = “medium” y returnclass = “sf” genera información espacial de todos los países del mundo. Aquí hay una descripción de lo que cada argumento hace:

type = “countries”: Indica que se deben devolver datos de países.
scale = “medium”: Especifica la escala de los datos que se deben devolver. En este caso, “medium” probablemente se refiere a una resolución media de los datos geoespaciales.
returnclass = “sf”: Solicita que los datos se devuelvan en una clase especial llamada “sf”, que generalmente se utiliza con objetos espaciales en R. “sf” es una abreviatura de “simple features”, un estándar de representación de datos geoespaciales en R.

Entonces, en resumen, esta función devolverá datos espaciales en formato “sf” (simple features) para todos los países del mundo, con una resolución de escala media. Para este ejemplo usemos el mapa de Italia.

library(rnaturalearth)
map1 <- ne_countries(type = "countries", country = "Italy",
                     scale = "medium", returnclass = "sf")

#SALIDA DEL OBJETO ESPACIAL

##Simple feature collection with 1 feature and 168 fields
##Geometry type: MULTIPOLYGON
##Dimension:     XY
##Bounding box:  xmin: 6.627734 ymin: 36.68784 xmax: 18.48584 ymax: 47.08213
##Geodetic CRS:  WGS 84
##         featurecla scalerank labelrank sovereignt sov_a3 adm0_dif level              type tlc admin adm0_a3
## Admin-0 country         1         2      Italy    ITA        0     2 Sovereign country   1 Italy     ITA
##  geou_dif geounit gu_a3 su_dif subunit su_a3 brk_diff  name name_long brk_a3 brk_name brk_group abbrev postal
##        0   Italy   ITA      0   Italy   ITA        0 Italy     Italy    ITA    Italy      <NA>  Italy      I 
## formal_en formal_fr name_ciawf note_adm0 note_brk name_sort name_alt mapcolor7 mapcolor8 mapcolor9
## Italian Republic      <NA>      Italy      <NA>     <NA>     Italy     <NA>         6         7         8

names(map1)

##   [1] "featurecla" "scalerank"  "labelrank"  "sovereignt" "sov_a3"    
##   [6] "adm0_dif"   "level"      "type"       "tlc"        "admin"     
##  [11] "adm0_a3"    "geou_dif"   "geounit"    "gu_a3"      "su_dif"    
##  [16] "subunit"    "su_a3"      "brk_diff"   "name"       "name_long" 
##  [21] "brk_a3"     "brk_name"   "brk_group"  "abbrev"     "postal"    
##  [26] "formal_en"  "formal_fr"  "name_ciawf" "note_adm0"  "note_brk"  
##  [31] "name_sort"  "name_alt"   "mapcolor7"  "mapcolor8"  "mapcolor9" 
##  [36] "mapcolor13" "pop_est"    "pop_rank"   "pop_year"   "gdp_md"    
##  [41] "gdp_year"   "economy"    "income_grp" "fips_10"    "iso_a2"    
##  [46] "iso_a2_eh"  "iso_a3"     "iso_a3_eh"  "iso_n3"     "iso_n3_eh" 
##  [51] "un_a3"      "wb_a2"      "wb_a3"      "woe_id"     "woe_id_eh" 
##  [56] "woe_note"   "adm0_iso"   "adm0_diff"  "adm0_tlc"   "adm0_a3_us"
##  [61] "adm0_a3_fr" "adm0_a3_ru" "adm0_a3_es" "adm0_a3_cn" "adm0_a3_tw"
##  [66] "adm0_a3_in" "adm0_a3_np" "adm0_a3_pk" "adm0_a3_de" "adm0_a3_gb"
##  [71] "adm0_a3_br" "adm0_a3_il" "adm0_a3_ps" "adm0_a3_sa" "adm0_a3_eg"
##  [76] "adm0_a3_ma" "adm0_a3_pt" "adm0_a3_ar" "adm0_a3_jp" "adm0_a3_ko"
##  [81] "adm0_a3_vn" "adm0_a3_tr" "adm0_a3_id" "adm0_a3_pl" "adm0_a3_gr"
##  [86] "adm0_a3_it" "adm0_a3_nl" "adm0_a3_se" "adm0_a3_bd" "adm0_a3_ua"
##  [91] "adm0_a3_un" "adm0_a3_wb" "continent"  "region_un"  "subregion" 
##  [96] "region_wb"  "name_len"   "long_len"   "abbrev_len" "tiny"      
## [101] "homepart"   "min_zoom"   "min_label"  "max_label"  "label_x"   
## [106] "label_y"    "ne_id"      "wikidataid" "name_ar"    "name_bn"   
## [111] "name_de"    "name_en"    "name_es"    "name_fa"    "name_fr"   
## [116] "name_el"    "name_he"    "name_hi"    "name_hu"    "name_id"   
## [121] "name_it"    "name_ja"    "name_ko"    "name_nl"    "name_pl"   
## [126] "name_pt"    "name_ru"    "name_sv"    "name_tr"    "name_uk"   
## [131] "name_ur"    "name_vi"    "name_zh"    "name_zht"   "fclass_iso"
## [136] "tlc_diff"   "fclass_tlc" "fclass_us"  "fclass_fr"  "fclass_ru" 
## [141] "fclass_es"  "fclass_cn"  "fclass_tw"  "fclass_in"  "fclass_np" 
## [146] "fclass_pk"  "fclass_de"  "fclass_gb"  "fclass_br"  "fclass_il" 
## [151] "fclass_ps"  "fclass_sa"  "fclass_eg"  "fclass_ma"  "fclass_pt" 
## [156] "fclass_ar"  "fclass_jp"  "fclass_ko"  "fclass_vn"  "fclass_tr" 
## [161] "fclass_id"  "fclass_pl"  "fclass_gr"  "fclass_it"  "fclass_nl" 
## [166] "fclass_se"  "fclass_bd"  "fclass_ua"  "geometry"

1. Mapa de un país (Italia)

Inmersamente el objeto en formato sf tiene información adicional muy util para el desarrollo de análisis más sofisticados, entre los atributos más importantes se encuentran:

featurecla: Tipo de característica geográfica (por ejemplo, país, océano, etc.).
scalerank: Rango de la escala del país.
labelrank: Rango del etiquetado.
sovereignt: Nombre del país.
adm0_a3: Código del país según el estándar ISO A3.
pop_est: Estimación de la población.
gdp_md: Producto Interno Bruto (PIB) en millones de dólares.
continent: Continente al que pertenece el país.
region_un: Región según la ONU.
subregion: Subregión geográfica.
name: Nombre del país (usualmente el mismo que “sovereignt”).
name_long: Nombre largo del país.
iso_a2: Código del país según el estándar ISO A2 (código de dos letras).
iso_a3: Código del país según el estándar ISO A3 (código de tres letras).
iso_n3: Código del país según el estándar ISO N3 (código de tres dígitos).
geometry: Datos geométricos que representan la forma y límites del país en formato espacial.

Estos son algunos de los campos más destacados y útiles para trabajar con datos geográficos de países. Por supuesto, la relevancia de cada campo dependerá del contexto específico del análisis o proyecto.

Para acceder a alguna de las variables anteriormente mensionadas:

map1[,1]

## Simple feature collection with 1 feature and 1 field
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 6.627734 ymin: 36.68784 xmax: 18.48584 ymax: 47.08213
## Geodetic CRS:  WGS 84
##          featurecla                       geometry
## 138 Admin-0 country MULTIPOLYGON (((7.021094 45...

map1[,4]

## Simple feature collection with 1 feature and 1 field
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 6.627734 ymin: 36.68784 xmax: 18.48584 ymax: 47.08213
## Geodetic CRS:  WGS 84
##     sovereignt                       geometry
## 138      Italy MULTIPOLYGON (((7.021094 45...

Observemos la primera visualización de la información (Mapa de Italia):

plot(map1[,4])

2. Mapa de continentes

Una de las grandes ventajas de esta libreria es poder graficar cualquier continente, para esto se usa la función ne_countries() la cual permite seleccionar el continente al cual se quiere hacer enfasis.

africa  = ne_countries( continent = "Africa", returnclass = "sf")
class(africa)

## [1] "sf"         "data.frame"

africa$geometry      #Geometría del objeto

## Geometry set for 51 features 
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -17.62504 ymin: -34.81917 xmax: 51.13387 ymax: 37.34999
## Geodetic CRS:  WGS 84
## First 5 geometries:

La anterior salida es el resumen de los atributos geométricos del conjunto de datos de países de África:

Geometry set for 51 features: Indica que hay 51 polígonos (es decir, características geográficas, en este caso, países) en el conjunto de datos.
Geometry type: MULTIPOLYGON: Se refiere al tipo de geometría presente en el conjunto de datos. En este caso, se trata de polígonos múltiples, lo que significa que algunos países pueden tener múltiples partes o regiones.
Dimension: XY: Indica que las coordenadas geográficas están representadas en un sistema de coordenadas bidimensional (X e Y), lo que sugiere que es un mapa plano.
Bounding box: xmin, ymin, xmax, ymax: Estos valores representan las coordenadas mínimas y máximas en los ejes X (longitud) e Y (latitud) que definen el área total cubierta por los polígonos. Es una forma de describir el “rectángulo” que encierra todos los polígonos.
Geodetic CRS: WGS 84: Indica el sistema de referencia geodésico utilizado para las coordenadas geográficas. WGS 84 es un sistema comúnmente utilizado en sistemas de información geográfica y representa la forma de la Tierra como un elipsoide.
First 5 geometries: Muestra las geometrías de los primeros 5 polígonos en el conjunto de datos. Esto te dará una idea de cómo están representados los países en términos de sus límites y formas en el mapa.

terra::plot(africa)

2.1 Población estimada en Africa

Accediendo a la variable 10 que hace referencia a división política del continente:

africa_pob = africa[,10]
plot(africa_pob)

Ahora usando la librería ggplot2()

require(ggplot2)

# Crear un mapa básico
mapa <- ggplot(africa) + 
  geom_sf()

# Personalizar colores y estilos
mapa <- mapa +
  theme_bw() + # Utilizar un tema de fondo blanco y negro
  theme(axis.text = element_text(size = 8), # Ajustar el tamaño de las etiquetas del eje
        plot.title = element_text(hjust = 0.5)) # Centrar el título del gráfico

# Añadir título y etiquetas
mapa <- mapa +
  labs(title = "Mapa de países de África", 
       x = "Longitud", y = "Latitud")

# Mostrar el mapa
print(mapa)

2.2 Población estimada en Africa

El código que se muestra a continuación específicamente está usando el paquete rnaturalearth en R para cargar datos espaciales de los países de América en un objeto de clase sf. Luego, selecciona los datos de población de estos países y los almacena en un vector. Finalmente, visualiza estos datos de población en un gráfico de dispersión utilizando la función plot(). Este proceso permite explorar y comparar la población de los países sudamericanos de manera visual, lo que puede ser útil para análisis geográficos y demográficos.

library(ggplot2)
library(viridis)   #es como una escala de colores
ggplot(africa) + geom_sf(aes(fill = pop_est)) +  #pop_est para que coloree esta variable
  scale_fill_viridis() + theme_bw()

Ahora el de América del sur pero con una la nueva paleta de colores:

s_america=ne_countries( continent = "South America", returnclass = "sf")
ggplot(s_america) + geom_sf(aes(fill = pop_est)) +
  scale_fill_viridis() + theme_bw()

2.3 Cambiando etiquetas del mapa de América del Sur

#Transfomando las escalas de etiquetas
s_america2=s_america[,35]  #(Col 35 es la de poblacion estimada)
s_america_df=data.frame(s_america)
s_america2$log_pop_est=round(s_america_df$pop_est/1000000,1)
#s_america2$log_pop_est=log(s_america$pop_est/1000000)

g1=ggplot(s_america2) + geom_sf(aes(fill = log_pop_est)) +
  scale_fill_viridis() + theme_bw()
g1

2.4 Mapas dinámicos con ggplotly()

require(plotly)
ggplotly(g1)

2.5 Creando Nuevas Medidas (Producto Interno Bruto Percapita)

s_america$gdp_per_dolar = (s_america$gdp_md)/(s_america$pop_est)*1000000   #gdp_md_est

g2 = ggplot(s_america) + geom_sf(aes(fill = gdp_per_dolar)) +
  scale_fill_viridis() + theme_bw()

ggplotly(g2)

2.5 Quitando países en un mapa (Islas Malvinas)

## Quitando países
s_america$sovereignt

##  [1] "Argentina"      "Chile"          "United Kingdom" "Uruguay"       
##  [5] "Brazil"         "Bolivia"        "Peru"           "Colombia"      
##  [9] "Venezuela"      "Guyana"         "Suriname"       "Ecuador"       
## [13] "Paraguay"

pos=which(s_america$sovereignt=="United Kingdom")
pos

## [1] 3

Nuevo mapa sin Islas Malvinas

s_america3=s_america[-pos,]

g3=ggplot(s_america3) + geom_sf(aes(fill = gdp_per_dolar)) +
  scale_fill_viridis() + theme_bw()

ggplotly(g3)

3. Mapas Binarios

s_america3$gdp_valor=as.numeric(s_america3$gdp_per_dolar>10000)     #Valores 0 y 1nos
s_america3$gdp_valor2="malo"
pos2=which(s_america3$gdp_per_dolar>10000)
s_america3$gdp_valor2[pos2]="bueno"

ggplot(s_america3) + geom_sf(aes(fill = gdp_valor2))  + theme_bw()

# Cargar paquetes necesarios
library(ggplot2)
library(viridis)

# Crear un mapa basado en el PIB per cápita con una paleta de colores continua
ggplot(s_america3) + 
  geom_sf(aes(fill = gdp_per_dolar)) +
  scale_fill_viridis_c(name = "PIB per cápita", labels = scales::dollar_format()) + # Utilizar una paleta de colores continua
  theme_bw()

4. Graficando puntos en mapas.

library(paqueteMODELOS)
data("vivienda")
head(vivienda)

## # A tibble: 6 × 13
##      id zona    piso  estrato preciom areaconst parqueaderos banios habitaciones
##   <dbl> <chr>   <chr>   <dbl>   <dbl>     <dbl>        <dbl>  <dbl>        <dbl>
## 1  1147 Zona O… <NA>        3     250        70            1      3            6
## 2  1169 Zona O… <NA>        3     320       120            1      2            3
## 3  1350 Zona O… <NA>        3     350       220            2      2            4
## 4  5992 Zona S… 02          4     400       280            3      5            3
## 5  1212 Zona N… 01          5     260        90            1      2            3
## 6  1724 Zona N… 01          5     240        87            1      3            3
## # ℹ 4 more variables: tipo <chr>, barrio <chr>, longitud <dbl>, latitud <dbl>

datos<-data.frame(vivienda)
datos_sin_na <- na.omit(datos)

datos_sin_na <- na.omit(datos)

library("viridis")  
ggplot(data = datos_sin_na) +
  geom_point(mapping = aes(x = longitud, y = latitud, color = tipo)) +
  labs(title = "Viviendas en Cali", 
      x = "Longitud", 
      y = "Latitud", 
      color = "Tipo de Vivienda",
      hjust = 0.5) +
  scale_color_brewer(type = "qual", palette = 2) + 
  theme_bw() + 
  theme(plot.title = element_text(hjust = 0.5))

  #scale_color_viridis(discrete = TRUE, option = "E")   #D, E
  #theme(legend.position = "top")

library(leaflet)

crear_mapa <- function(datos, longitud, latitud) {
  leaflet() %>% 
    addTiles() %>% 
    addCircleMarkers(data = datos,
                     lng = ~ longitud,
                     lat = ~ latitud,
                     radius = 0.2,
                     color = "blue")
}

mapa <- crear_mapa(datos_sin_na, "longitud" , "latitud")

# Metadatos en HTML con estilos CSS para centrar el texto
metadatos <- '<div style="text-align: justify; font-size: 8px;margin-top: 1px;"><h2>Distribución de Viviendas en Cali</h2></div>'
metadatos_html <- htmltools::HTML(metadatos) # Convertir a objeto HTML
htmlwidgets::prependContent(mapa, metadatos_html) # Agregar metadatos al mapa

Distribución de Viviendas en Cali

Otra forma de poder visualizar puntos es por medio de las líneas de código:

library(sp)
viviendas_sp <- SpatialPointsDataFrame(coords = datos_sin_na[,12:13],  data = datos_sin_na)


#para poder graficar puntos (EJEMPLO2)
m <- leaflet(viviendas_sp) %>% addTiles() %>% addCircleMarkers(radius = 0.3,color = "red") 
m %>% addProviderTiles(providers$Stadia.AlidadeSmooth)

Cambiando el tipo de mapa OpenStreetMap.France

m <- leaflet(viviendas_sp) %>% addTiles() %>% addCircleMarkers(radius = 0.3,color = "red") 
m %>% addProviderTiles(providers$OpenStreetMap.France)

datos2=na.omit(datos)

leaflet() %>% addTiles() %>% 
  addCircleMarkers(lng = datos2[,12],lat = datos2[,13],clusterOptions = markerClusterOptions())

Para poder graficar puntos dados en Latitudes y Longitudes, se debe transformar los objetos que se encuentran en DataFrames en objetos espaciales tipo SpatialPointsDataFrame. Para hacer este proceso de transformar a SpatialPPoint se usa la librería sp de R. Obsérvese:

library(sp)
viviendas_sp <- SpatialPointsDataFrame(coords = datos_sin_na[,12:13],  data = datos_sin_na[, 1:5])
#viviendas_sp=SpatialPointsDataFrame(coords = datos[,12:13],data = datos[,1:5])  "nunca corrio
plot(viviendas_sp)

Carga de archivos shapefiles

require(raster)
#comunas = shapefile()   --> llamando archivo shape

## class       : SpatialPolygonsDataFrame 
## features    : 22 
## extent      : -76.59284, -76.46125, 3.331802, 3.505871  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +no_defs 
## variables   : 8
## names       : casos, comuna,   nombre, zona_recol,         area, perimetro, covid, viajes 
## min values  :     1,      1, Comuna 1,         NA,  2329397.941,  7983.949,     5,    348 
## max values  :    39,     22, Comuna 9,         NA, 12555929.024, 26480.361,    54,   2968

class : indica que es un archivo espacial con geometria de poligonoy dataframe por que intermanete hay una base de datos
features : 22 # indica que son 22 poligonos
extent : Indica la extension
crs : #el CRS indica que está en coordenadas geográficas
variables : 4 # la base interna de la tabla tiene 4 columnas

A menudo utilizan archivos con formato shape, los cuales son comúnmente empleados en sistemas de información geográfica (GIS) para representar datos geoespaciales. Estos archivos constan de un conjunto de archivos que contienen información sobre elementos geográficos, tales como puntos, líneas o polígonos, junto con sus atributos asociados. Para este caso se usará un archivo shape de las comunas de Cali. Por tal motivo se procede a cargar en consola para hacer los análisis pertinentes.

plot(comunas)

Obsérvese como se puede graficar el mapa con la librería ggplot2.

# Convertir el objeto SpatialPolygonsDataFrame a un objeto sf
library(sf)
library(ggplot2)
comunas_sf <- st_as_sf(comunas)

# Graficar el objeto sf con ggplot2
ggplot() +
  geom_sf(data = comunas_sf, fill = "lightblue", color = "black") + 
  labs(title = "Mapa de Comunas en Cali", x = "", y = "") +  # Etiquetas de ejes vacías
  theme_void()  # Usa theme_void() en lugar de theme_minimal()

Para poder ingresar a la tabla de atribtos inmersa en el shapefile, se usa @:

#ingresando a la tabla de atributos
tabla = comunas@data
tabla

##    casos comuna    nombre zona_recol     area perimetro covid viajes
## 1     20      6  Comuna 6       <NA>  5385422 12496.285    19    998
## 2     10      4  Comuna 4       <NA>  4524983 11430.279     8   1493
## 3      8      5  Comuna 5       <NA>  4197624  8441.174    21    664
## 4      7      7  Comuna 7       <NA>  5107267 12547.823     9    734
## 5      3      8  Comuna 8       <NA>  5266743 12178.376    20   1124
## 6      2      9  Comuna 9       <NA>  2899409  7983.949    13    952
## 7      1     21 Comuna 21       <NA>  4828927 16149.224     9    818
## 8     10     13 Comuna 13       <NA>  4737262 10030.224    13   1238
## 9     15      1  Comuna 1       <NA>  3842243 15518.149     5    794
## 10    20      3  Comuna 3       <NA>  3704463 11003.318    18   2121
## 11    18     19 Comuna 19       <NA> 11318091 26480.361    41   2850
## 12    24     12 Comuna 12       <NA>  2329398  8254.095    14    348
## 13    32     10 Comuna 10       <NA>  4297730  9404.276    26   1401
## 14    32     20 Comuna 20       <NA>  2439498  8448.577    14    733
## 15    32     16 Comuna 16       <NA>  4275834 11700.410     9    993
## 16    31     15 Comuna 15       <NA>  4060433 10425.092    10   1171
## 17    30     17 Comuna 17       <NA> 12555929 16841.538    54   2377
## 18    39     18 Comuna 18       <NA>  5428611 11681.301    33   1536
## 19    10     14 Comuna 14       <NA>  4543322  8937.073     7    685
## 20     8     11 Comuna 11       <NA>  3699601  9315.019    20    925
## 21     7      2  Comuna 2       <NA> 11401109 26477.419    39   2968
## 22     6     22 Comuna 22       <NA> 10589125 15552.137    25   1153

Obsérvese como se ingresa a uno de los atributos del objeto shape SpatialPolygonsDataFrame por ejemplo a la comuna número 22.

#Trayendo la comuna 22 solamente
comuna22=comunas[22,]
plot(comuna22)
plot(comuna22,add=TRUE,col="red")

plot(comunas)
plot(comuna22,add=TRUE,col="red")

5. Graficando puntos en shapes.

Shape con puntos espaciales.

viviendas_df <- as.data.frame(viviendas_sp)
library(ggplot2)
# Graficar el shapefile de comunas y los puntos de viviendas
ggplot() +
  geom_polygon(data = comunas, aes(x = long, y = lat, group = group), fill = "lightblue", color = "black") +
  geom_point(data = viviendas_df, aes(x = longitud, y = latitud), color = "red", size = 1) +
  labs(title = "Viviendas en Comunas de Cali", x = "Longitud", y = "Latitud") +
  theme_minimal() +
  coord_sf(xlim = c(min(viviendas_sp$longitud), max(viviendas_sp$longitud)), 
           ylim = c(min(viviendas_sp$latitud), max(viviendas_sp$latitud)))

plot(comunas)
points(viviendas_sp,col="red")

#ingresando a la tabla de atributos
tabla = comunas@data
tabla

##    casos comuna    nombre zona_recol     area perimetro covid viajes
## 1     20      6  Comuna 6       <NA>  5385422 12496.285    19    998
## 2     10      4  Comuna 4       <NA>  4524983 11430.279     8   1493
## 3      8      5  Comuna 5       <NA>  4197624  8441.174    21    664
## 4      7      7  Comuna 7       <NA>  5107267 12547.823     9    734
## 5      3      8  Comuna 8       <NA>  5266743 12178.376    20   1124
## 6      2      9  Comuna 9       <NA>  2899409  7983.949    13    952
## 7      1     21 Comuna 21       <NA>  4828927 16149.224     9    818
## 8     10     13 Comuna 13       <NA>  4737262 10030.224    13   1238
## 9     15      1  Comuna 1       <NA>  3842243 15518.149     5    794
## 10    20      3  Comuna 3       <NA>  3704463 11003.318    18   2121
## 11    18     19 Comuna 19       <NA> 11318091 26480.361    41   2850
## 12    24     12 Comuna 12       <NA>  2329398  8254.095    14    348
## 13    32     10 Comuna 10       <NA>  4297730  9404.276    26   1401
## 14    32     20 Comuna 20       <NA>  2439498  8448.577    14    733
## 15    32     16 Comuna 16       <NA>  4275834 11700.410     9    993
## 16    31     15 Comuna 15       <NA>  4060433 10425.092    10   1171
## 17    30     17 Comuna 17       <NA> 12555929 16841.538    54   2377
## 18    39     18 Comuna 18       <NA>  5428611 11681.301    33   1536
## 19    10     14 Comuna 14       <NA>  4543322  8937.073     7    685
## 20     8     11 Comuna 11       <NA>  3699601  9315.019    20    925
## 21     7      2  Comuna 2       <NA> 11401109 26477.419    39   2968
## 22     6     22 Comuna 22       <NA> 10589125 15552.137    25   1153

5.1 Intersecciones.

crs(viviendas_sp)=crs(comuna22) 
viviendas22=intersect(viviendas_sp,comuna22)  #para que esto pueda pasar, crs = crs
plot(comunas)
points(viviendas22,col="red")

plot(viviendas22)

comuna22=comunas[22,]
plot(comuna22)
points(viviendas22,col="red")

5.2 Intersecciones - Comunas multiples.

#filtrar solo las comunas 20, 21 y 22 podemos usar lo siguiente:
comuna202122=comunas[comunas$comuna>=20,]
plot(comunas)
plot(comuna202122,col="red",add=TRUE)

Filtrando casas estrato 6.

Datos_Vivienda <-data.frame(vivienda)


#Ahora vamos a crear el archivo shapefile con geometria de puntos con base en estos datos (tener en cuenta el archivo no puede tener valores faltantes)


Datos_Vivienda2=na.omit(Datos_Vivienda)
vivienda_map=SpatialPointsDataFrame(coords = Datos_Vivienda2[,12:13],
                                    data = Datos_Vivienda2,
                                    proj4string =crs(comunas))

plot(vivienda_map)

#De manera similar podemos hacer consultas sobre este nuevo archivo shapefile y filtrar solo casas de estrato 6:
pos_casas=which(vivienda_map$tipo == "Casa" & vivienda_map$estrato == 6)

vivienda_casas=vivienda_map[pos_casas,]
plot(vivienda_casas)

Ahora usando Leaflet

#   visalizarlo de manera interactiva utilizando la libreria leaflet y recordar que esta solo funciona 
#con shapes que esten con sistema de referencia de coordenadas geografico:

require(leaflet)
leaflet() %>% addTiles() %>% 
  addCircleMarkers(lng = vivienda_casas$longitud,
                   lat =vivienda_casas$latitud ,
                   radius = 0.1)

require(leaflet)
leaflet() %>% addTiles() %>% 
  addCircleMarkers(lng = vivienda_casas$longitud,
                   lat =vivienda_casas$latitud ,
                   radius = 0.1,clusterOptions = markerClusterOptions() )

5.3 ¿Como poder integrar datos a nivel espacial?

crs(comunas) = "+proj=longlat +ellps=WGS84 +datum=WGS84"

comunas2=spTransform(comunas,CRS("+proj=longlat +ellps=WGS84 +datum=WGS84"))

vivienda_map2=vivienda_map
crs(vivienda_map2)=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84")

total_ofertas=over(x =comunas2 ,y = vivienda_map2[,1],fn =length)
comunas2$ofertas=total_ofertas[,1]
head(comunas2@data)

##   casos comuna   nombre zona_recol    area perimetro covid viajes ofertas
## 1    20      6 Comuna 6       <NA> 5385422 12496.285    19    998      31
## 2    10      4 Comuna 4       <NA> 4524983 11430.279     8   1493      77
## 3     8      5 Comuna 5       <NA> 4197624  8441.174    21    664     104
## 4     7      7 Comuna 7       <NA> 5107267 12547.823     9    734      26
## 5     3      8 Comuna 8       <NA> 5266743 12178.376    20   1124      48
## 6     2      9 Comuna 9       <NA> 2899409  7983.949    13    952      91

spplot(comunas2[,5])

Datos Espaciales - Visualizaciones

Oscar Julian Layton

2024-04-19