Santander

library(tidyverse)

-- Attaching packages --------------------------------------- tidyverse 1.3.0 --
v ggplot2 3.3.0     v purrr   0.3.3
v tibble  2.1.3     v dplyr   0.8.4
v tidyr   1.0.2     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)

Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3

lo primero a realizar es abrir el archivo shape file que nos servira como base de datos geograficos, para ello usamos la funcion read_sf y almacenar estos datos en una variable

deptos<-read_sf("./santander/COL_adm/COL_adm1.shp")
head(deptos)

Simple feature collection with 6 features and 9 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -77.149 ymin: -4.228429 xmax: -69.36835 ymax: 11.10792
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs

##luego usamos la funcion ggplot para graficar los datos almacenados en la variable deptos

ggplot()+geom_sf(data=deptos)

esta linea muestra el resultado del ploteo con el nuevo sistema de coordenadas, el cual es usado en Canada

ggplot() + geom_sf(data = deptos) + coord_sf(crs=st_crs(3978))

esta linea muestra un cambio de sistema de coordenadas de la variable deptos, ya que ahora se trabajara con el sistema UTM 18N

deptos.UTM<-st_transform(deptos,crs=st_crs(32618))
deptos.UTM

Simple feature collection with 32 features and 9 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -245935.3 ymin: -469204.3 xmax: 1407491 ymax: 1763314
epsg (SRID):    32618
proj4string:    +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs

esta linea muestra un grafico de Colombia bajo estas nuevas coordenadas

ggplot() + geom_sf(data = deptos.UTM)

aqui se realiza un ploteo del departamento a estudiar, en este caso Santander

Santander <-  deptos %>%   filter(NAME_1 == "Santander")
ggplot()+geom_sf(data=Santander)

estas lineas muestran como obtener no solo el grafico dela region del departamento, sino que nos muestra la division administativa del mismo

municipios <-  read_sf("./santander/COL_adm/COL_adm2.shp")
mun_santander <- municipios %>% filter(NAME_1 == "Santander")
ggplot() + geom_sf(data = mun_santander) 

##aqui vemos los datos almacenados en esta variable

mun_santander

Simple feature collection with 87 features and 11 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -74.55179 ymin: 5.705701 xmax: -72.5161 ymax: 8.120555
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs

Esta instruccion es previa para un nuevo ploteo, en la cual se le dira donde debe colocar los textos aqui asignados

puntos_santander <- cbind(mun_santander, st_coordinates(st_centroid(mun_santander$geometry)))

Etsta instruccion nos muestra un ploteo del departamento y sus municipios, solo que esta vez los municipios llevan el numero identificador de la base de datos ubicado segun la linea anterior

ggplot(Santander) +
    geom_sf() +
    geom_sf(data = puntos_santander, fill = "antiquewhite") + 
    geom_text(data = puntos_santander, aes(x=X, y=Y,label = ID_2), size = 2) +
    coord_sf(xlim = c(-74.7, -72.4), ylim = c(5.6, 8.15), expand = FALSE)+labs(title="Santander")

##esta nueva libreria, servira para darle color al mapa en base a una escala especifica

library(scales)

##este ploteo nos muestra una variante donde los diferentes municipios estan identificados con su codigo y pintados de colores en base a una escala determinada

ggplot(Santander) + 
  geom_sf(data=puntos_santander, aes(x=X, y=Y, fill =                           ID_2), color = "black", size = 0.3) +
  geom_text(data = puntos_santander, aes(x=X, y=Y,label = ID_2), size = 2) +
  theme(aspect.ratio=1)+
  scale_fill_distiller(name="ID_2", palette = "YlGn", breaks = pretty_breaks(n = 10))+
  labs(title="otro estilo de mapa para Santander")

Ignoring unknown aesthetics: x, y

con esto se guarda el mapa antes realizado en formato pdf y png

ggsave("Santander_municipios.pdf")
ggsave("map_santander.png", width = 6, height = 6, dpi = "screen")

esta libreria sera util para realizar un nuevo mapa

install.packages("leaflet")
library(leaflet)

para usar esta libreria es necesario que los datos sean de tipo espacial, para ello usamos esta linea

ant_points <- as(puntos_santander, 'Spatial')

con esto veremos que hay almacenado en esta variable

head(ant_points)

esta libreria sera util para conocer algunos valores de cada dato

install.packages("lwgeom")

library(lwgeom)

con esto podemos hallar el area de cada uno de los datos

mun_santander$area <- st_area(mun_santander)
mun_santander$km2 <- mun_santander$area/(1000000)
mun_santander$km2

para realizar el mapa, es necesario que la variable mun_santander sea de tipo espacial, para ello, usamos esta instruccion

sant_municipios <- as(mun_santander, 'Spatial')
head(sant_municipios)

esta instruccion es util para el proximo ploteo, servira para establecer la escala de colores

bins <- c(0, 50, 100, 200, 300, 500, 1000, 2000, Inf)
pal <- colorBin("YlOrRd", domain = ant_municipios$km2, bins = bins)


labels <- mun_santander$NAME_2

labels

 [1] "Aguada"                 "Albania"               
 [3] "Aratoca"                "Barbosa"               
 [5] "Barichara"              "Barrancabermeja"       
 [7] "Betulia"                "Bolívar"               
 [9] "Bucaramanga"            "Cabrera"               
[11] "California"             "Capitanejo"            
[13] "Carcasí"                "Cepitá"                
[15] "Cerrito"                "Charalá"               
[17] "Charta"                 "Chimá"                 
[19] "Chipatá"                "Cimitarra"             
[21] "Concepción"             "Confines"              
[23] "Contratación"           "Coromoro"              
[25] "Curití"                 "El Carmen de Chucurí"  
[27] "El Guacamayo"           "El Peñon"              
[29] "El Playón"              "Encino"                
[31] "Enciso"                 "Florián"               
[33] "Floridablanca"          "Gámbita"               
[35] "Güepsa"                 "Galán"                 
[37] "Girón"                  "Guaca"                 
[39] "Guadalupe"              "Guapotá"               
[41] "Guavatá"                "Hato"                  
[43] "Jesús María"            "Jordán"                
[45] "La Belleza"             "La Paz"                
[47] "Landázuri"              "Lebrija"               
[49] "Los Santos"             "Málaga"                
[51] "Macaravita"             "Matanza"               
[53] "Mogotes"                "Molagavita"            
[55] "Ocamonte"               "Oiba"                  
[57] "Onzaga"                 "Páramo"                
[59] "Palmar"                 "Palmas del Socorro"    
[61] "Piedecuesta"            "Pinchote"              
[63] "Puente Nacional"        "Puerto Parra"          
[65] "Puerto Wilches"         "Rionegro"              
[67] "Sabana de Torres"       "San Andrés de Cuerquia"
[69] "San Benito"             "San Gil"               
[71] "San Joaquín"            "San José de Miranda"   
[73] "San Miguel"             "San Vicente de Chucurí"
[75] "Santa Bárbara"          "Santa Helena del Opón" 
[77] "Simacota"               "Socorro"               
[79] "Suaita"                 "Sucre"                 
[81] "Suratá"                 "Tona"                  
[83] "Vélez"                  "Valle de San José"     
[85] "Vetas"                  "Villanueva"            
[87] "Zapatoca"              

este ploteo muestra otra forma de ver el departamento

m <- leaflet(sant_municipios) %>%
  setView(-73.5, 6.5, 7.45)  %>% addPolygons(
  fillColor = ~pal(km2),
  weight = 2,
  opacity = 1,
  color = "white",
  dashArray = "3",
  fillOpacity = 0.7,
  highlight = highlightOptions(
    weight = 5,
    color = "#555",
    dashArray = "",
    fillOpacity = 0.7,
    bringToFront = TRUE),
  label = labels) %>%
  addLegend(pal = pal, values = ~km2, opacity = 0.7, title = NULL,
    position = "bottomright")
m

estas ultimas lineas de codigo muestran las lineas administrativas del mapa sobre una proyeccion del continente, esto para ubicar la zona de trabajo

a<-leaflet() %>%
  addProviderTiles(providers$Esri.WorldImagery, options= providerTileOptions(opacity = 0.99)) %>%
  addPolygons(data = sant_municipios, popup= sant_municipios$NAME_2,
    stroke = TRUE, fillOpacity = 0.1, smoothFactor = 0.25)
    a