A demo of how to generate a map in R that can be edited in Adobe Illustrator.
This content was presented to Nelson\Nygaard Staff at a GIS Hangout on Wednesday, June 23rd, 2021, and is available as a recording here and embedded below.
A downloadable package of the project files Esther presented is available here, and the main script is also embedded below.
## ---------------------------
##
## Script name:
##
## Purpose of script:
##
## Project Name:
##
## Project Number:
##
## Author: Esther Needham
##
## Date Created: 2021-06-22
##
## ---------------------------
##
## Notes:
##
##
## ---------------------------
## set working directory here or make sure there is a .Rproj file
## ---------------------------
options(scipen = 6, digits = 4) # view outputs in non-scientific notation
## ---------------------------
## load up the packages we will need: (uncomment as required)
require(tidyverse)
require(sf)
require(ggthemes) # preset "themes" for plots
require(classInt) # for creating breaks
require(ggsn) # north arrows and scale bars
require(scales) # for formatting of text and numbers
## ---------------------------
## load up functions from other scripts into memory
# source("functions/summarise_data.R") #example
## ---------------------------
# Read in Data ------------------------------------------------------------
neighborhoods <- st_read('input/neighborhoods_sf.shp')
streets <- st_read('input/streets_sf.shp')
inc <- st_read('input/ACS18_5YR_IncomeRatio_BG.shp')
# Colors ------------------------------------------------------------------
# Colors
text_color <- '#27323d'
cumulative_colors <- c('#ffda09', '#ccd84d', '#9ad27e', '#67cb9b', '#34bfa4', '#01b199')
cumulative_colors_districts <- c('#ffda09', '#9ad27e', '#67cb9b', '#01b199')
risk_colors <- c('#fff3c3', '#ffeea6', '#ffe889', '#ffe36a', '#ffde43', '#ffda09')
risk_colors_districts <- c('#fff3c3', '#ffe889', '#ffe36a', '#ffda09')
asset_colors <- c('#01b199', '#0eb9a4', '#58c3b1', '#7eccbf', '#9fd8cd', '#bee4dc')
asset_colors_districts <- c('#01b199', '#58c3b1', '#7eccbf', '#bee4dc')
# Functions ---------------------------------------------------------------
map_style_legend <- function(){
theme_void() + theme(legend.position = "bottom",
legend.text=element_text(size=10, color = text_color),
legend.title=element_text(size=15, color = text_color),
panel.grid.major = element_blank(),
plot.margin = unit(c(0, 0, 0.5, 0), "cm"))#top, #right #bottom #left
}
map_style_without_legend <- function(){
theme_void() + theme(legend.position = "none",
panel.grid.major = element_line(colour = 'white'))
}
# this needs to be edited to match the crs being used
# review documentation for options https://cran.r-project.org/web/packages/ggsn/ggsn.pdf
scale_bar <- function(data, dist){
scalebar(data, st.bottom = FALSE, location="bottomright", dist = dist, dist_unit = "mi",
transform = FALSE, model = "WGS84", # ellipsoid model
st.color=text_color, box.fill=c('#929090', 'white'),
box.color='#929090', border.size = 0.25, st.size=3.5, height=0.02)
}
fill_color_discrete <- function(colorList, name, labels){
scale_fill_manual(values = c(colorList[1:length(colorList)]),
name = name,
labels = labels,
na.value="#929090",
guide = guide_legend(
default.unit="inch",
keywidth=1,
keyheight=0.2,
direction = "horizontal",
nrow = 1,
byrow = T,
draw.ulim = F,
title.position = 'top',
title.hjust = 0.5,
label.position = 'bottom')
)
}
fill_color_stretch <- function(low_color, high_color, name){
scale_fill_gradient(low = low_color, high = high_color,
name = name,
na.value="#929090",
guide = guide_colorbar(
direction = "horizontal",
barheight = unit(3, units = "mm"),
barwidth = unit(60, units = "mm"),
draw.ulim = F,
title.position = 'top',
title.hjust = 0.5,
label.hjust = 0.5)
)
}
createBreaks <- function(data, variable, num, style){
breaks <- classIntervals(data[[variable]], num, style = style)$brks
return (breaks)
}
createRoundedLabels <- function(data, variable, num, style){
breaks <- createBreaks(data, variable, num, style)
breaks_rounded <- comma(round(breaks,0))
these_labels <- vector()
counter <- 0
for(i in breaks[1:length(breaks)-1]){
counter <- counter + 1
these_labels <- c(these_labels, paste(breaks_rounded[counter], '-',breaks_rounded[counter + 1])) # change to paste instead of paste0 for spaces between the dash
}
if(any(is.na(data[[variable]]))){
these_labels <- c(these_labels, 'No value')
}
return (these_labels)
}
createBins <- function(data, variable, num, style){
breaks <- createBreaks(data, variable, num, style)
data <- data %>%
mutate(bins = findInterval(data[[variable]], breaks)) %>%
mutate(bins = ifelse(bins == max(bins, na.rm=T), bins - 1, bins))
data$bins <- as.factor(data$bins)
return(data)
}
# A note on coordinate systems ----------------------------------
# take a look at default NC dataset
nc = st_read(system.file("shape/nc.shp", package="sf"))
st_crs(nc)
ggplot() +
geom_sf(data = nc)
# you can reproject and set the coord sf to match the crs
# meters
nc_2163 = st_transform(nc, 2163)
ggplot() +
geom_sf(data = nc_2163) #+
# coord_sf(datum = st_crs(2163))
# Maps --------------------------------------------------------------------
# With default ggplot settings
ggplot() +
geom_sf(data = inc, aes(fill = Under200PC), colour = NA) +# colour = NA remove the lines width, change to size = 1 (or whatever number) to get a line width
coord_sf(datum = st_crs(2227))
ggplot() +
geom_sf(data = inc %>% st_transform(4326), aes(fill = Under200PC), colour = NA)
# Removing the background junk
# With theme_void
ggplot() +
geom_sf(data = inc, aes(fill = Under200PC), colour = NA) +
theme_void()
# Testing exporting
# this combines neighborhoods and block groups into one clip group
# setting obviously different line widths or colors makes it easy in illustrator to group them back together as their own layers
map <-
ggplot() +
geom_sf(data = inc, aes(fill = Under200PC), colour = NA) +
geom_sf(data = neighborhoods, size = 1, fill = NA, colour = "white") +
theme_void()
ggsave('output/graphics/block_groups_neighborhoods.svg', plot = map, device = 'svg', width = 8, height = 10)
# You can set the colors for the choropleth
ggplot() +
geom_sf(data = inc, aes(fill = Under200PC), colour = NA) +
scale_fill_gradient(low = 'yellow', high = 'purple')
# and decide where the legend goes
ggplot() +
geom_sf(data = inc, aes(fill = Under200PC), colour = NA) +
scale_fill_gradient(low = 'yellow', high = 'purple') +
theme(legend.position = "bottom") # and move the legend
# Using bins for a choropleth instead of continuous (default)
# you have to tell ggplot what the bins should be
# this makes the data categorical, so can use similar idea for other types of data
# this is a function from above, output is a dataframe with a variable called bins
binned_data <- createBins(inc, 'Under200PC', 6, 'jenks')
ggplot() +
geom_sf(data = binned_data, aes(fill = bins), colour = NA)
# But you can manually set the colors
# Also notice here I'm just using the function to create the bins directly in the ggplot call
# this works because the output is a dataframe, no real difference than naming the output above and feeding that in
ggplot() +
geom_sf(data = createBins(inc, 'Under200PC', 6, 'jenks'), aes(fill = bins), colour = NA) +
scale_fill_manual(values = cumulative_colors)
# Adding pretty labels that state what the bins represent
# createRoundedLabels is a function I wrote, see it above. It takes the same inputs as the createBins function
# finds the breaks and rounds the data.
ggplot() +
geom_sf(data = createBins(inc, 'Under200PC', 6, 'jenks'), aes(fill = bins), colour = NA) +
scale_fill_manual(values = cumulative_colors,
labels = createRoundedLabels(inc, 'Under200PC', 6, 'jenks'))
# Bringing it all together and making it look nice
# These maps use a function call map_style_legend which is created above,
# it uses the theme_void and sets a few other parameters, like moving the legend to the bottom.
# You could just use theme_void if you don't want the other settings, or set them manually
# The scale_fill_manual and scale_fill_gradient have also been wrapped in two functions, created above
# they both set some setting for the legends, like how high or how wide the legend bars are
# these settings are helpful to tweak, especially if you have very large numbers and need more space
# the guide settings can be placed directly in the ggplot call and do not need to be a part of function if you want to set them manually
# Categorical/binned
ggplot() +
geom_sf(data = createBins(inc, 'Under200PC', 6, 'jenks'), aes(fill = bins), colour = NA) +
fill_color_discrete(cumulative_colors, "Low Income Residents", createRoundedLabels(inc, 'Under200PC', 6, 'jenks')) +
# scale_bar(inc, 2) + # we can also add a scalebar, this is a function set above the wraps the ggsn package scalebar function and create pre-set settings
# north(data = inc, symbol = 10) + # setting a north arrow, see ggsn documentation, northSymbols() for style types
map_style_legend()
# Continuous
ggplot() +
geom_sf(data = inc, aes(fill = Under200PC), colour = NA) +
fill_color_stretch(low_color = cumulative_colors[1], high_color = cumulative_colors[length(cumulative_colors)], "Low Income Residents") +
map_style_legend()
# One main issue remains -- the extent of the map
# ggplot defaults to the extent of the data, but you can manually set this with coord_sf()
# coord_sf wants and x limit and a y limit, to create those, use the following piece of code based on
# whatever dataset (or subset) you want the map extent to be
# in this case we'll use the neighborhoods file
# create coordinates for the map range, using the neighborhoods file
mapRange <- c(range(st_coordinates(neighborhoods)[,1]),
range(st_coordinates(neighborhoods)[,2]))
# Setting the map to be centered around the neighborhoods file
# bg <-
ggplot() +
geom_sf(data = createBins(inc, 'Under200PC', 6, 'jenks'), aes(fill = bins), colour = NA) +
coord_sf(xlim = mapRange[c(1:2)], ylim = mapRange[c(3:4)]) +
map_style_legend() +
fill_color_discrete(cumulative_colors, "Low Income Residents", createRoundedLabels(inc, 'Under200PC', 6, 'jenks'))
ggsave('output/graphics/block_groups_pretty.svg', plot = bg, device = 'svg', width = 12, height = 12)
ggsave('output/graphics/block_groups_pretty.eps', plot = bg, device = 'eps', width = 12, height = 12)
ggsave('output/graphics/block_groups_pretty.pdf', plot = bg, device = 'pdf', width = 12, height = 12)
# to create labels, you need an x y
# centroids aren't perfect, but a good approximation
neigh_centroids <- st_centroid(neighborhoods) %>%
mutate(x = st_coordinates(.)[,1],
y = st_coordinates(.)[,2])
bg_w_labs <-
ggplot() +
geom_sf(data = createBins(inc, 'Under200PC', 6, 'jenks'), aes(fill = bins), colour = NA) +
geom_sf(data = neighborhoods, colour = "white", size = 0.5, fill = NA) +
geom_text(data = neigh_centroids, aes(label = name, x = x, y = y), # label all or you could make a selection
size = 4, color = 'white', hjust = 'right', check_overlap = TRUE) +
coord_sf(xlim = mapRange[c(1:2)], ylim = mapRange[c(3:4)]) +
scale_bar(neighborhoods, 1) + # scale bar location on the map is based on the input file, changed to neighborhoods (location of block groups puts the scale bar out of range)
map_style_legend() +
fill_color_discrete(cumulative_colors, "Low Income Residents", createRoundedLabels(inc, 'Under200PC', 6, 'jenks'))
ggsave('output/graphics/block_groups_pretty_labs.svg', plot = bg_w_labs, device = 'svg', width = 12, height = 12)
ggsave('output/graphics/block_groups_pretty_labs.eps', plot = bg_w_labs, device = 'eps', width = 12, height = 12)
ggsave('output/graphics/block_groups_pretty_labs.pdf', plot = bg_w_labs, device = 'pdf', width = 12, height = 12)
# Making a few more maps to demonstrate how you could generate a map set with matching extents
# Castro neighborhood
# create coordinates for the map range, using the neighborhoods file
mapRange_castro <- c(range(st_coordinates(neighborhoods %>% filter(name == 'Castro'))[,1]),
range(st_coordinates(neighborhoods %>% filter(name == 'Castro'))[,2]))
ggplot() +
geom_sf(data = createBins(inc, 'Under200PC', 6, 'jenks'), aes(fill = bins), colour = NA) + # colour = NA remove the lines width, change to size = 1 (or whatever number) to get a line width
geom_sf(data = neighborhoods, size = 1, fill = NA) +
coord_sf(xlim = mapRange_castro[c(1:2)], ylim = mapRange_castro[c(3:4)]) +
scale_bar(neighborhoods, 2) + # scale bar location on the map is based on the input file, changed to neighborhoods (location of block groups puts the scale bar out of range)
map_style_legend() +
fill_color_discrete(cumulative_colors, "Low Income Residents", createRoundedLabels(inc, 'Under200PC', 6, 'jenks'))
# Inner sunset neighborhood
# create coordinates for the map range, using the neighborhoods file
mapRange_inner_sunset <- c(range(st_coordinates(neighborhoods %>% filter(name == 'Inner Sunset'))[,1]),
range(st_coordinates(neighborhoods %>% filter(name == 'Inner Sunset'))[,2]))
ggplot() +
geom_sf(data = createBins(inc, 'Under200PC', 6, 'jenks'), aes(fill = bins), colour = NA) + # colour = NA remove the lines width, change to size = 1 (or whatever number) to get a line width
geom_sf(data = neighborhoods, size = 1, fill = NA) +
coord_sf(xlim = mapRange_inner_sunset[c(1:2)], ylim = mapRange_inner_sunset[c(3:4)]) +
scale_bar(neighborhoods, 2) + # scale bar location on the map is based on the input file, changed to neighborhoods (location of block groups puts the scale bar out of range)
map_style_legend() +
fill_color_discrete(cumulative_colors, "Low Income Residents", createRoundedLabels(inc, 'Under200PC', 6, 'jenks'))
# Notice two things about the neighborhood maps
# 1. The legend does not update
# 2. The shape of the map is based on the shape of the neighborhood and does not stay consistent
# The solution to #1 is to clip or filter the block group data. If anyone knows of another solution, I would love to know
# The solution to #2 is to manually set a consistent map shape with the code below. Again, if anyone knows of another solution, please share
# Notice below that exporting with a set height and width makes the shape of the image the same but not the shape of the map
# Fixing issue # 1
# Filtering block group data based on intersection with selected neighborhood
# there are fewer data points so you may want to reduce the number of bins
# c <-
ggplot() +
geom_sf(data = createBins(inc %>% filter(st_intersects(., neighborhoods %>% filter(name == 'Castro'), sparse = FALSE)), 'Under200PC', 4, 'jenks'), aes(fill = bins), colour = NA) + # colour = NA remove the lines width, change to size = 1 (or whatever number) to get a line width
geom_sf(data = neighborhoods, size = 1, fill = NA) +
coord_sf(xlim = mapRange_castro[c(1:2)], ylim = mapRange_castro[c(3:4)]) +
scale_bar(neighborhoods, 2) + # scale bar location on the map is based on the input file, changed to neighborhoods (location of block groups puts the scale bar out of range)
map_style_legend() +
fill_color_discrete(cumulative_colors, "Low Income Residents",
createRoundedLabels(inc %>% filter(st_intersects(., neighborhoods %>% filter(name == 'Castro'), sparse = FALSE)), 'Under200PC', 4, 'jenks'))
# ggsave('output/graphics/castro.png', plot = c, device = 'png', width = 8, height = 8)
# is <-
ggplot() +
geom_sf(data = createBins(inc %>% filter(st_intersects(., neighborhoods %>% filter(name == 'Inner Sunset'), sparse = FALSE)), 'Under200PC', 4, 'jenks'), aes(fill = bins), colour = NA) + # colour = NA remove the lines width, change to size = 1 (or whatever number) to get a line width
geom_sf(data = neighborhoods, size = 1, fill = NA) +
coord_sf(xlim = mapRange_inner_sunset[c(1:2)], ylim = mapRange_inner_sunset[c(3:4)]) +
scale_bar(neighborhoods, 2) + # scale bar location on the map is based on the input file, changed to neighborhoods (location of block groups puts the scale bar out of range)
map_style_legend() +
fill_color_discrete(cumulative_colors, "Low Income Residents",
createRoundedLabels(inc %>% filter(st_intersects(., neighborhoods %>% filter(name == 'Inner Sunset'), sparse = FALSE)), 'Under200PC', 4, 'jenks'))
# ggsave('output/graphics/inner_sunset.png', plot = is, device = 'png', width = 8, height = 8)
# Fixing issue # 2
# Setting all maps to be square, despite shape of neighborhood
# The following sequence is used to create the bounds of the map. Because each neighborhood is a different shape,
# this needs to be carried out in order to create a buffer around each neighborhood and to make the map extent a square
# regardless of neighborhood shape.
local_crs <- 2227
# find bounding coorindates
bbox <- st_bbox(neighborhoods %>% filter(name == "Inner Sunset"))
# create four corners of actual bounding polygon
points = st_sfc(st_point(c(bbox$xmin, bbox$ymax)), st_point(c(bbox$xmin, bbox$ymin)),
st_point(c(bbox$xmax, bbox$ymin)), st_point(c(bbox$xmax, bbox$ymax)), crs=local_crs) #make sure to appropriately set crs
# test
plot(points)
# find diagonal distance of coordinates of bounding polygon
diagonal <- st_distance(points[1], points[3])
# create the buffer distance, which is half of the diagonal distance of coordinates
buffer_dist <- diagonal/2
# convert points into a data.frame
points_df <- data.frame(st_coordinates(points))
# add the first point coordinates again to the end, this will make the polygon able to close
points_df <- rbind(points_df, points_df[1,])
# create an actual bounding polygon, not just coordinates
poly <- st_sf(
st_sfc(
st_polygon(
list(as.matrix(points_df)))), crs = local_crs)
#test
plot(poly)
# find centroid of polygon
poly_centroid <- st_centroid(poly)
#test
ggplot() +
geom_sf(data = points, color = 'red') +
geom_sf(data = poly, color = "blue", fill = NA) +
geom_sf(data = poly_centroid, color = "green")
# create buffer of centroid using the the half diagonal length as distance
the_shape_buffer = st_buffer(poly_centroid, dist = buffer_dist)
#test
ggplot() +
geom_sf(data = points, color = 'red') +
geom_sf(data = poly, color = "blue", fill = NA) +
geom_sf(data = poly_centroid, color = "green") +
geom_sf(data = the_shape_buffer, color = "purple", fill = NA)
# create coordinates for the map range, using the buffer
mapRange_inner_sunset_v2 <- c(range(st_coordinates(the_shape_buffer)[,1]),
range(st_coordinates(the_shape_buffer)[,2]))
ggplot() +
geom_sf(data = createBins(inc %>% filter(st_intersects(., neighborhoods %>% filter(name == 'Inner Sunset'), sparse = FALSE)), 'Under200PC', 4, 'jenks'), aes(fill = bins), colour = NA) + # colour = NA remove the lines width, change to size = 1 (or whatever number) to get a line width
geom_sf(data = neighborhoods, size = 1, fill = NA) +
coord_sf(xlim = mapRange_inner_sunset_v2[c(1:2)], ylim = mapRange_inner_sunset_v2[c(3:4)]) +
map_style_legend() +
fill_color_discrete(cumulative_colors, "Low Income Residents",
createRoundedLabels(inc %>% filter(st_intersects(., neighborhoods %>% filter(name == 'Inner Sunset'), sparse = FALSE)), 'Under200PC', 4, 'jenks'))
# The above code to set the square map range could (and probably should) be converted
# to a function if you're using it for more than one map
# One thing to nice too is that the block groups now look a little weird because the map range is
# bigger than the block groups that intersect with the neighborhood boundary, you can fix this
# by setting them to equal the intersection with the buffer or by some more manual solution (intersecting neighborhoring neighborhoods for example)
# not perfect, but looks better
# is2 <-
ggplot() +
geom_sf(data = createBins(inc %>% filter(st_intersects(., the_shape_buffer, sparse = FALSE)), 'Under200PC', 4, 'jenks'), aes(fill = bins), colour = NA) + # colour = NA remove the lines width, change to size = 1 (or whatever number) to get a line width
geom_sf(data = neighborhoods, size = 1, fill = NA) +
coord_sf(xlim = mapRange_inner_sunset_v2[c(1:2)], ylim = mapRange_inner_sunset_v2[c(3:4)]) +
map_style_legend() +
fill_color_discrete(cumulative_colors, "Low Income Residents",
createRoundedLabels(inc %>% filter(st_intersects(., neighborhoods %>% filter(name == 'Inner Sunset'), sparse = FALSE)), 'Under200PC', 4, 'jenks'))
# ggsave('output/graphics/inner_sunset.svg', plot = is2, device = 'svg', width = 8, height = 8)
# ggsave('output/graphics/inner_sunset.eps', plot = is2, device = 'eps', width = 8, height = 8)
# mapsf -------------------------------------------------------------------
library(mapsf)
# examples from mapsf -----------------------------------------------------
mtq <- mf_get_mtq()
# Plot the base map
mf_map(x = mtq)
# Plot proportional symbols
mf_map(x = mtq, var = "POP", type = "prop")
# Plot a map layout
mf_layout(title = "Population in Martinique",
credits = "T. Giraud; Sources: INSEE & IGN, 2018")
# Plot choropleth
mf_map(x = mtq, var = "POP", type = "choro")
mf_layout(title = "Population in Martinique",
credits = "T. Giraud; Sources: INSEE & IGN, 2018")
# More detailed example
# Export a map figure with a theme and extra margins
mf_export(x = mtq, export = "svg", filename = "output/graphics/mtq.svg", width = 5,
theme = "dark", expandBB = c(0,0,0,.3))
# Plot a shadow
mf_shadow(mtq, col = "grey10", add = TRUE)
# Plot a choropleth map
mf_map(x = mtq, var = "MED", type = "choro",
pal = "Dark Mint",
breaks = "quantile",
nbreaks = 6,
leg_title = "Median Income\n(euros)",
leg_val_rnd = -2,
add = TRUE)
# Start an inset map
mf_inset_on(x = "worldmap", pos = "right")
# Plot the position of the sample dataset on a worlmap
mf_worldmap(mtq, col = "#0E3F5C")
# Close the inset
mf_inset_off()
# Plot a title
mf_title("Wealth in Martinique, 2015")
# Plot credits
mf_credits("T. Giraud\nSources: INSEE & IGN, 2018")
# Plot a scale bar
mf_scale(size = 5)
# Plot a north arrow
mf_arrow('topleft')
dev.off()
# SF Map ------------------------------------------------------------------
# recreating our inner sunset map
# define a target
target <- neighborhoods %>% filter(name == 'Inner Sunset')
# OR we can define the extent of the map using the shape buffer we created above
# Export a map figure with a theme and extra margins
# mf_export(x = the_shape_buffer, export = "svg", filename = "output/graphics/inner_sunset_mapsf.svg", theme = "default")
mf_init(the_shape_buffer, theme = "default")
# display target
mf_map(inc %>% filter(st_intersects(., neighborhoods %>% filter(name == 'Inner Sunset'), sparse = FALSE)),
var = "Under200PC", type = "choro",
pal = 'ag_GrnYl',
border = NA,
leg_pos = "topleft",
leg_title = "Low Income Residents",
add = TRUE)
# display all neighborhoods
mf_map(neighborhoods, col = NA,
border = "white", add = TRUE)
# display target
mf_map(target, col = NA,
border = "white", lwd = 5, add = TRUE)
# display streets
mf_map(streets, col = "gray", lwd = 1, add = TRUE)
# notice that drawing order matters here too and the legend is linked to the choropleth layer,
# so the streets draw on top of the legend --- not sure what the solution is here yet
# dev.off()
# This is a new package for me and there are a lot of settings I don't know about yet
# As you can see, it has a wider default than ggplot2, so even initializing the map at the extent
# of the "square" buffer of inner sunset results in a rectangular map because it adds a margin
# to the left and right. This can likely be removed by tweaking some settings (not sure which yet)
# Also notice that the choro setting automatically bins the data rather than creating a
# continuous stretch like ggplot2
# there are settings to manually set breaks, see below
mf_init(the_shape_buffer, theme = "default")
# display target
mf_map(x = inc %>% filter(st_intersects(., neighborhoods %>% filter(name == 'Inner Sunset'), sparse = FALSE)),
var = "Under200PC", type = "choro",
pal = 'ag_GrnYl',
breaks = 'jenks',
nbreaks = 4,
border = NA,
leg_pos = "topleft",
leg_title = "Low Income Residents",
add = TRUE)
# display all neighborhoods
mf_map(neighborhoods, col = NA,
border = "white", add = TRUE)
# display target
mf_map(target, col = NA,
border = "white", lwd = 5, add = TRUE)