Data Science II

Session 4: Spatial Joins, Measurement, and Buffers

Bogdan G. Popescu
bogdan.popescu@johncabot.edu

John Cabot University

What You’ll Learn Today

  1. st_join — combining layers by location instead of by a shared code
  2. Measurement: st_area, st_centroid, st_distance, st_nearest_feature
  3. Buffers, unions, and intersections — the overlay toolkit
  4. A complete access-to-services analysis, end to end

The question we answer today: how far is each municipio in Yucatán from its nearest hospital?

Before We Start

#Step1: Libraries for today
library(sf)
library(dplyr)
library(ggplot2)
library(patchwork)

sf_use_s2(FALSE)

Data, as before, in ../data/:

  • denue/denue_salud_yucatan.csv — 6,460 health establishments
  • gadm41_MEX_shp/ — municipio boundaries

The complete, ordered code is in session4_code.R.

Where We Left Off

In Session 3 we joined homicides to boundaries on CVEGEO — a shared code present in both tables.

That is an ordinary left_join. Nothing spatial happened.

But often there is no shared code. You have coordinates on one side and polygons on the other, and the only thing connecting them is where they are.

Part A: Joining by Location

Join by Location

A spatial join matches rows from two layers by a geometric relationship — typically “this point falls inside that polygon”.

The sf Operations Toolkit

Everything today comes from this family. They all begin st_, they all take sf objects, and they all return sf objects.

The Setup

#Step1: The points -- health establishments with coordinates
d   <- read.csv("../data/denue/denue_salud_yucatan.csv")
pts <- st_as_sf(d, coords = c("longitud", "latitud"), crs = 4326)

#Step2: The polygons -- Yucatán municipios
mex <- st_read("../data/gadm41_MEX_shp/gadm41_MEX_2.shp", quiet = TRUE)
yuc <- mex %>% filter(NAME_1 == "Yucatán") %>% select(NAME_2)

nrow(pts); nrow(yuc)
[1] 6460
[1] 106

6,460 points, 106 polygons, and no shared key — the DENUE file has a municipio name, but we are going to ignore it and use geometry instead.

st_join

#Step1: Attach to each point the polygon it falls inside
joined <- st_join(pts, yuc, join = st_within)

#Step2: Every point now carries a NAME_2 from the polygon layer
joined %>% st_drop_geometry() %>% select(nom_estab, NAME_2) %>% head(3)
  nom_estab NAME_2
1           Mérida
2           Mérida
3           Mérida

st_join(x, y) keeps every row of x and adds columns from y. It is a left join, in spatial clothing.

Choosing the Predicate

Predicate Matches when
st_within x is entirely inside y
st_intersects x and y touch or overlap at all (the default)
st_contains x entirely contains y
st_nearest_feature the closest y, whatever the distance

For points in polygons, st_within and st_intersects give the same answer except exactly on a boundary.

Tip

st_nearest_feature never returns NA — it always finds something. That is useful, and occasionally dangerous.

Counting Points per Polygon

#Step1: Drop the geometry and count -- ordinary dplyr from here
spatial_count <- joined %>%
  st_drop_geometry() %>%
  filter(!is.na(NAME_2)) %>%
  count(NAME_2, name = "spatial")

head(spatial_count, 4)
   NAME_2 spatial
1   Abalá       5
2 Acanceh      19
3    Akil      32
4    Baca      22

This is the same number we computed in Session 3 — but derived from coordinates rather than from a text label.

Do the Two Methods Agree?

#Step1: The attribute-based count, as in Session 3
attr_count <- d %>% count(municipio, name = "attribute")

#Step2: Compare them side by side
cmp <- full_join(spatial_count, attr_count,
                 by = c("NAME_2" = "municipio")) %>%
  mutate(diff = spatial - attribute)

sum(cmp$diff != 0, na.rm = TRUE)
[1] 4

Four municipios out of 106 disagree. This is the interesting part of the lesson.

Where They Disagree

cmp %>% filter(diff != 0) %>% arrange(desc(abs(diff)))
      NAME_2 spatial attribute diff
1    Kanasín      78        85   -7
2     Mérida    4218      4212    6
3 San Felipe       4         5   -1
4       Umán      95        94    1

Seven establishments that DENUE labels Kanasín have coordinates inside Mérida. One labelled Mérida falls in Umán. And San Felipe loses one entirely.

The Metropolitan Fringe

#Step1: Which label-versus-geometry mismatches are there?
joined %>% st_drop_geometry() %>%
  filter(!is.na(NAME_2), municipio != NAME_2) %>%
  count(municipio, NAME_2, name = "n")
  municipio NAME_2 n
1   Kanasín Mérida 7
2    Mérida   Umán 1

Eight points in total, all on the Mérida metropolitan boundary — exactly where an address-based label and a GPS coordinate are most likely to part company.

The Orphan

#Step1: Points that fell inside no polygon at all
orphan <- joined %>% filter(is.na(NAME_2))
orphan %>% st_drop_geometry() %>% select(nom_estab, municipio)
                nom_estab  municipio
1 CONSULTORIO MEDICO 24/7 San Felipe
#Step2: How far outside is it?
round(as.numeric(min(st_distance(st_transform(orphan, 6372),
                                 st_transform(yuc, 6372)))))
[1] 34

Thirty-four metres. A coastal clinic in San Felipe, sitting just outside GADM’s simplified shoreline.

Which Answer Do You Trust?

Neither is simply “right”:

  • The label reflects the administrative address — which is what a registry records
  • The geometry reflects where the GPS reading fell — and GADM’s boundary is itself simplified

Important

Report the discrepancy rather than hiding it. Eight points out of 6,460 is 0.1% and changes nothing; if it had been 8%, your whole analysis would rest on which boundary file you happened to download.

Sharing Attributes the Other Way

#Step1: Reverse the arguments -- polygons gain a column from points
#       Each polygon matches MANY points, so rows multiply
hospitals <- d %>% filter(substr(as.character(codigo_act), 1, 4) == "6221") %>%
  st_as_sf(coords = c("longitud", "latitud"), crs = 4326)

yuc_h <- st_join(yuc, hospitals["nom_estab"], join = st_intersects)
nrow(yuc); nrow(yuc_h)
[1] 106
[1] 167

106 polygons became 167 rows — one per polygon-point pair, plus a row of NA for each polygon with no hospital.

Warning

A spatial join can increase your row count. Always check nrow() afterwards. If you want one row per polygon, aggregate.

The CRS Error, Again

#Step1: Two layers, two different coordinate systems
st_join(pts, st_transform(yuc, 6372))
Error in `st_geos_binop()`:
! st_crs(x) == st_crs(y) is not TRUE

The error from Session 1, exactly as promised. st_transform one of them and it goes away.

Logical Predicates

st_intersects and friends return a sparse list, not a matrix — one element per feature in x, holding the indices of matching features in y.

#Step1: For each municipio, which hospitals fall inside it?
hits <- st_intersects(yuc, hospitals)
hits[1:3]
[[1]]
[1] 97

[[2]]
[1] 35

[[3]]
integer(0)

Municipio 1 has no hospital; municipio 2 has one, and it is hospital number 3 in the hospitals layer.

The lengths() Idiom

#Step1: lengths() turns the sparse list into a count per polygon
yuc$n_hospitals <- lengths(st_intersects(yuc, hospitals))

#Step2: How many municipios have at least one hospital?
sum(yuc$n_hospitals > 0)
[1] 45
table(yuc$n_hospitals)

 0  1  2  3  5  7 35 
61 28 11  3  1  1  1

Only 45 of 106 municipios contain a hospital at all. Sixty-one have none.

Tip

lengths(st_intersects(a, b)) is the fastest way to count b-features per a-feature. Learn it — it replaces a whole st_join plus group_by plus summarise pipeline.

Spatial Filtering

#Step1: Keep only the municipios that contain a hospital
with_hosp <- yuc %>% filter(lengths(st_intersects(., hospitals)) > 0)
nrow(with_hosp)
[1] 45
#Step2: st_filter does the same thing more directly
nrow(st_filter(yuc, hospitals))
[1] 45

st_filter(x, y) keeps the rows of x that satisfy a predicate against y — like st_join, but it adds no columns and never duplicates rows.

Aggregating with a Spatial Join

#Step1: Join, then summarise -- the general pattern when you need
#       something other than a count
by_mun <- st_join(pts, yuc["NAME_2"], join = st_within) %>%
  st_drop_geometry() %>%
  filter(!is.na(NAME_2)) %>%
  group_by(NAME_2) %>%
  summarise(n = n(),
            share_with_phone = mean(telefono != "" & !is.na(telefono)),
            .groups = "drop")

head(by_mun, 4)
# A tibble: 4 × 3
  NAME_2      n share_with_phone
  <chr>   <int>            <dbl>
1 Abalá       5            0    
2 Acanceh    19            0.632
3 Akil       32            0.656
4 Baca       22            0.591

Once the join is done, everything is ordinary dplyr. The spatial part is one line.

Part B: Measurement

Three Kinds of Calculation

  • Numeric — returns a number: st_area, st_distance, st_length
  • Logical — returns TRUE/FALSE: st_intersects, st_within, st_disjoint
  • Spatial — returns new geometry: st_centroid, st_buffer, st_union

All three need a projected CRS if you want answers in metres.

st_area

#Step1: Project first -- EPSG:6372 has units of metres
yuc_p <- st_transform(yuc, 6372)

#Step2: Area comes back with units attached
head(st_area(yuc_p), 3)
Units: [m^2]
[1] 290237841 136063311  75871595
#Step3: Convert to km2 and strip the units for plotting
yuc_p$area_km2 <- as.numeric(units::set_units(st_area(yuc_p), km^2))
round(sum(yuc_p$area_km2))
[1] 39102

39,102 km² — the published area of Yucatán is about 39,600. The 1% gap is boundary simplification.

st_centroid

cen <- st_centroid(yuc_p)
nrow(cen)
[1] 106

One point per polygon — the geometric centre of mass.

Centroids Can Lie

An L-shaped municipio.

Centroids Can Lie

st_centroid() — the centre of mass falls outside the polygon.

Centroids Can Lie

st_point_on_surface() — guaranteed to fall inside.

When That Matters

The centroid of a crescent-shaped or very concave polygon can fall outside the polygon itself.

#Step1: How many centroids fall outside their own municipio?
sum(!as.logical(st_within(cen, yuc_p, sparse = FALSE) %>% diag()))
[1] 2

None here — Yucatán’s municipios are compact. But for a coastal state with islands, use st_point_on_surface(), which guarantees a point inside.

st_distance

#Step1: The hospitals, projected
hosp <- hospitals %>% st_transform(6372)
nrow(hosp)
[1] 106
#Step2: By default st_distance returns a full MATRIX
dim(st_distance(cen[1:3, ], hosp[1:4, ]))
[1] 3 4

106 centroids × 106 hospitals would be 11,236 distances. Usually you do not want them all.

st_nearest_feature

#Step1: For each centroid, the index of the closest hospital
nearest <- st_nearest_feature(cen, hosp)
head(nearest)
[1]  97  35  62  53  72 102
#Step2: by_element gives a vector, not a matrix -- one distance per pair
cen$km_to_hospital <- as.numeric(
  st_distance(cen, hosp[nearest, ], by_element = TRUE)) / 1000

round(range(cen$km_to_hospital), 1)
[1]  0.5 44.2

Between 0.5 km and 44.2 km to the nearest hospital.

Who Is Furthest?

cen %>% st_drop_geometry() %>%
  arrange(desc(km_to_hospital)) %>%
  transmute(NAME_2, km = round(km_to_hospital, 1)) %>%
  head(6)
        NAME_2   km
1 Chikindzonot 44.2
2 Río Lagartos 43.5
3   San Felipe 41.2
4      Tizimín 32.7
5       Chemax 32.5
6     Buctzotz 31.7

Chikindzonot, Río Lagartos, San Felipe — the eastern interior and the north coast.

Mapping Accessibility

#Step1: Put the distance back on the polygons
yuc_p$km_to_hospital <- cen$km_to_hospital

#Step2: Map it, hospitals on top
ggplot() +
  geom_sf(data = yuc_p,
          aes(fill = km_to_hospital),
          colour = "white",
          linewidth = 0.1) +
  geom_sf(data = hosp, size = 0.5,
          colour = "white") +
  scale_fill_viridis_c(direction = -1,
                       name = "km") +
  labs(title = "Distance to the nearest hospital") +
  theme_void(base_size = 9)

Part C: Buffers, Unions, and Intersections

Buffers

st_buffer(x, dist) grows a geometry outward by dist — turning a point into a circle, a line into a corridor, a polygon into a larger polygon.

Buffers Need Metres

#Step1: The layer is already in EPSG:6372, so dist is in METRES
buf10 <- st_buffer(hosp, 10000)   # 10 km
nrow(buf10)
[1] 106

Warning

The number you pass to dist is in the units of the layer’s CRS. In EPSG:6372 that is metres, so 10000 is 10 km. In EPSG:4326 it would be 10,000 degrees — twenty-seven times around the planet.

Modern sf will quietly do the geodesically correct thing on a geographic CRS, which is worse in a way: you get a plausible answer without ever learning you were sloppy.

st_union: Dissolving

#Step1: 106 overlapping circles become one multipolygon
covered <- st_union(buf10)
length(covered)
[1] 1

Without the union, overlapping buffers would double-count the area they share.

st_intersection

#Step1: Clip the coverage to the state -- circles spill into the sea
state    <- st_union(yuc_p)
coverage <- st_intersection(state, covered)

st_difference

#Step1: The complement -- everywhere NOT within 10 km of a hospital
gap <- st_difference(state, covered)

intersection is “in both”, difference is “in the first but not the second”, union is “in either”. Those three cover most overlay work.

How Much of the State Is Covered?

#Step1: Areas in km2
a_state <- as.numeric(st_area(state))   / 1e6
a_cov   <- as.numeric(st_area(coverage)) / 1e6

#Step2: The share
round(c(state_km2 = a_state, covered_km2 = a_cov,
        percent = 100 * a_cov / a_state), 1)
  state_km2 covered_km2     percent 
    39102.0     11685.9        29.9

Only 29.9% of Yucatán is within 10 km of a hospital.

The Coverage Map

ggplot() +
  geom_sf(data = state,
          fill = "#f2d7d5",
          colour = NA) +
  geom_sf(data = coverage,
          fill = "#4c9f70",
          colour = NA,
          alpha = 0.85) +
  geom_sf(data = hosp, size = 0.4,
          colour = "white") +
  labs(title = "Within 10 km of a hospital",
       subtitle = "Green: covered. Pink: not.") +
  theme_void(base_size = 9)

Other Shapes You Can Derive

st_convex_hull

The tightest convex shape containing the points — a crude “service area”.

st_sym_difference

In one or the other, but not both.

Geometry Casting

#Step1: The state is one MULTIPOLYGON
st_geometry_type(state)
[1] MULTIPOLYGON
18 Levels: GEOMETRY POINT LINESTRING POLYGON MULTIPOINT ... TRIANGLE
#Step2: st_cast breaks it into its component POLYGONs
pieces <- st_cast(state, "POLYGON")
length(pieces)
[1] 6

Casting is how you split a multi-part feature apart — to find the mainland and discard the islands, say — or how you go from polygons to their boundary lines with st_cast(x, "LINESTRING").

How Sensitive Is the Answer?

#Step1: Recompute coverage at three thresholds
coverage_at <- function(d) {
  cov <- st_intersection(state, st_union(st_buffer(hosp, d)))
  100 * as.numeric(st_area(cov)) / as.numeric(st_area(state))
}

round(sapply(c(`5km` = 5000, `10km` = 10000, `20km` = 20000),
             coverage_at), 1)
 5km 10km 20km 
10.2 29.9 60.1

10.2% at 5 km, 29.9% at 10 km, 60.1% at 20 km.

Important

The headline number is entirely a function of the threshold you chose. State the threshold, justify it, and show the sensitivity — otherwise you are picking the number that suits your argument.

Part D: Putting It Together

The Whole Analysis

Five operations, in order:

  1. st_as_sf — coordinates become geometry (Session 2)
  2. st_transform — into metres, EPSG:6372
  3. st_centroid — one representative point per municipio
  4. st_nearest_feature + st_distance — how far to the nearest hospital
  5. st_buffer + st_union + st_intersection — what share of the territory is served

That is a complete accessibility study, in about twenty lines.

What We Found

  • 106 hospitals for 106 municipios — but only 45 municipios contain one
  • Distance from a municipio centroid to its nearest hospital runs from 0.5 to 44.2 km
  • 29.9% of the state’s territory lies within 10 km of a hospital
  • The worst-served are Chikindzonot, Río Lagartos, and San Felipe

What We Have Not Found

Warning

Territory is not people. 70% of the area is more than 10 km from a hospital, but far less than 70% of the population — people cluster near services. Do not report the area figure as if it were a population figure.

Other limits worth stating in your reto:

  • A centroid is not where people live
  • Straight-line distance is not travel time
  • DENUE records establishments, not capacity — a two-bed clinic counts the same as a regional hospital

What Have We Learned?

  • st_join combines layers by location; the predicate (st_within, st_intersects, st_nearest_feature) decides what “matches” means
  • A spatial join can change your row count — always check
  • Joining by code and joining by location disagreed for 4 of 106 municipios, all on the metropolitan fringe, plus one point 34 m offshore
  • Measurement needs a projected CRS: st_area, st_distance, and st_buffer all take their units from it
  • st_nearest_feature + by_element = TRUE is how you get one distance per row instead of a matrix
  • union → intersection → area is the standard recipe for “what share is covered”
  • Area covered is not population covered

Next session: rasters — continuous surfaces, and extracting them to your polygons.