Chapter 2 Measuring and constructing shapes

Let’s practice using spatial functions to measure and characterize the shapes of polygons using the example of watersheds. All of these measurements are dependent upon the accuracy of the map projection in which the data is stored.

Also, different from normal GIS, measurements produce a value with a unit in R! This can be helpful, but when it isn’t, it’s easy enough to convert back to a simple number with as.numeric().

2.1 Area and Perimeter

What is the area of each watershed? You’ll need to know which column stores the geographic geometry data.

watersheds_shp <- mutate(watersheds, 
                     warea = st_area(geom),
                     wperim = st_perimeter(geom))
## Warning: There was 1 warning in `stopifnot()`.
## ℹ In argument: `wperim = st_perimeter(geom)`.
## Caused by warning in `st_perimeter()`:
## ! 'st_perimeter' is deprecated.
## Use 'sf::st_perimeter or lwgeom::st_perimeter_lwgeom' instead.
## See help("Deprecated")
watersheds_shp |> st_drop_geometry() |> knitr::kable()
name warea wperim
Moshassuck 1031250 [m^2] 4560.660 [m]
Narragansett 593750 [m^2] 3353.553 [m]
Blackstone 2375000 [m^2] 7914.214 [m]

We can also calculate the compactness as a function of perimeter and area. The Polsby-Popper isoperimetric ratio is 4 * pi * A / P^2

watersheds_shp <- mutate(watersheds_shp, 
                     compact = round(
                       as.numeric(
                         (4 * pi * warea) / wperim^2),
                       2))

watersheds_shp |> st_drop_geometry() |> knitr::kable()
name warea wperim compact
Moshassuck 1031250 [m^2] 4560.660 [m] 0.62
Narragansett 593750 [m^2] 3353.553 [m] 0.66
Blackstone 2375000 [m^2] 7914.214 [m] 0.48

Do the compactness scores make sense? Let’s map them.

pointgrid_map +
  watersheds_shp |> 
    tm_shape() +
    tm_polygons(fill = "compact",
                fill_alpha = 0.5,
                fill.legend = tm_legend_hide()) +
    tm_text("compact")

Looks good to me!

2.2 Centroids

We also have a few methods for converting polygons to points.

Centroids are located at the center of gravity of a polygon.

watersheds_centroids <- watersheds_shp |> st_centroid()

Map the centroid points as blue dots.

watersheds_centroids_map <- watersheds_centroids |> 
  tm_shape() +
  tm_symbols(
    size = 0.8,
    fill = "blue",
    fill_alpha = 0.8
  )

pointgrid_map +
  watersheds_map +
  watersheds_centroids_map

2.3 Point on surface

Centroids sometimes fall outside of the polygon, which might produce undesirable results for cartography or topological analysis. In those cases point_on_surface may be better.

watersheds_pt_on_surf <- watersheds_shp |> st_point_on_surface()

Let’s see the points on surface in relation to the centroids.

watersheds_points_on_surface_map <- watersheds_pt_on_surf |> 
  tm_shape() +
  tm_symbols(
    size = 0.8,
    fill = "white",
    fill_alpha = 0.8
  )

pointgrid_map +
  watersheds_map +
  watersheds_centroids_map +
  watersheds_points_on_surface_map

2.4 Convex hull

Some measures of gerrymandering rely on the ratio of the area to the area of the convex hull or of the minimum bounding circle (Niemi et al. 1990).

Let’s calculate the convex hull.

watersheds_hull <- watersheds_shp |> st_convex_hull()

And visualize the resulting polygons.

pointgrid_map +
  watersheds_map +
  tm_shape(watersheds_hull) +
  tm_polygons(
    col = "white",
    lwd = 3,
    lty = "dotted",
    fill = "red",
    fill_alpha = 0.2
  )

How would the ratios of area to area of the convex hull come out? The watersheds contributing to more overlapping areas whould result in lower compactness scores.

watersheds_hull <- watersheds_hull |> 
  mutate(hullarea = st_area(geom),
         compact_hull = round(as.numeric(warea / hullarea), 2))

watersheds_hull |> st_drop_geometry() |> knitr::kable()
name warea wperim compact hullarea compact_hull
Moshassuck 1031250 [m^2] 4560.660 [m] 0.62 1031250 [m^2] 1.00
Narragansett 593750 [m^2] 3353.553 [m] 0.66 656250 [m^2] 0.90
Blackstone 2375000 [m^2] 7914.214 [m] 0.48 2906250 [m^2] 0.82

2.5 Minimum bounding circle

Another method of measuring compactness is to find the ratio of the area of the polygon and divide it by the area of its minimum bounding circle (Reock 1961).

Let’s find the minimum bounding circles of watersheds. This function is not yet included in the current sf package on CRAN, but it is in lwgeom.

watersheds_boundcircle <- watersheds_shp |> st_minimum_bounding_circle()

Let’s map the resulting circle polygons.

pointgrid_map +
  watersheds_map +
  tm_shape(watersheds_boundcircle) +
  tm_polygons(
    col = "white",
    lwd = 3,
    lty = "dotted",
    fill = "red",
    fill_alpha = 0.2
  )

Finally, let’s calculate the ratio of area of the polygon to the area of its minimum bounding circle.

watersheds_boundcircle <- watersheds_boundcircle |> 
  mutate(watersheds_boundcircle,
         mbcarea = st_area(geom),
         compact_circ = round(as.numeric(warea / mbcarea), 2))

watersheds_boundcircle |> st_drop_geometry() |> knitr::kable()
name warea wperim compact mbcarea compact_circ
Moshassuck 1031250 [m^2] 4560.660 [m] 0.62 2845768 [m^2] 0.36
Narragansett 593750 [m^2] 3353.553 [m] 0.66 1226624 [m^2] 0.48
Blackstone 2375000 [m^2] 7914.214 [m] 0.48 6280315 [m^2] 0.38

2.6 Summarize compactness

And summarize the various metrics into one table by joining. First, drop geometry data for the table to be joined, and select the columns to be used by the join. Then, join the two tables with a left_join. Left joins keep all of the records in the left table and only the matching records from the right table to be joined.

watersheds_boundcircle_compactness <- watersheds_boundcircle |> 
  st_drop_geometry() |> 
  select(name, compact_circ)

watersheds_compactness_summary <- watersheds_hull |> 

  select(name, compact_shp = compact, compact_hull) |> 
  left_join(watersheds_boundcircle_compactness, by = "name")

watersheds_compactness_summary |> st_drop_geometry() |> kable()
name compact_shp compact_hull compact_circ
Moshassuck 0.62 1.00 0.36
Narragansett 0.66 0.90 0.48
Blackstone 0.48 0.82 0.38