Austin Full + Limited Purpose Census Tracts

CRG CHeck – Notes below = KW

1 Sources

Downloaded the City of Austin jurisdictions shapefile here: https://catalog.data.gov/dataset/boundaries-city-of-austin-jurisdiction-and-regulatory-boundaries

Downloaded Texas census tracts shapefile from TIGER/line here: https://www.census.gov/geographies/mapping-files/time-series/geo/tiger-line-file.html/

2 Notes

I am filtering for only full and limited purpose jurisdictions.

I decided to do this analysis in R because performing an intersecting overlay or spatial join in ArcGIS would have cost 84 credits.

3 Data Loading and Prep

library(sf)
library(readr)
library(dplyr)

setwd("C:/Users/rayo-garza/OneDrive - Center for Public Policy Priorities/Documents/Data 2025")
getwd()

[1] "C:/Users/rayo-garza/OneDrive - Center for Public Policy Priorities/Documents/Data 2025"

library(sf)

austin_jur <- st_read(
  "C:/Users/rayo-garza/OneDrive - Center for Public Policy Priorities/Documents/Data 2025/geo_export_1413bf5f-ae81-4181-9e51-2348a7e3b5ea.shp",
  quiet = TRUE
)

names(austin_jur)

 [1] "objectid"   "jurisdicti" "city_name"  "jurisdic_2" "modified_f"
 [6] "jurisdic_3" "jurisdic_4" "shape_area" "shape_leng" "globalid"  
[11] "geometry"  

unique(austin_jur$jurisdic_2)

[1] "AUSTIN 2 MILE ETJ"                   "AUSTIN LTD"                         
[3] "AUSTIN 5 MILE ETJ"                   "AUSTIN ETJ AG DEVELOPMENT AGREEMENT"
[5] "AUSTIN FULL PURPOSE"                

#filter to Austin LTD + Full Purpose
austin_filtered <- austin_jur %>%
  filter(jurisdic_2 %in% c("AUSTIN LTD", "AUSTIN FULL PURPOSE")) %>%
  st_make_valid()

#union into ONE geometry (prevents duplicates)
austin_union <- st_union(austin_filtered)

#read Texas tracts
tx_tracts <- st_read(
  "C:/Users/rayo-garza/OneDrive - Center for Public Policy Priorities/Documents/Data 2025/tl_2025_48_tract.shp",
  quiet = TRUE
) %>%
  st_transform(st_crs(austin_union)) %>%
  st_make_valid()


#select tracts whose centroid is W/IN Austin boundaries
selected_tracts <- tx_tracts[
  st_within(st_centroid(tx_tracts), austin_union, sparse = FALSE),
]


st_write(selected_tracts,
         "C:/Users/rayo-garza/OneDrive - Center for Public Policy Priorities/Documents/Data 2025/Austin_LTD_FULL_Tracts.shp",
         delete_dsn = TRUE)

Deleting source `C:/Users/rayo-garza/OneDrive - Center for Public Policy Priorities/Documents/Data 2025/Austin_LTD_FULL_Tracts.shp' using driver `ESRI Shapefile'
Writing layer `Austin_LTD_FULL_Tracts' to data source 
  `C:/Users/rayo-garza/OneDrive - Center for Public Policy Priorities/Documents/Data 2025/Austin_LTD_FULL_Tracts.shp' using driver `ESRI Shapefile'
Writing 218 features with 13 fields and geometry type Polygon.

nrow(selected_tracts)

[1] 218

selected_tracts %>%
  count(COUNTYFP, name = "n_tracts")

Simple feature collection with 3 features and 2 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -98.0104 ymin: 30.08853 xmax: -97.45043 ymax: 30.52177
Geodetic CRS:  WGS 84
  COUNTYFP n_tracts                       geometry
1      209        1 POLYGON ((-98.01009 30.1316...
2      453      204 POLYGON ((-97.62795 30.2912...
3      491       13 POLYGON ((-97.80904 30.4914...

#total tracts
nrow(selected_tracts)

[1] 218

#county breakdown
county_counts <- selected_tracts %>%
  count(COUNTYFP, name = "n_tracts")

county_counts

Simple feature collection with 3 features and 2 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -98.0104 ymin: 30.08853 xmax: -97.45043 ymax: 30.52177
Geodetic CRS:  WGS 84
  COUNTYFP n_tracts                       geometry
1      209        1 POLYGON ((-98.01009 30.1316...
2      453      204 POLYGON ((-97.62795 30.2912...
3      491       13 POLYGON ((-97.80904 30.4914...

selected_df <- selected_tracts %>% st_drop_geometry()

write_csv(
  selected_df,
  "C:/Users/rayo-garza/OneDrive - Center for Public Policy Priorities/Documents/Data 2025/Austin_LTD_FULL_Tracts.csv"
)

library(sf)
library(dplyr)


tracts_intersect <- tx_tracts[st_intersects(tx_tracts, austin_union, sparse = FALSE), ]

tracts_centroid_within <- tx_tracts[st_within(st_centroid(tx_tracts), austin_union, sparse = FALSE), ]

nrow(tracts_intersect)

[1] 271

nrow(tracts_centroid_within)

[1] 218

tracts_intersect %>% count(COUNTYFP, name="n_tracts")

Simple feature collection with 4 features and 2 fields
Geometry type: GEOMETRY
Dimension:     XY
Bounding box:  xmin: -98.08075 ymin: 30.00357 xmax: -97.37302 ymax: 30.53111
Geodetic CRS:  WGS 84
  COUNTYFP n_tracts                       geometry
1      021        1 POLYGON ((-97.568 30.14146,...
2      209        2 MULTIPOLYGON (((-98.01009 3...
3      453      249 POLYGON ((-97.76561 30.4324...
4      491       19 POLYGON ((-97.72115 30.4558...

tracts_centroid_within %>% count(COUNTYFP, name="n_tracts")

Simple feature collection with 3 features and 2 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -98.0104 ymin: 30.08853 xmax: -97.45043 ymax: 30.52177
Geodetic CRS:  WGS 84
  COUNTYFP n_tracts                       geometry
1      209        1 POLYGON ((-98.01009 30.1316...
2      453      204 POLYGON ((-97.62795 30.2912...
3      491       13 POLYGON ((-97.80904 30.4914...

ids_intersect <- tracts_intersect$GEOID
ids_centroid  <- tracts_centroid_within$GEOID

only_in_intersect <- setdiff(ids_intersect, ids_centroid)

length(only_in_intersect)

[1] 53

tx_tracts %>%
  filter(GEOID %in% only_in_intersect) %>%
  count(COUNTYFP, name="n_border_tracts")

Simple feature collection with 4 features and 2 fields
Geometry type: GEOMETRY
Dimension:     XY
Bounding box:  xmin: -98.08075 ymin: 30.00357 xmax: -97.37302 ymax: 30.53111
Geodetic CRS:  WGS 84
  COUNTYFP n_border_tracts                       geometry
1      021               1 POLYGON ((-97.568 30.14146,...
2      209               1 POLYGON ((-97.82383 30.0792...
3      453              45 MULTIPOLYGON (((-97.63627 3...
4      491               6 MULTIPOLYGON (((-97.78192 3...

library(sf)
library(dplyr)
library(ggplot2)

ids_intersect <- tracts_intersect$GEOID
ids_centroid  <- tracts_centroid_within$GEOID
only_in_intersect <- setdiff(ids_intersect, ids_centroid)

border_tracts <- tx_tracts %>%
  filter(GEOID %in% only_in_intersect)


ggplot() +
  geom_sf(data = austin_union, fill = NA, linewidth = 1) +
  geom_sf(data = tracts_centroid_within, fill = NA, linewidth = 0.2) +
  geom_sf(data = border_tracts, fill = NA, linewidth = 0.7, linetype = "dashed") +
  ggtitle("Austin LTD + Full Purpose: Centroid-within tracts vs Intersect-only border tracts") +
  theme_minimal()

THe above code rendered 218 Census tracts – this is inaccurate. Likely due to the strict bounds I used (has to fall within).

4 TRY AGAIN CRG

library(sf)
library(dplyr)
library(readr)

# --- Read Austin jurisdictions ---
austin_jur <- st_read(
  "C:/Users/rayo-garza/OneDrive - Center for Public Policy Priorities/Documents/Data 2025/geo_export_1413bf5f-ae81-4181-9e51-2348a7e3b5ea.shp",
  quiet = TRUE
)

# Filter to Austin LTD + Full Purpose + make valid
austin_filtered <- austin_jur %>%
  filter(jurisdic_2 %in% c("AUSTIN LTD", "AUSTIN FULL PURPOSE")) %>%
  st_make_valid()

# Union into a single boundary (prevents duplicates)
austin_union <- st_union(austin_filtered)

# --- Read Texas tracts ---
tx_tracts <- st_read(
  "C:/Users/rayo-garza/OneDrive - Center for Public Policy Priorities/Documents/Data 2025/tl_2025_48_tract.shp",
  quiet = TRUE
) %>%
  st_make_valid() %>%
  st_transform(st_crs(austin_union))

# --- Use an equal-area CRS for overlap area test ---
# EPSG:3083 (Texas Centric Albers) is good for area calculations in Texas
austin_a <- st_transform(austin_union, 3083)
tracts_a <- st_transform(tx_tracts, 3083)

# --- Identify tracts with positive overlap area (excludes touch-only) ---
# st_intersection returns only true overlaps (area > 0)
int <- st_intersection(
  tracts_a %>% select(GEOID, COUNTYFP),
  st_sf(geometry = austin_a)
)

# Summarize to unique tract GEOIDs
overlap_ids <- int %>%
  st_drop_geometry() %>%
  distinct(GEOID)

# Final selected tracts (keep full tract geometry)
selected_tracts <- tx_tracts %>%
  filter(GEOID %in% overlap_ids$GEOID)

# --- Counts ---
cat("Total selected tracts:", nrow(selected_tracts), "\n")

Total selected tracts: 271 

county_counts <- selected_tracts %>% count(COUNTYFP, name = "n_tracts")
print(county_counts)

Simple feature collection with 4 features and 2 fields
Geometry type: GEOMETRY
Dimension:     XY
Bounding box:  xmin: -98.08075 ymin: 30.00357 xmax: -97.37302 ymax: 30.53111
Geodetic CRS:  WGS 84
  COUNTYFP n_tracts                       geometry
1      021        1 POLYGON ((-97.568 30.14146,...
2      209        2 MULTIPOLYGON (((-98.01009 3...
3      453      249 POLYGON ((-97.76561 30.4324...
4      491       19 POLYGON ((-97.72115 30.4558...

# --- Export CSV (attributes only) ---
selected_df <- selected_tracts %>% st_drop_geometry()

write_csv(
  selected_df,
  "C:/Users/rayo-garza/OneDrive - Center for Public Policy Priorities/Documents/Data 2025/Austin_LTD_FULL_Tracts.csv"
)

# --- Export shapefile (optional) ---
st_write(
  selected_tracts,
  "C:/Users/rayo-garza/OneDrive - Center for Public Policy Priorities/Documents/Data 2025/Austin_LTD_FULL_Tracts.shp",
  delete_dsn = TRUE
)

Deleting source `C:/Users/rayo-garza/OneDrive - Center for Public Policy Priorities/Documents/Data 2025/Austin_LTD_FULL_Tracts.shp' using driver `ESRI Shapefile'
Writing layer `Austin_LTD_FULL_Tracts' to data source 
  `C:/Users/rayo-garza/OneDrive - Center for Public Policy Priorities/Documents/Data 2025/Austin_LTD_FULL_Tracts.shp' using driver `ESRI Shapefile'
Writing 271 features with 13 fields and geometry type Polygon.

Interactive Map to visualize

library(sf)
library(dplyr)
library(tmap)
library(htmlwidgets)

tmap_mode("view")

# helper: forces leaflet to resize after load (fixes blank widgets on RPubs)
fix_widget <- function(w) {
  htmlwidgets::onRender(
    w,
    "function(el, x) {
       setTimeout(function(){ window.dispatchEvent(new Event('resize')); }, 250);
       setTimeout(function(){ window.dispatchEvent(new Event('resize')); }, 1000);
     }"
  )
}

tiny_threshold <- 0.01  # 1% overlap cutoff

#project to equal-area CRS for overlap calcs
austin_a   <- st_transform(austin_union, 3083)
tracts_a   <- st_transform(tx_tracts, 3083)
selected_a <- st_transform(selected_tracts, 3083)

#calculate overlap share
int <- st_intersection(
  selected_a %>% select(GEOID, COUNTYFP),
  st_sf(geometry = austin_a)
)

overlap_share <- int %>%
  mutate(int_area = as.numeric(st_area(.))) %>%
  st_drop_geometry() %>%
  group_by(GEOID) %>%
  summarise(int_area = sum(int_area), .groups = "drop") %>%
  left_join(
    selected_a %>%
      mutate(tract_area = as.numeric(st_area(.))) %>%
      st_drop_geometry() %>%
      select(GEOID, tract_area, COUNTYFP),
    by = "GEOID"
  ) %>%
  mutate(
    share_in_austin = int_area / tract_area,
    pct_in_austin   = share_in_austin * 100,
    pct_label       = sprintf("%.2f%%", pct_in_austin)
  )

selected_a <- selected_a %>%
  left_join(overlap_share %>% select(GEOID, COUNTYFP, share_in_austin, pct_in_austin, pct_label),
            by = c("GEOID", "COUNTYFP")) %>%
  mutate(
    overlap_class = ifelse(
      share_in_austin < tiny_threshold,
      paste0("Tiny overlap (<", tiny_threshold * 100, "%)"),
      paste0("Meaningful overlap (≥", tiny_threshold * 100, "%)")
    ),
    county_name = recode(
      COUNTYFP,
      "021" = "Bastrop",
      "209" = "Hays",
      "453" = "Travis",
      "491" = "Williamson",
      .default = paste0("CountyFP ", COUNTYFP)
    )
  )

#split layers
main_tracts <- selected_a %>% filter(share_in_austin >= tiny_threshold)
tiny_tracts <- selected_a %>% filter(share_in_austin <  tiny_threshold)
tiny_centroids <- st_centroid(tiny_tracts)


austin_boundary <- st_make_valid(austin_a)


austin_m        <- st_transform(austin_boundary, 3857)
main_m          <- st_transform(main_tracts, 3857)
tiny_m          <- st_transform(tiny_tracts, 3857)
tiny_cent_m     <- st_transform(tiny_centroids, 3857)

#build map object (tmap) 
m_toggle <-
  tm_basemap("OpenStreetMap") +

  #boundary
  tm_shape(austin_m, name = "Austin LTD + Full Purpose Boundary") +
  tm_borders(col = "black", lwd = 3) +

  #main tracts 
  tm_shape(main_m, name = "Selected tracts (≥ threshold overlap)") +
  tm_polygons(
    col = "county_name",
    fill_alpha = 0.25,
    popup.vars = c(
      "GEOID" = "GEOID",
      "County" = "county_name",
      "% of tract in Austin" = "pct_label",
      "Overlap class" = "overlap_class"
    )
  ) +
  tm_shape(main_m, name = "Selected tracts borders") +
  tm_borders(col = "dodgerblue3", lwd = 1) +

  #tiny overlap tracts 
  tm_shape(tiny_m, name = "Tiny-overlap tracts (< threshold)") +
  tm_polygons(
    col = "county_name",
    fill_alpha = 0.10,
    popup.vars = c(
      "GEOID" = "GEOID",
      "County" = "county_name",
      "% of tract in Austin" = "pct_label",
      "Overlap class" = "overlap_class"
    )
  ) +
  tm_shape(tiny_m, name = "Tiny-overlap borders") +
  tm_borders(col = "red3", lwd = 2, lty = "dashed") +

  #tiny-overlap centroids
  tm_shape(tiny_cent_m, name = "Tiny-overlap centroids") +
  tm_dots(
    col = "red3",
    size = 0.07,
    popup.vars = c(
      "GEOID" = "GEOID",
      "County" = "county_name",
      "% of tract in Austin" = "pct_label"
    )
  ) +

  tm_title(
    paste0(
      "Austin LTD + Full Purpose: Selected Census Tracts\n",
      "Tiny overlap threshold = ", tiny_threshold * 100, "% of tract area"
    )
  )

#convert to leaflet widgete 
w_toggle <- fix_widget(tmap_leaflet(m_toggle))
w_toggle

library(sf)
library(dplyr)
library(tmap)
library(htmlwidgets)

tmap_mode("view")

fix_widget <- function(w) {
  htmlwidgets::onRender(
    w,
    "function(el, x) {
       setTimeout(function(){ window.dispatchEvent(new Event('resize')); }, 250);
       setTimeout(function(){ window.dispatchEvent(new Event('resize')); }, 1000);
     }"
  )
}

austin_a   <- st_transform(austin_union, 3083)
tracts_a   <- st_transform(selected_tracts, 3083)

inside_pieces <- st_intersection(
  tracts_a %>% select(GEOID, COUNTYFP),
  st_sf(geometry = austin_a)
)

share_tbl <- inside_pieces %>%
  mutate(piece_area = as.numeric(st_area(.))) %>%
  st_drop_geometry() %>%
  group_by(GEOID) %>%
  summarise(in_austin_area = sum(piece_area), .groups = "drop") %>%
  left_join(
    tracts_a %>%
      mutate(tract_area = as.numeric(st_area(.))) %>%
      st_drop_geometry() %>%
      select(GEOID, tract_area),
    by = "GEOID"
  ) %>%
  mutate(pct_in_austin = 100 * in_austin_area / tract_area)

tracts_a <- tracts_a %>%
  left_join(share_tbl %>% select(GEOID, pct_in_austin), by = "GEOID") %>%
  mutate(split_flag = ifelse(pct_in_austin > 5 & pct_in_austin < 95, "Split tract", "Mostly in/out"))

inside_pieces <- inside_pieces %>%
  left_join(share_tbl %>% select(GEOID, pct_in_austin), by = "GEOID")

austin_m <- st_transform(austin_a, 3857)
tracts_m <- st_transform(tracts_a, 3857)
inside_m <- st_transform(inside_pieces, 3857)

m_overlap <-
  tm_basemap("OpenStreetMap") +
  tm_shape(austin_m, name = "Austin LTD + Full Purpose boundary") +
  tm_borders(col = "black", lwd = 3) +
  tm_shape(tracts_m, name = "Full selected tracts") +
  tm_borders(col = "dodgerblue3", lwd = 1) +
  tm_shape(inside_m, name = "Portion of tract inside Austin") +
  tm_polygons(
    col = "pct_in_austin",
    fill_alpha = 0.45,
    popup.vars = c(
      "GEOID" = "GEOID",
      "% of tract in Austin" = "pct_in_austin"
    )
  ) +
  tm_shape(tracts_m %>% filter(split_flag == "Split tract"), name = "Split tracts (pop out)") +
  tm_borders(col = "red3", lwd = 3) +
  tm_title("Austin LTD + Full Purpose: Where each tract overlaps the boundary")

w_overlap <- fix_widget(tmap_leaflet(m_overlap))
w_overlap

5 PERC IN ATX

library(sf)
library(dplyr)

# Equal-area CRS for Texas for accurate area calculations
austin_a   <- st_transform(austin_union, 3083)
selected_a <- st_transform(selected_tracts, 3083)

# Intersection pieces (portion of each selected tract inside Austin)
inside_pieces <- st_intersection(
  selected_a %>% select(GEOID, COUNTYFP),
  st_sf(geometry = austin_a)
)

# Percent of tract area inside Austin
share_tbl <- inside_pieces %>%
  mutate(piece_area = as.numeric(st_area(.))) %>%
  st_drop_geometry() %>%
  group_by(GEOID) %>%
  summarise(in_austin_area = sum(piece_area), .groups = "drop") %>%
  left_join(
    selected_a %>%
      mutate(tract_area = as.numeric(st_area(.))) %>%
      st_drop_geometry() %>%
      select(GEOID, tract_area, COUNTYFP),
    by = "GEOID"
  ) %>%
  mutate(pct_in_austin = 100 * in_austin_area / tract_area)

# Attach pct to selected tracts
selected_a <- selected_a %>%
  left_join(share_tbl %>% select(GEOID, pct_in_austin), by = "GEOID")


library(ggplot2)

df_share <- selected_a %>% st_drop_geometry()

ggplot(df_share, aes(x = pct_in_austin)) +
  geom_histogram(bins = 30) +
  labs(
    title = "How much of each selected tract lies inside Austin LTD + Full Purpose",
    x = "% of tract area inside Austin boundary",
    y = "Number of tracts"
  ) +
  theme_minimal()

MAP of Tracts split by % inside ATX jurisdiction

library(tmap)
tmap_mode("view")

split_tracts <- selected_a %>%
  filter(pct_in_austin > 5, pct_in_austin < 95)

tm_shape(st_transform(austin_a, 3857)) +
  tm_borders(col = "black", lwd = 3) +
  tm_shape(st_transform(split_tracts, 3857)) +
  tm_polygons(
    col = "pct_in_austin",
    alpha = 0.55,
    border.col = "red3",
    lwd = 2,
    popup.vars = c("GEOID"="GEOID", "% in Austin"="pct_in_austin"),
    title = "% of tract inside Austin"
  ) +
  tm_layout(
    title = "Split tracts (5%–95% of tract area inside Austin)",
    legend.outside = TRUE,
    frame = FALSE
  )

df_share %>%
  count(COUNTYFP, name = "n_tracts") %>%
  ggplot(aes(x = COUNTYFP, y = n_tracts)) +
  geom_col() +
  labs(
    title = "Selected tracts by county (Austin LTD + Full Purpose overlap)",
    x = "County FIPS (COUNTYFP)",
    y = "Number of tracts"
  ) +
  theme_minimal()

library(dplyr)

# Centroid-within set (you already have this in your workflow, but here it is again)
centroid_within <- tx_tracts[
  st_within(st_centroid(tx_tracts), austin_union, sparse = FALSE),
]

# Any overlap set (you already have selected_tracts = 271)
any_overlap <- selected_tracts

# Area-threshold sets using pct_in_austin
area_10 <- selected_a %>% filter(pct_in_austin >= 10)
area_25 <- selected_a %>% filter(pct_in_austin >= 25)
area_50 <- selected_a %>% filter(pct_in_austin >= 50)

sensitivity <- tibble::tibble(
  rule = c("Any overlap (area>0)", "Centroid within", "≥10% area in Austin", "≥25% area in Austin", "≥50% area in Austin"),
  n_tracts = c(nrow(any_overlap), nrow(centroid_within), nrow(area_10), nrow(area_25), nrow(area_50))
)

sensitivity

# A tibble: 5 × 2
  rule                 n_tracts
  <chr>                   <int>
1 Any overlap (area>0)      271
2 Centroid within           218
3 ≥10% area in Austin       250
4 ≥25% area in Austin       234
5 ≥50% area in Austin       217

WE need to look at tracts carefully and decide which to include in the index. We can do two things: 1)use census blocks to weight population inside Austin or 2) use % of tract area inside Austin as proxy

Here I try the first option:

library(dplyr)

tract_classification <- selected_a %>%
  st_drop_geometry() %>%
  mutate(
    inclusion_class = case_when(
      pct_in_austin >= 25 ~ "Included (>=25% in Austin)",
      pct_in_austin >= 10 ~ "Borderline (10–25% in Austin)",
      pct_in_austin < 10  ~ "Excluded (<10% in Austin)"
    )
  )

tract_classification %>%
  count(inclusion_class)

                inclusion_class   n
1 Borderline (10–25% in Austin)  16
2     Excluded (<10% in Austin)  21
3    Included (>=25% in Austin) 234

#this is the table we can use in appendix 
excluded_tracts_table <- tract_classification %>%
  filter(inclusion_class == "Excluded (<10% in Austin)") %>%
  arrange(pct_in_austin) %>%
  select(
    GEOID,
    COUNTYFP,
    pct_in_austin,
    inclusion_class
  )

excluded_tracts_table

         GEOID COUNTYFP pct_in_austin           inclusion_class
1  48021950303      021   0.004450826 Excluded (<10% in Austin)
2  48491020347      491   0.020754077 Excluded (<10% in Austin)
3  48453002217      453   0.027582693 Excluded (<10% in Austin)
4  48453036600      453   0.368559843 Excluded (<10% in Austin)
5  48453046400      453   0.405443205 Excluded (<10% in Austin)
6  48453001916      453   0.511745425 Excluded (<10% in Austin)
7  48491020333      491   0.512965327 Excluded (<10% in Austin)
8  48453035600      453   0.542336500 Excluded (<10% in Austin)
9  48453045600      453   0.811250725 Excluded (<10% in Austin)
10 48453031400      453   0.871103424 Excluded (<10% in Austin)
11 48209010923      209   1.019364783 Excluded (<10% in Austin)
12 48453044700      453   1.171145237 Excluded (<10% in Austin)
13 48453031100      453   1.195452089 Excluded (<10% in Austin)
14 48453035700      453   1.199655826 Excluded (<10% in Austin)
15 48453001918      453   1.717760603 Excluded (<10% in Austin)
16 48453002434      453   3.964166748 Excluded (<10% in Austin)
17 48453002216      453   4.489320305 Excluded (<10% in Austin)
18 48491020356      491   6.649132617 Excluded (<10% in Austin)
19 48453002215      453   8.517314922 Excluded (<10% in Austin)
20 48453002221      453   8.710728051 Excluded (<10% in Austin)
21 48453031500      453   8.756173990 Excluded (<10% in Austin)

write.csv(
  excluded_tracts_table,
  "C:/Users/rayo-garza/OneDrive - Center for Public Policy Priorities/Documents/Data 2025/Excluded_Tracts_Less_Than_10pct_Austin.csv",
  row.names = FALSE
)


excluded_tracts_table %>%
  count(COUNTYFP, name = "n_excluded_tracts")

  COUNTYFP n_excluded_tracts
1      021                 1
2      209                 1
3      453                16
4      491                 3

tract_classification %>%
  group_by(inclusion_class) %>%
  summarise(
    n_tracts = n(),
    mean_pct_in_austin = mean(pct_in_austin),
    min_pct = min(pct_in_austin),
    max_pct = max(pct_in_austin)
  )

# A tibble: 3 × 5
  inclusion_class               n_tracts mean_pct_in_austin  min_pct max_pct
  <chr>                            <int>              <dbl>    <dbl>   <dbl>
1 Borderline (10–25% in Austin)       16              17.9  10.4       24.1 
2 Excluded (<10% in Austin)           21               2.45  0.00445    8.76
3 Included (>=25% in Austin)         234              91.1  26.3      100.  

Mapping Tracts Likely to Exclude based on cut-off’s

library(sf)
library(dplyr)
library(tmap)
library(htmlwidgets)

tmap_mode("view")

fix_widget <- function(w) {
  htmlwidgets::onRender(
    w,
    "function(el, x) {
       setTimeout(function(){ window.dispatchEvent(new Event('resize')); }, 250);
       setTimeout(function(){ window.dispatchEvent(new Event('resize')); }, 1000);
     }"
  )
}

excluded_ids <- tract_classification %>%
  filter(inclusion_class == "Excluded (<10% in Austin)") %>%
  pull(GEOID)

excluded_tracts <- selected_a %>%
  filter(GEOID %in% excluded_ids)

name_field <- dplyr::case_when(
  "NAMELSAD" %in% names(excluded_tracts) ~ "NAMELSAD",
  "NAME" %in% names(excluded_tracts) ~ "NAME",
  TRUE ~ NA_character_
)

excluded_tracts <- excluded_tracts %>%
  mutate(
    tract_name = if (!is.na(name_field)) .data[[name_field]] else GEOID,
    county_name = recode(
      COUNTYFP,
      "021" = "Bastrop",
      "209" = "Hays",
      "453" = "Travis",
      "491" = "Williamson",
      .default = paste0("CountyFP ", COUNTYFP)
    )
  ) %>%
  select(GEOID, tract_name, COUNTYFP, county_name, pct_in_austin)

excluded_inside <- inside_pieces %>%
  filter(GEOID %in% excluded_ids) %>%
  left_join(
    excluded_tracts %>%
      st_drop_geometry() %>%
      select(GEOID, COUNTYFP, tract_name, county_name, pct_in_austin),
    by = c("GEOID", "COUNTYFP")
  ) %>%
  select(GEOID, tract_name, COUNTYFP, county_name, pct_in_austin)

austin_m      <- st_transform(austin_union, 3857)
excluded_m    <- st_transform(excluded_tracts, 3857)
excluded_in_m <- st_transform(excluded_inside, 3857)

m_excluded <-
  tm_basemap("OpenStreetMap") +
  tm_shape(austin_m, name = "Austin boundary") +
  tm_fill(col = "grey75", alpha = 0.4) +
  tm_borders(col = "black", lwd = 3) +
  tm_shape(excluded_m, name = "Excluded tracts (full)") +
  tm_borders(col = "blue", lwd = 3.5) +
  tm_fill(alpha = 0, popup.vars = c(
    "Tract name" = "tract_name",
    "County" = "county_name",
    "GEOID" = "GEOID",
    "% in Austin" = "pct_in_austin"
  )) +
  tm_shape(excluded_in_m, name = "Portion of tract inside Austin") +
  tm_fill(col = "red3", alpha = 0.7, popup.vars = c(
    "Tract name" = "tract_name",
    "County" = "county_name",
    "GEOID" = "GEOID",
    "% of tract in Austin" = "pct_in_austin"
  )) +
  tm_borders(col = "red4", lwd = 1) +
  tm_title("Excluded Census Tracts: Only Small Portions Lie Inside Austin")

w_excluded <- fix_widget(tmap_leaflet(m_excluded))
w_excluded

kept_tracts <- tract_classification %>%
  filter(inclusion_class != "Excluded (<10% in Austin)")

nrow(kept_tracts)

[1] 250

After excluding census tracts with less than 10 percent of their land area inside Austin’s Full or Limited Purpose jurisdiction, the final analysis universe includes 250 census tracts

KAITLAN BELOW

# head(austin_jur)

Here I am filtering to include only limited and full-purpose jurisdictions in Austin.

# austin_jur_filtered <- austin_jur %>%
#   filter(jurisdic_2 %in% c("AUSTIN LTD", "AUSTIN FULL PURPOSE"))
# 
# print(unique(austin_jur_filtered$jurisdic_2))

# # read the Texas Census tracts shapefile
# tx_census_tracts <- st_read("tl_2025_48_tract/tl_2025_48_tract.shp")
# 
# # check columns
# print(names(tx_census_tracts))

# # rebuild a clean GEOID (I was getting mostly NA values for some reason for GEOID)
# tx_census_tracts <- tx_census_tracts %>%
#   mutate(
#     GEOID_clean = paste0(STATEFP, COUNTYFP, TRACTCE)
#   )

# check
# head(tx_census_tracts$GEOID_clean)

# # double check
# #unique(census_tracts$GEOID_clean)
# 
# # inspect new GEOIDs in Excel
# tract_ids <- unique(tx_census_tracts$GEOID_clean)
# 
# write.csv(tract_ids, "All_Census_Tract_GEOIDs.csv", row.names = FALSE)

# there are 6,896 tracts in Texas according to this list

# # make sure both layers have valid geometry
# austin_jur_filtered <- st_make_valid(austin_jur_filtered)
# tx_census_tracts <- st_make_valid(tx_census_tracts)

# # match coordinate reference systems
# if (st_crs(tx_census_tracts) != st_crs(austin_jur_filtered)) {
#   tx_census_tracts <- st_transform(tx_census_tracts, st_crs(austin_jur_filtered))
# }

6 Spatial Join

# # spatially join tracts to any Austin jurisdiction polygon they intersect
# selected_tracts <- st_join(tx_census_tracts, austin_jur_filtered, join = st_intersects)

# # keep only those that actually intersect Austin (non-NA jurisdiction label)
# selected_tracts <- selected_tracts %>%
#   filter(!is.na(jurisdic_2))

7 Check the Work

# check for NA values
# table(is.na(selected_tracts$GEOID_clean))

# head(selected_tracts$GEOID_clean, 20)

# # quick check that the selected tracts look reasonable
# print(paste("# of tracts intersecting Austin:", nrow(selected_tracts)))
# print(paste("# of unique GEOIDs:", length(unique(selected_tracts$GEOID_clean))))

There are 271 intersecting census tracts, meaning there are 271 tracts in Austin full and limited purpose jurisdictions. This makes sense, as Travis County alone has 291 census tracts, though not all of Travis County is in Austin. There are a total of 493 census tracts in all of Bastrop, Hays, Travis, and Williamson Counties combined.

The reason there are 383 rows, but only 271 distinct tracts, could be because some tracts may touch multiple Austin jurisdiction polygons since Austin’s boundary is split into parts.

# optional quick view of unique GEOIDs
# print(sort(unique(selected_tracts$GEOID_clean))[1:20])

# # to save only distinct GEOIDs with no duplicates
# selected_unique <- selected_tracts %>%
#   distinct(GEOID_clean, .keep_all = TRUE)

8 Save Final Files

# # save the intersected tracts shapefile
# st_write(selected_tracts, "Austin_LTD_FULL_Tracts.shp")
# 
# # Also export a CSV list of tract GEOIDs and names
# tract_list <- selected_tracts %>%
#   st_drop_geometry() %>%
#   arrange(GEOID_clean)
# 
# write_csv(tract_list, "Austin_LTD_FULL_Tract_List.csv")
# 
# ###################################################################################
# 
# # files without duplicate GEOIDs
# st_write(selected_unique, "Austin_LTD_FULL_Tracts_Unique.shp")
# 
# unique_tract_list <- selected_unique %>%
#   st_drop_geometry() %>%
#   arrange(GEOID_clean)
# 
# write_csv(unique_tract_list, "Austin_LTD_FULL_Tract_List_Unique.csv")