Libraries

library(readr)
library(dplyr)
library(leaflet)
library(sf)

About

An interactive map could be a great product to support a spatial analysis. This project shows how to create an interactive map using the leaflet package (library(leaflet)) in R.The leaflet package use data that could be read as spatial data, for example, coordinates to add markers or polygons to delimit an area.

Data

I use data from the National Statistical Directory of Economic Units (DENUE in spanish) of the National Institute of Statistics and Geography (INEGI in spanish), this directory contains information on all the establishments in Mexico.

To download the data you can visit the website of DENUE, the data could be filter by Economic Activity, Size of Establishments and/or Geographic Area.

For this project, I focus on veterinary services in the Mexico City. I use the exact location provided by latitude and longitude, and the name of establishments.

Manage data

Read the data from DENUE. This file contains all the establishments in Mexico City.

denue_inegi <- read_csv("denue_inegi_09_.csv", show_col_types = F, locale = locale(encoding = "LATIN1"))
head(denue_inegi[c("nom_estab", "nombre_act","entidad", "municipio", "latitud", "longitud")])
## # A tibble: 6 × 6
##   nom_estab                   nombre_act      entidad municipio latitud longitud
##   <chr>                       <chr>           <chr>   <chr>       <dbl>    <dbl>
## 1 AJOLOTARIO APATLACO         Piscicultura e… Ciudad… Xochimil…    19.3    -99.1
## 2 AJOLOTARIO CARRIZAL         Otra acuicultu… Ciudad… Xochimil…    19.3    -99.1
## 3 AJOLOTARIO CHINAMPA ONKALI  Otra acuicultu… Ciudad… Xochimil…    19.3    -99.1
## 4 AJOLOTARIO TLAZOCAMA TONANA Otra acuicultu… Ciudad… Xochimil…    19.3    -99.1
## 5 AJOLOTERIA APANTLI          Otra acuicultu… Ciudad… Xochimil…    19.3    -99.1
## 6 AMILLI TLAPALLI             Piscicultura e… Ciudad… Xochimil…    19.3    -99.1

Filter by veterinary services. Based on the North American Industry Classification System, Mexico SCIAN 2023 used by INEGI. The codes for these activities are 541941, 541942, 541943 and 541944.

data_vet <- denue_inegi %>% 
  filter(codigo_act >= 541941 & codigo_act <= 541944) 

head(data_vet[c("nom_estab", "nombre_act","entidad", "municipio", "latitud", "longitud")])
## # A tibble: 6 × 6
##   nom_estab                   nombre_act      entidad municipio latitud longitud
##   <chr>                       <chr>           <chr>   <chr>       <dbl>    <dbl>
## 1 A PATA DE PERRO             Servicios vete… Ciudad… Venustia…    19.4    -99.1
## 2 ACUA MASCOTA                Servicios vete… Ciudad… Iztapala…    19.3    -99.0
## 3 ACUAMASCOTAS VETERINARIA    Servicios vete… Ciudad… Azcapotz…    19.5    -99.2
## 4 ADOPTANDO UN CORAZON CANINO Servicios vete… Ciudad… Miguel H…    19.4    -99.2
## 5 ALFA PETS                   Servicios vete… Ciudad… Gustavo …    19.5    -99.1
## 6 AMALTEA                     Servicios vete… Ciudad… Coyoacán     19.3    -99.1

For a map by geographic area, I aggregate the data by municipality.

data_vet$vet <- 1 
data_vet_mun <- data_vet %>%
  group_by(municipio) %>%
  summarise(total_vet = sum(vet))

head(data_vet_mun)
## # A tibble: 6 × 2
##   municipio             total_vet
##   <chr>                     <dbl>
## 1 Azcapotzalco                107
## 2 Benito Juárez               199
## 3 Coyoacán                    170
## 4 Cuajimalpa de Morelos        40
## 5 Cuauhtémoc                  130
## 6 Gustavo A. Madero           286

Polygons

Read the shapefile with the municipality polygons. You can download the polygons in the Data Portal of CDMX.

shp_mun <- st_read("alcaldias/poligonos_alcaldias_cdmx.shp", quiet = TRUE)

Now, we merge the data with the polygons by municipality.

mun_sf <- merge(shp_mun , data_vet_mun, by.x = "NOMGEO", by.y = "municipio", all.x = TRUE)

Get the centroids, this step help with the elements in the map.

centroides <- st_centroid(mun_sf)
coords <- st_coordinates(centroides)
mun_sf$lon <- coords[,1]
mun_sf$lat <- coords[,2]

Finally the data by municipality is the next.

head(mun_sf)
## Simple feature collection with 6 features and 7 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -99.36492 ymin: 19.23257 xmax: -99.09908 ymax: 19.51514
## Geodetic CRS:  WGS 84
##                  NOMGEO CVEGEO CVE_ENT CVE_MUN total_vet
## 1        Álvaro Obregón  09010      09     010       155
## 2          Azcapotzalco  09002      09     002       107
## 3         Benito Juárez  09014      09     014       199
## 4              Coyoacán  09003      09     003       170
## 5 Cuajimalpa de Morelos  09004      09     004        40
## 6            Cuauhtémoc  09015      09     015       130
##                         geometry       lon      lat
## 1 POLYGON ((-99.18906 19.3955... -99.24683 19.33617
## 2 POLYGON ((-99.18231 19.5074... -99.18211 19.48533
## 3 POLYGON ((-99.14762 19.4040... -99.16113 19.38064
## 4 POLYGON ((-99.13427 19.3565... -99.15038 19.32667
## 5 POLYGON ((-99.25738 19.4011... -99.31094 19.32431
## 6 POLYGON ((-99.12951 19.4626... -99.14906 19.43137

Interactive map with Leaflet

I present the total of veterinarian services by municipality. This map shows the total of veterinarian services by municipality in Mexico City. Next, I explain how the code works.

  • leaflet package use %>% as pipe to create a complete map.

  • addProviderTiles(providers$CartoDB.Positron add the basemap

  • addPolygons() add the area of the municipalities and its appearance. popup = ~ NOMGEO is the pop up in the map, it is the municipality label that appears.

  • addCircles() add the circles in the map. This circles are centered in each municipality using the centroids calculated before, and its radius is based in the total of veterinarian in the area. Here, popup = ~paste0(NOMGEO, ": ", total_vet) is used as in addPolygons() but just for the area of the circle, note that this pop up contains more information.

  • addControl() add the title in the top right of the map.

leaflet(mun_sf) %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  addPolygons( color = "darkblue", weight = 2,
    fillColor = "white", fillOpacity = 0.5,
    popup = ~ NOMGEO) %>%
  addCircles(lng = ~lon, lat = ~lat,radius = ~total_vet*10,      
    color = "white",fillColor = "blue",fillOpacity = 0.7,
    popup = ~paste0(NOMGEO, ": ", total_vet) ) %>%
  addControl(html = "<div style='text-align: center; font-size: 16px;'> <b> Total of Veterinarian Services in Mexico City </b> </div>", position = "topright" )

This map shows the distribution of all veterinarian establishments in Mexico city.Next, I explain how the code works, note that some functions are used before.

  • makeIcon() upload a custom icon to use as marker in the map.

  • addMarkers add markers in the map using coordinates (latitude and longitude). Each marker represent a veterinarian establishment.

  • setView() configure the map view with a specific zoom level and coordinates..

icon_gato <- makeIcon(iconUrl = "pets.png",  iconWidth = 15, iconHeight = 15)

leaflet(data = data_vet) %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  addMarkers(lng = ~longitud, lat = ~latitud, icon = icon_gato, label = ~nom_estab) %>%
 addControl(html = "<div style='text-align: center; font-size: 16px;'> <b> Veterinarian Services in Mexico City </b> </div>", position = "topright" ) %>% 
  setView(lng = -99.140487, lat = 19.430072, zoom = 13)