Chapter 5 Advanced Aesthetics

In this chapter, we will continue to use default dataset MafiaNodes and MafiaEdges in SSNtools to illustrate examples.

This chapter covers the following topics:

  • How to export tmap objects as static images
  • How to create small multiples with tmap
  • How to create an interactive map and publish it online
  • How to add a basemap to static tmap object in plot mode
  • How to add inset map to static tmap object in plot mode
  • How to apply edge bundling to non-planar networks in tmap

You may need one or more of these libraries to complete the codes in this chapter.

library(sf) #required for any functions start with`st_`
library(tidyverse) #required for any process that involves `%>%` 
library(tmap) #required for any map visualization 
library(SSNtools) #required for any steps that involved MafiaNodes or MafiaEdges dataset 
library(igraph) #required for creating network from graph_from_data_frame and calculating network metrics 
library(basemaps) #required for loading basemaps 
library(tigris) #required for downloading boundary shapefiles 
library(grid) #required for inset map
library(stplanr) #optional for od2line function used to convert points to lines
library(edgebundle) #required for edge bundling

5.1 Export Maps

tmap allows users to export tmap object as a static map in the format of pdf, eps, svg, wmf, png, jpg, bmp, or tiff. You can define the resolution, height, and width of the export in the tmap_save arguments.

tmap_save(tm = TMAP_OBJECT, filename='PATH/FILENAME.png')

With these formats, you can import the tmap maps into Adobe Illustrators or Photoshop to further edit the details. You can also export the shapefiles you created in R to experiment in GIS softwares through sf package.

#MafiaNodes is a built-in dataset from SSNtools
MafiaSpatial = MafiaNodes %>% 
  st_as_sf(coords=c("LonX", "LatY"), crs = 4326) 

st_write(MafiaSpatial, "MafiaSpatial.shp") 

Unfortunately, there isn’t a way to export tmap object to GIS softwares without losing the aesthetics yet (e.g., scale, size, color etc., see updates on this issue here).

Another thing to note is that graphics exported through tmap_save may look different from the Plots panel in RStudio. This is because the size of the graphic elements may be scaled dependent on the attributes you put in tmap_save.

5.2 Small Multiples

In our case of the Mafia connectivity map, the power of the Mafia members can be represented by different network metrics. For example, degree centrality (i.e. degree connections) measures power as the number of connections to a Mafia member. Eigenvector centrality (i.e., PageRank centrality) measures power by not only the count of connections, but also the quality of connections so that connecting to other powerful members increases the weight of the link. Betweeness centrality measures power by the ability of a Mafia member to be the “middle-man/middle-woman” and connect groups of other members that would not be connected otherwise. (Local) Clustering coefficient measures whether a Mafia member’s friends tend to connect with each others (i.e., whether your friends are also your friends’ friends).

Noted that closeness centrality is also a common measure of power, but the metric is not well-defined for disconnected graphs. In our dataset, there exists small Mafia cliques that do not interact with others. Thus, we did not use this metric. To contrast these different notions of power in the Mafia map, small multiple comes in handy.

Small multiple is a group of charts or graphs that use the same scale or axes for easy comparison. Here we show how to create a small multiple of the Mafia member map with various network metrics through tmap tm_facets.

Before we map the small multiple, there are some pre-processing. Firstly, when we add network metrics into the MafiaSpatial, each metric is in one column. However, the tm_facet only create multiples based on a column of value. Thus, we need to convert MafiaSpatial into a long table using gather function in tidyverse package (more accuratly, the tidyr package in tidyverse) so that column metric stores the type of metrics (e.g., “degree”, “eigen”, “betweenness”, “cluster”), and the column value stores the metric value.

Secondly, we convert the metric column into a factor column, so that the order of the variable to be mapped is specified. Otherwise, the order of mapping for the small multiple will follow the alphebatical order of the value string in the metric column.

Lastly, we also retrieve the state boundary shapefile from tigris as the background for the map.

library(sf) 
library(tidyverse) 
library(tmap) 
library(SSNtools) 
library(igraph)
library(tigris)

#MafiaNodes and MafiaEdges are built-in dataset from SSNtools 
MafiaSpatial = MafiaNodes %>% 
  st_as_sf(coords=c("LonX", "LatY"), crs = 4326)

g = graph_from_data_frame(MafiaEdges, directed = FALSE, vertices=MafiaSpatial)

#Calculate network metrics
MafiaSpatial$degree = degree(g)
MafiaSpatial$eigen = eigen_centrality(g)$vector
MafiaSpatial$betweenness = betweenness(g, v=V(g))
MafiaSpatial$cluster = transitivity(g, type='local')

MafiaSpatial = MafiaSpatial %>% 
  gather(key='metric', value='value', degree:cluster) %>% 
  mutate(metric = factor(metric, levels=c('degree', 'eigen', 'betweenness', 'cluster')))

#states is a function from tigris package
us_states = states(cb=TRUE, progress_bar = FALSE) %>%
  filter(!STUSPS %in% c('PR','AS', 'AK', 'GU','MP','VI', 'HI'))

Now we can map the small multiple! We specify the column metric in the by argument in tm_facets, which stores variables that the split will use. free.scales=T enables each facet (subplot) creates its own legend with its own scale. free.coords = F enables each map to have its own coordinate range, which is necessary if you have a multi-layer small multiple. You can specify the style of each facet in tm_layout with panel related arguments. In our example, we specify the name for each facet.

tm_shape(us_states) +  
  tm_polygons(alpha=0) + 
  tm_shape(MafiaSpatial) +
  tm_facets(by='metric', free.coords=F, free.scales = T) + 
  tm_symbols(size="value", scale=1, 
             col='lightblue', border.col='#6698CC', 
             title.size=c('Measure of Power')) +
  tm_layout(panel.labels = c('Degree Centrality', 'Eigenvector Centrality',
                             'Betweenness Centrality', 'Clustering Coefficient'))

From the maps above we can see that the most powerful mafia members are located around the New York City region. They are relatively more “powerful” in eigenvector centrality and betweenness centrality measures, and less so in degree centrality and clustering coefficients. Despite the concentration of power is in the Northeast, local mafia members can still be powerful as being well-embedded in a local social network. Many of the mafia families in the isolated regions (e.g., Colorado, Midwest, Texas) are very closely connected with each other.

If you found the different value range across the metrics confusing or misleading (some metrics tend to follow power law), you can also standardize the metric values into respective percentile by adding the following lines before the gather. For example, a metric with value c(4, 1, 2, 3, 5) will be converted to c(0.75, 0, 0.25, 0.5, 1).

MafiaSpatial = MafiaNodes %>% 
  st_as_sf(coords=c("LonX", "LatY"), crs = 4326)

g = graph_from_data_frame(MafiaEdges, directed = FALSE, vertices=MafiaSpatial)

#Convert the metrics to percentile 
MafiaSpatial$degree = rank(degree(g))/length(degree(g))
MafiaSpatial$eigen = rank(eigen_centrality(g)$vector)/length(eigen_centrality(g)$vector)
MafiaSpatial$betweenness = rank(betweenness(g, v=V(g)))/length(betweenness(g, v=V(g)))
MafiaSpatial$cluster = rank(transitivity(g, type='local'))/length(transitivity(g, type='local'))

#same as above 
MafiaSpatial = MafiaSpatial %>% 
  gather(key='metric', value='value', degree:cluster) %>% 
  mutate(metric = factor(metric, levels=c('degree', 'eigen', 'betweenness', 'cluster'))) 

tm_shape(us_states) +  
  tm_polygons(alpha=0) + 
  tm_shape(MafiaSpatial) +
  tm_facets(by='metric', free.coords=F, free.scales = T) + 
  tm_symbols(size="value", scale=1, 
             col='lightblue', border.col='#6698CC', 
             title.size=c('Measure of Power')) +
  tm_layout(panel.labels = c('Degree Centrality', 'Eigenvector Centrality',
                             'Betweenness Centrality', 'Clustering Coefficient'))

5.3 Interactive Maps

tmap allows users to view their maps in interactive modes with pre-defined Leaflet basemap styles, such as OpenStreetMap, Esri.WorldGrayCanvas, and Esri.WorldTopoMap. To view tmap object in interactive mode, set tmap_mode('view'). You also can change basemap styles through the layer widget in the RStudio Plots panel, or defined through tm_basemap() argument. This is a quick and easy way to generate preliminary maps with base maps for presentations and illustrations. The interactive mode can also help you examine your data in a spatial context by adding pop up labels or hover labels.

MafiaSpatial = MafiaNodes %>% 
  st_as_sf(coords=c("LonX", "LatY"), crs = 4326) 
g = graph_from_data_frame(MafiaEdges, directed = FALSE, vertices=MafiaSpatial)
MafiaSpatial$degree = degree(g)

tmap_mode('view')
tm_shape(MafiaSpatial) + 
  tm_symbols(size="degree", scale=2,  
             col='lightblue', border.col='#6698CC', 
             title.size=c('Degree of Connections')) + 
  tm_basemap(server = "OpenTopoMap")

You can also publish your tmap interactive map to web through RPub under your account. As an example, copy and paste the following code in a Rmd file. You can create Rmd file from Rstudio File->New File->R Markdown. In R Markdown, Knit the file and the option to Publish should appear. You can select the RPub option and create a free account. See this post for more details. At the moment, the interactive map from tmap did not have the function to adapt the node size with the zoom level. This function exists in Leaflet package with addCircleMarkers() function (more details see this tutorial)

#### ---- Uncomment the commented lines below in R markdown ----- #### 

#---
#output: html_document
#---

#```{r warning=FALSE, message=FALSE, echo=TRUE, eval=FALSE}

library(tidyverse)
library(sf)
library(tmap)
library(SSNtools)
library(igraph)

MafiaSpatial = MafiaNodes %>% 
  st_as_sf(coords=c("LonX", "LatY"), crs = 4326) 
g = graph_from_data_frame(MafiaEdges, directed = FALSE, vertices=MafiaSpatial)
MafiaSpatial$degree = degree(g)

#set tmap_mode to 'view' will enable interactive mode in Rstudio Plots window
tmap_mode('view')
tm_shape(MafiaSpatial) +
  tm_symbols(size="degree", scale=5, 
             col='lightblue', border.col='#6698CC', 
             title.size=c('Degree of Connections'), 
             #determine the hover label
             id = 'NiceLabel',  
             #determine the pop up contents
             popup.vars = c('Name: '='NiceLabel', 'Mafia Family: '='Family', 
                            'Degree of Connections: '='degree')) +
  #add a function where you can see coordinates on mouseover
  tm_mouse_coordinates()

#```

This is an example publish to RPub and what the output will look like. Click and hover your mouse in the map to explore the pop up labels and see the changing coordinates on the top left corner. You can also switch the base map style in the box below the zoom functions.

You may notice in the tmap interactive map above that when you zoom in, the node size did not change, which makes it difficult to exam areas with clustered nodes (e.g., New York City). For more advanced interactive map functions, we suggest users swtich to Leaflet, which is an open-source library for web-based interactive maps. The following codes generate an interactive map in R and allows users to see the name of the Mafia members when hover over the nodes. The node size dynamically adjusts when you zoom and out of the map. For ideas to reduce node cluttering in the static map, refer to Inset Map.

library(leaflet)

leaflet(MafiaSpatial) %>% addTiles() %>% 
  addCircleMarkers(radius = ~ sqrt(degree), fillColor='#6698CC', fillOpacity = 0.8, label=~htmlEscape(NiceLabel)) 

5.4 Base Map

While interactive map comes with base maps, you can manually add base map into static maps in tmap. This section shows you how to use basemaps R package to add an Carto base map to tmap map in the plot mode.

We use st_bbox to find the bounding box of U.S. states and use st_as_sfc() to transform it into a sfc object for retrieving stars object. Then we use the basemap_stars function in basemaps package to retrieve Mapbox base map based on the input bounding box information. basemaps support basemaps from OpenStreetMaps, Mapbox, Carto, and ESRI. Please refer to basemaps document for more examples on using other basemaps. The base map raster can be visualized by tmap through tm_rgb().

library(tigris)
library(sf)
library(tidyverse)
library(basemaps)
library(tmap)

#states is a function in tigris to download U.S. state boundary shapefile
us_states = states(cb=TRUE, progress_bar = FALSE) %>%
  filter(!STUSPS %in% c('PR','AS', 'AK', 'GU','MP','VI', 'HI')) 

MafiaSpatial = MafiaNodes %>% 
  st_as_sf(coords=c("LonX", "LatY"), crs = 4326) 
g = graph_from_data_frame(MafiaEdges, directed = FALSE, vertices=MafiaSpatial)
MafiaSpatial$degree = degree(g)

#Decide which basemap service to use and what map styles. 
set_defaults(map_service = "carto", map_type = "light")
bgbox = st_bbox(us_states, crs=4326) %>% st_as_sfc()
bg = basemap_stars(bgbox) #process basemap into a raster
## Loading basemap 'light' from map service 'carto'...
tmap_mode('plot')
map = tm_shape(bg) + #add base map
  tm_rgb() + 
  tm_shape(us_states) + 
  tm_polygons(alpha=0) + 
  tm_shape(MafiaSpatial) +
  tm_symbols(size="degree", scale=1, 
             col='lightblue', border.col='#6698CC', 
             title.size=c('Degree of Connections')) + 
  tm_layout(legend.position = c('left', 'bottom'))

map

Here are all the available map_service and map_type:

get_maptypes()
## $osm
##  [1] "streets"      "streets_de"   "streets_fr"   "humanitarian" "topographic" 
##  [6] "hike"         "hillshade"    "grayscale"    "no_labels"    "mtb"         
## 
## $osm_stamen
## [1] "toner"      "toner_bg"   "toner_lite" "terrain"    "terrain_bg"
## [6] "watercolor"
## 
## $osm_thunderforest
##  [1] "cycle"          "transport"      "landscape"      "outdoors"      
##  [5] "transport_dark" "spinal"         "pioneer"        "mobile_atlas"  
##  [9] "neighbourhood"  "atlas"         
## 
## $carto
##  [1] "light"                "light_no_labels"      "light_only_labels"   
##  [4] "dark"                 "dark_no_labels"       "dark_only_labels"    
##  [7] "voyager"              "voyager_no_labels"    "voyager_only_labels" 
## [10] "voyager_labels_under"
## 
## $mapbox
## [1] "streets"   "outdoors"  "light"     "dark"      "satellite" "hybrid"   
## [7] "terrain"  
## 
## $esri
##  [1] "natgeo_world_map"                     
##  [2] "usa_topo_maps"                        
##  [3] "world_imagery"                        
##  [4] "world_physical_map"                   
##  [5] "world_shaded_relief"                  
##  [6] "world_street_map"                     
##  [7] "world_terrain_base"                   
##  [8] "world_topo_map"                       
##  [9] "world_dark_gray_base"                 
## [10] "world_dark_gray_reference"            
## [11] "world_light_gray_base"                
## [12] "world_light_gray_reference"           
## [13] "world_hillshade_dark"                 
## [14] "world_hillshade"                      
## [15] "world_ocean_base"                     
## [16] "world_ocean_reference"                
## [17] "antarctic_imagery"                    
## [18] "arctic_imagery"                       
## [19] "arctic_ocean_base"                    
## [20] "arctic_ocean_reference"               
## [21] "world_boundaries_and_places_alternate"
## [22] "world_boundaries_and_places"          
## [23] "world_reference_overlay"              
## [24] "world_transportation"                 
## [25] "delorme_world_base_map"               
## [26] "world_navigation_charts"

Sometimes you may want a customized bounding box that is not based on data layer. A helpful tool to find the coordinates around your data is to visualize your data in the interactive mode with tm_mouse_coordinates() (see example above). You can also specify the bounding box of your base map using the following code:

#remember to add crs argument for specialized coordinates
box = st_bbox(c(xmin=-77.25586, xmax=41.82353, ymin=-71.10352, ymax=38.94553), crs=4326) %>% st_as_sfc()

5.5 Inset Map

From the current mafia member map we can see that a lot of the mafia members cluster around the New York City region. We are not able to see the details in that areas due to cluttered nodes. This is a common problem at mapping spatial social networks when you have a both regional and local social connections.

To decide on the Inset Map bounding box, you can either use the measurement from other shapefiles (e.g., NYC five boroughs) or use tm_mouse_coordinates() (see example in Interactive Map) to identify ideal bounding box. The example below uses information picked up by the tm_mouse_coordinates().

First, we need to create the inset map. We pick a bounding box around the central New York City including where the most connected mafia member (in degree of connections) locates. We also calculate aspect_ratio for the inset map so that when plugged into the background map, we can preserve the inset map’s natural ratio of height and width.

set_defaults(map_service = "carto", map_type = "light")
# create the bounding box for the inset map and retrieve corresponding base map 
bgbox = st_bbox(c(xmin=-74.065141, xmax=-73.537727, ymin=40.558198, ymax=40.882353), crs=4326) 
inset_bg = basemap_stars(bgbox)

# calculate aspect ratio.  
aspect_ratio = unname((bgbox$ymax - bgbox$ymin)/(bgbox$xmax - bgbox$xmin))

# putting base map as the first layer restricts the bounding box of the map, so no additional filtering for MafiaSpatial is needed.
tmap_mode('plot')
inset = tm_shape(inset_bg) +
  tm_rgb() + 
  tm_shape(MafiaSpatial) + 
  tm_symbols(size="degree", scale=2, 
             col='lightblue', border.col='#6698CC', 
             title.size=c('Degree of Connections'),
             legend.size.show = FALSE) + 
  tm_layout(main.title='Inset of New York',
            main.title.fontface = 2,
            main.title.position = c('center', 'top'))

Then, we use tmap_save to export the background map map created in the Base Map section with the inset map. viewport function comes from the grid package and can be used to position the inset map on the background map. We set an arbitrary value for width and adjust the height according to the aspect ratio of the inset map to preserve its shape. You also can create multiple inset maps and pass to the insets_tm argument as a list (e.g., insets_tm = list(inset1, inset2)). Same for insets_vp. You can pass a list of viewport (e.g., insets_vp = list(viewport(X), viewport(Y))).

#library(grid)
tmap_save(map, insets_tm=inset, 
          insets_vp = viewport(0.17, 0.75, width = 1, height = aspect_ratio*0.75),
          filename='PATH', dpi=600)

We can also add a small black bounding box on the background map to indicate where the inset map locates. We also increase the line width of the frame of the inset map to visually match with the aesthetic of the bounding box.

inset = inset + tm_layout(frame=TRUE, frame.lwd = 2)

map = map + tm_shape(st_bbox(inset_bg) %>% st_as_sfc()) + 
  tm_polygons(alpha=0, border.col = 'black', lwd=2)

tmap_save(map, insets_tm=inset, 
          insets_vp = viewport(0.17, 0.75, width = 1, height = aspect_ratio*0.75),
          filename='PATH', dpi=600)

5.6 Edge Bundling

Edge bundling is a technique to reduce visual cluttering of the edges. The algorithm behind edge bundling essentially bundles edges that go in the same direction and split them apart when they are close to their destinations. The algorithm of a complete edge bundling technique has three steps: first, it clusters edges based on their origin and destination; second, it calculates how to bend the geometry of a straight line given the user input of how “strong” you want the bundle to be. A stronger bundle means that you want more edges and edges far from the bundle to be bundled together; third, the algorithm may summarise the weights of the bundled edges so that the weight of the resultant bundle is larger (and thus visually thicker). You can implement edge bundling in QGIS2 or D3 or Python.

This technique may help reveal structures in your edges when there are many lines, especially if you believe your edges follow a hierarchical structure, yet it is not a savior to your flow map. You should use this technique with caution because it has the following drawbacks:

  • It takes a long time to do edge bundling. We recommend running our codes on a network with 1000-3000 edges to get quick feedback on the visual output.
  • It may bend the edges in unnatural ways.
  • It may be difficult to tell exactly where the edge come from and go to after you bend the edges.
  • It may be difficult to visualize directionality in the flow map.
  • You have to manually tune the parameters (thres in our case) to find the best visual outcome. There is no shortcut to find the optimized parameters at once.
  • Some edge bundling package may only be able to bend the lines, but not summarise the edge weights. It means that it is difficult to visualize the edge weight of the bundle.

The edgebundle package we used is still in development. We showed an example of applying the edgebundle package on the mafia map. To reduce the processing time, we filtered mafia connections to those that have more than 50km in euclidean distance, which results in 1668 edges. To filter by distance, we have to convert the coordinate system to North America Equidistant Conic (crs=102010), use st_distance to compute the distance between two point geometry, and lastly convert the meter object from st_distance into numeric numbers for filtering. It is important to note that you set by_element=TRUE in st_distance so that the output of the function is a vector rather than a matrix. The unit of the distance also follows the coordinate projection system you use.

This is a mafia connectivity map without the edge bundling.

#MafiaEdges and MafiaNodes is a built-in dataset of SSNtools package 

MafiaSpatial = MafiaNodes %>% 
  st_as_sf(coords=c("LonX", "LatY"), crs = 4326)

#Filter MafiaEdges to edges that have more than 50km. 
FilteredEdges = MafiaEdges %>% 
  left_join(MafiaSpatial %>% select(c(NODE, geometry)) %>% st_transform(crs=102010), by = c('Source' = 'NODE')) %>%
  left_join(MafiaSpatial %>% select(c(NODE, geometry)) %>% st_transform(crs=102010), by = c('Target' = 'NODE')) %>%
  mutate(distance = st_distance(geometry.x, geometry.y, by_element=TRUE)) %>% 
  select(-c(geometry.x, geometry.y)) %>%
  mutate(distance = as.numeric(distance)) %>% 
  filter(distance > 50000)

g = graph_from_data_frame(MafiaEdges, directed = FALSE)
MafiaSpatial$degree = degree(g)

#Convert FilteredEdges dataframe to line geometry. 
cnty_centroid = GA_county %>% st_centroid() %>% select(c(GEOID, geometry))
edges = od2line(ctc_commutes, cnty_centroid)

EdgeSpatial = od2line(FilteredEdges, MafiaSpatial %>% select(c(NODE, geometry)))

tmap_mode('plot')
map = tm_shape(bg) + 
  tm_rgb() + 
  tm_shape(us_states) + 
  tm_polygons(alpha=0) + 
  tm_shape(EdgeSpatial) + 
  tm_lines(lwd=0.5, alpha=0.3, col='#1B665D') + 
  tm_shape(MafiaSpatial) +
  tm_symbols(size="degree", scale=2, 
             col='lightblue', border.col='#6698CC', 
             title.size=c('Degree of Connections')) + 
  tm_layout(legend.position = c('left', 'bottom'))

map

The map shows clearly that many of the connections are going to New York City. There are also a lot of connections that go stragith from west coast to east coast. Let’s try edge bundling and see if it makes the map better.

The edge_bundle_force function from the edgebundle pacakage takes into three arguments: 1) g which is the network constructed through igraph, 2) xy which is the longitude and latitude of all the nodes in the network, and 3) compatibility_threshold which indicates the strength of the bundling. A higher value means that the bending and the bundling will be less intense. This package can only bend the edges into bundles but not beyond that.

The result of the edge_bundle_force returns four columns. x and y are the coordinates of the points that consist of a bundled edge. group indicates that coordinates in the same group are for one bundled edge and index indicates the order of points for that particular bundled edge. We can see that for group 1, the fbundle result starts with the first name in the Source of FilteredEdges, BLANDA-CHARLES (-104.6270, 38.2476), and after 34 points, it will end with the Target, SMALDONE-EUGENE (-104.9479, 39.7681).

library(edgebundle)
#FilteredEdges %>% slice(1:3)
#               Source             Target   distance
# 1     BLANDA-CHARLES    SMALDONE-EUGENE   97252.56
# 2     DIVARCO-JOSEPH        SICA-JOSEPH 6054967.18
# 3 DEMARTINO-BENJAMIN DEMARTINO-THEODORE  182989.39

#The current EdgeSpatial is line geometry 
bundle_g = graph_from_data_frame(FilteredEdges, directed=FALSE)
node = data.frame(id = V(bundle_g)$name) %>% mutate(id = as.character(id))
#>                   id
#> 1     BLANDA-CHARLES
#> 2     DIVARCO-JOSEPH
#> 3 DEMARTINO-BENJAMIN

node = node %>% 
  left_join(MafiaNodes %>% select(c(NODE, LonX, LatY)), by=c('id' = 'NODE'), copy=FALSE)

#>                   id      LonX    LatY
#> 1     BLANDA-CHARLES -104.6270 38.2476
#> 2     DIVARCO-JOSEPH  -87.7354 42.0154
#> 3 DEMARTINO-BENJAMIN  -72.9263 40.9530

xy = cbind(node$LonX, node$LatY)
fbundle = edge_bundle_force(bundle_g,xy,compatibility_threshold = 0.85)

#>           x        y      index group
#> 1 -104.6270 38.24760 0.00000000     1
#> 2 -104.6506 38.30726 0.03030303     1
#> 3 -104.6641 38.33735 0.06060606     1

Thus, to create lines acceptable to tmap from the fbundle result, we need to convert coordinates to points, group by the group column, and summarise the points into lines. It is important that do_union=FALSE argument is set in the summarise() to preserve the order of points (otherwise summarise will rearrange the point order).

fbundle2 = fbundle %>%
  st_as_sf(coords=c('x', 'y'), crs=4326) %>%
  group_by(group) %>%
  summarise(do_union=FALSE) %>% st_cast("LINESTRING")

#>   group geometry
#>   <dbl> <LINESTRING [°]>
#> 1     1 (-104.627 38.2476, -104.6506 38.30726, -104.6641 38.33735...)
#> 2     2 (-87.7354 42.0154, -88.8135 41.62707, -89.64513 41.35496 ...)
#> 3     3 (-72.9263 40.953, -72.95451 40.95169, -72.98361 40.94916 ...)

Now the fbundle2 is a line shapefile! We can either export to GIS for further editing or map it using tmap package. See Base Map on how to create bg background base map and us_states boundary shapefiles.

tmap_mode('plot')
tm_shape(bg) + 
  tm_rgb() + 
  tm_shape(us_states) + 
  tm_polygons(alpha=0) + 
  tm_shape(fbundle2) + 
  tm_lines(lwd=0.5, alpha=0.3, col='#1B665D') + 
  tm_shape(MafiaSpatial %>% arrange(desc(degree))) +
  tm_symbols(size="degree", scale=2, 
             col='lightblue', border.col='#6698CC', 
             title.size=c('Degree of Connections')) + 
  tm_layout(legend.position = c('left', 'bottom'))

The bundling effect is most obvious when you contrast lines going from Florida to New York. Several edges have been bundled together. However, you can also see some unnatural bends from lines going from Florida to New York City. That is because the lines were bended halfway. In our case of the mafia map, we are not sure the edge bundling improves the visual. Again, edge bundling should be applied with cautions and experiments.