Reading, filtering and visualizing geospatial vector data in R

deptos <-  read_sf("./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
CRS:            4326

ggplot() + geom_sf(data = deptos) 

# The CRS 3978 is used in Canada
ggplot() + geom_sf(data = deptos) + coord_sf(crs=st_crs(3978))

ggplot() + geom_sf(data = deptos_utm)

ggplot() + geom_sf(data = valledelcauca) 

munic <-  read_sf("./COL_adm2.shp")
mun_valledelcauca <- munic %>% filter(NAME_1 == "Valle del Cauca")
ggplot() + geom_sf(data = mun_valledelcauca) 

mun_valledelcauca

Simple feature collection with 42 features and 11 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -81.61708 ymin: 2.931806 xmax: -75.7098 ymax: 4.975899
CRS:            4326

ggplot(valledelcauca) +
    geom_sf() +
    geom_sf(data = valledelcauca_points, fill = "antiquewhite") + 
    geom_text(data = valledelcauca_points, aes(x=X, y=Y,label = ID_2), size = 2) +
    coord_sf(xlim = c(-78, -75), ylim = c(3, 5), expand = FALSE)

ggplot(valledelcauca) + 
  geom_sf(data=valledelcauca_points, aes(x=X, y=Y, fill =
                                       ID_2), color = "black", size = 0.25) +
  geom_text(data = valledelcauca_points, 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 = 5))+
  labs(title="Another  Map of Valle del Cauca")

Ignoring unknown aesthetics: x, y

mun_valledelcauca$km2

Units: [m^2]
 [1]   90.47328  217.03732  317.63325  193.84200  695.82052
 [6] 6789.96158  371.49077  131.56671  829.24901  260.71860
[11]  354.94225  797.82338  270.85422  145.15438  506.81757
[16]  327.78503  434.96694  232.89042  188.21736  975.69844
[21]  683.32275  205.68742  123.93291  309.70029  149.54531
[26] 1036.86196  427.67754  136.32020  238.46109  274.89575
[31]  300.76931  503.33924  920.25180  192.05739  241.24827
[36]  608.83412   65.84478  448.32781  170.36758  318.31794
[41]  158.85390  354.38819

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


labels <- mun_valledelcauca$NAME_2

labels

 [1] "Alcalá"              "Andalucía"          
 [3] "Ansermanuevo"        "Argelia"            
 [5] "Bolívar"             "Buenaventura"       
 [7] "Bugalagrande"        "Caicedonia"         
 [9] "Calima"              "Candelaria"         
[11] "Cartago"             "Dagua"              
[13] "El Águila"           "El Cairo"           
[15] "El Cerrito"          "El Dovio"           
[17] "Florida"             "Ginebra"            
[19] "Guacarí"             "Guadalajara de Buga"
[21] "Jamundí"             "La Cumbre"          
[23] "La Unión de Sucre"   "La Victoria"        
[25] "Obando"              "Palmira"            
[27] "Pradera"             "Restrepo"           
[29] "Riofrío"             "Roldanillo"         
[31] "San Pedro"           "Santiago de Cali"   
[33] "Sevilla"             "Toro"               
[35] "Trujillo"            "Tuluá"              
[37] "Ulloa"               "Versalles"          
[39] "Vijes"               "Yotoco"             
[41] "Yumbo"               "Zarzal"             

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

m

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