How to add multiple points to map with tmap in R?

Question

I'm learning how to make maps with R so I'm kinda new to all of this. I'm using tmap to plot a map (you can find the shp here https://www.data.rio/datasets/%C3%A1reas-program%C3%A1ticas-da-sa%C3%BAde/explore).

I also have a dataframe that looks like this:

data <- data.frame(COD_AP_SMS = c("AP 1.0", "AP 2.1", "AP 2.2", "AP 3.1", "AP 3.2", "AP 3.3", "AP 4.0", "AP 5.1", "AP 5.2", "AP 5.3"),
                   n = c(11, 21, 6, 32, 15, 32, 27, 26, 19, 18))

I'm trying to plot a map where the colors of the polygons are different for each COD_AP_SMS and add n points to each polygon. For example, AP 1.0 would plot 11 points in it's polygon. Is it possible? Here's my code

ap <- read_sf("shp/Limites AP.shp")

plot(ap)

  
ap <- full_join(ap, tb, by = "COD_AP_SMS")

  
# make some bbox magic
bbox_new <- st_bbox(ap) # current bounding box

xrange <- bbox_new$xmax - bbox_new$xmin # range of x values
yrange <- bbox_new$ymax - bbox_new$ymin # range of y values

# bbox_new[1] <- bbox_new[1] - (0.25 * xrange) # xmin - left
bbox_new[3] <- bbox_new[3] + (0.25 * xrange) # xmax - right
# bbox_new[2] <- bbox_new[2] - (0.25 * yrange) # ymin - bottom
bbox_new[4] <- bbox_new[4] + (0.2 * yrange) # ymax - top

bbox_new <- bbox_new %>%  # take the bounding box ...
  st_as_sfc() # ... and make it a sf polygon

   
tm_shape(ap, bbox = bbox_new)+
  tm_fill("COD_AP_SMS", auto.palette.mapping=FALSE, 
          title="AP")+
  tm_dots("n", size = 2) +
  tm_compass(position = c("left", "top"))+
  tm_scale_bar(position = c("left", "top"))+
  tm_borders(alpha=1, lwd = 1, col = "black")+
  tm_legend(legend.format = list(text.separator= "a")) +
  tm_layout(legend.position = c("right", "top"),
            frame = FALSE)

Tips on how to make laebls, scale bar and compass look better are welcome.

The desired output would look something like this:

Answer 1

You can add shapes from multiple objects to a tmap by adding a new tm_shape() for each. To generate a set of random points, sf provides st_sample() function which is also able to generate points by polygon.

library(dplyr)
library(tmap)
library(sf)
#> Linking to GEOS 3.11.2, GDAL 3.6.2, PROJ 9.2.0; sf_use_s2() is TRUE


data <- data.frame(COD_AP_SMS = c("AP 1.0", "AP 2.1", "AP 2.2", "AP 3.1", "AP 3.2", "AP 3.3", "AP 4.0", "AP 5.1", "AP 5.2", "AP 5.3"),
                   n = c(11, 21, 6, 32, 15, 32, 27, 26, 19, 18))

# shapes from https://www.data.rio/datasets/%C3%A1reas-program%C3%A1ticas-da-sa%C3%BAde/explore
# there's an invalid shape, in both GeoJSON and SHP formats, pass it through st_make_valid
ap <- 
  read_sf("/vsizip/81reas_ProgramA1ticas_da_SaBAde.zip") |> 
  st_make_valid() |>
  left_join(data, by = "COD_AP_SMS")
ap
#> Simple feature collection with 10 features and 6 fields
#> Geometry type: GEOMETRY
#> Dimension:     XY
#> Bounding box:  xmin: 623577.4 ymin: 7446298 xmax: 695396.8 ymax: 7483424
#> Projected CRS: SAD69 / UTM zone 23S
#> # A tibble: 10 × 7
#>    OBJECTID COD_AP_SMS GlobalID  Shape__Are Shape__Len                  geometry
#>       <int> <chr>      <chr>          <dbl>      <dbl>            <GEOMETRY [m]>
#>  1       11 AP 1.0     ca246675…  34395282.     91084. MULTIPOLYGON (((682126.9…
#>  2       12 AP 2.1     b1021547…  45267751.     91079. MULTIPOLYGON (((683530.6…
#>  3       13 AP 2.2     7c5a2fb0…  55165977.     49175. POLYGON ((683708.1 74653…
#>  4       14 AP 3.1     7ca38ee3…  85356489.    127247. MULTIPOLYGON (((682960.1…
#>  5       15 AP 3.2     b35f71c1…  41235758.     35087. POLYGON ((676067 7471346…
#>  6       16 AP 3.3     156993af…  76899060.     58051. POLYGON ((669091.1 74774…
#>  7       17 AP 4.0     bc09989a… 293782989.    112591. MULTIPOLYGON (((652583.1…
#>  8       18 AP 5.1     f0551e61… 116624827.     59923. POLYGON ((657109.4 74733…
#>  9       19 AP 5.2     7ad486bd… 291335884.    113941. MULTIPOLYGON (((646656.1…
#> 10       20 AP 5.3     f3e0d549… 164083566.     69515. MULTIPOLYGON (((633723.6…
#> # ℹ 1 more variable: n <dbl>

# generate sample points 
pnts <- st_sample(ap, ap$n, by_polygon = TRUE) |> st_as_sf() 
pnts
#> Simple feature collection with 207 features and 0 fields
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: 626147.9 ymin: 7448069 xmax: 689342.7 ymax: 7478607
#> Projected CRS: SAD69 / UTM zone 23S
#> First 10 features:
#>                           x
#> 1  POINT (685938.4 7466795)
#> 2  POINT (681973.6 7466671)
#> 3  POINT (686089.6 7463638)
#> 4  POINT (681518.6 7467800)
#> 5  POINT (685674.3 7464450)
#> 6  POINT (683766.8 7465139)
#> 7    POINT (684001 7462334)
#> 8  POINT (683373.8 7468592)
#> 9  POINT (683789.5 7470009)
#> 10 POINT (683143.1 7462329)

# make some bbox magic
bbox_new <- st_bbox(ap) # current bounding box
xrange <- bbox_new$xmax - bbox_new$xmin # range of x values
yrange <- bbox_new$ymax - bbox_new$ymin # range of y values
bbox_new[3] <- bbox_new[3] + (0.25 * xrange) # xmax - right
bbox_new[4] <- bbox_new[4] + (0.2 * yrange) # ymax - top

bbox_new <- bbox_new %>%  # take the bounding box ...
  st_as_sfc() # ... and make it a sf polygon

# tmap with 2 layers
tm_shape(ap, bbox = bbox_new)+
  tm_fill("COD_AP_SMS", auto.palette.mapping=FALSE, 
          title="AP")+
  tm_compass(position = c("left", "top"))+
  tm_scale_bar(position = c("left", "top"))+
  tm_borders(alpha=1, lwd = 1, col = "black")+
  tm_legend(legend.format = list(text.separator= "a")) +
  tm_layout(legend.position = c("right", "top"),
            frame = FALSE) +
tm_shape(pnts) +
  tm_dots(size = .1, col = "red")
#> Warning: The argument auto.palette.mapping is deprecated. Please use midpoint
#> for numeric data and stretch.palette for categorical data to control the
#> palette mapping.