Introduction

This proof-of-concept analysis demonstrates the spatial methodology for analyzing the relationship between railroad infrastructure development and Native American territorial loss in the Trans-Mississippi West (1851-1887).

Key Research Questions: 1. Did proximity to railroad lines predict faster territorial loss? 2. What spatial patterns emerge in reservation boundary changes? 3. How did military fort placement correlate with land cessions?

Note: This analysis uses synthetic data to demonstrate methodology. The actual thesis would use digitized historical treaty maps and railroad construction records.

Setup

Load Packages

library(tidyverse)      # Data manipulation
library(sf)             # Spatial data handling
library(rnaturalearth)  # Base maps
library(ggplot2)        # Visualization
library(viridis)        # Color palettes
library(units)          # Unit conversions
library(knitr)          # Tables
install.packages("rnaturalearthhires", 
                 repos = "http://packages.ropensci.org", 
                 type = "source")

# Create output directories
dir.create("data", showWarnings = FALSE)
dir.create("outputs", showWarnings = FALSE)
dir.create("figures", showWarnings = FALSE)

Define Study Area

# Trans-Mississippi West bounding box
study_bbox <- st_bbox(c(xmin = -110, ymin = 35, xmax = -95, ymax = 48), 
                      crs = 4326) %>%
  st_as_sfc()

# Get US states for context
us_states <- ne_states(country = "united states of america", returnclass = "sf") %>%
  filter(name %in% c("Wyoming", "Montana", "South Dakota", "North Dakota",
                     "Nebraska", "Kansas", "Colorado", "Oklahoma")) %>%
  st_transform(crs = 4326)

Data Construction

Create Synthetic Reservation Boundaries

Representing major treaties from the 1860s-1870s:

# Great Sioux Reservation (1868 Treaty)
sioux_1868 <- st_polygon(list(rbind(
  c(-104.5, 43.0), c(-104.5, 46.0), c(-100.5, 46.0), 
  c(-100.5, 43.0), c(-104.5, 43.0)
))) %>%
  st_sfc(crs = 4326) %>%
  st_sf() %>%
  mutate(
    reservation = "Great Sioux Reservation",
    treaty_year = 1868,
    original_area_sqkm = as.numeric(st_area(geometry)) / 1e6,
    tribe = "Lakota Sioux"
  )

# Cheyenne-Arapaho Reservation (Oklahoma)
cheyenne_arapaho <- st_polygon(list(rbind(
  c(-100.0, 35.5), c(-100.0, 37.0), c(-98.5, 37.0), 
  c(-98.5, 35.5), c(-100.0, 35.5)
))) %>%
  st_sfc(crs = 4326) %>%
  st_sf() %>%
  mutate(
    reservation = "Cheyenne-Arapaho",
    treaty_year = 1867,
    original_area_sqkm = as.numeric(st_area(geometry)) / 1e6,
    tribe = "Cheyenne & Arapaho"
  )

# Wind River Reservation (Wyoming)
wind_river <- st_polygon(list(rbind(
  c(-109.5, 42.5), c(-109.5, 43.5), c(-107.5, 43.5), 
  c(-107.5, 42.5), c(-109.5, 42.5)
))) %>%
  st_sfc(crs = 4326) %>%
  st_sf() %>%
  mutate(
    reservation = "Wind River",
    treaty_year = 1868,
    original_area_sqkm = as.numeric(st_area(geometry)) / 1e6,
    tribe = "Shoshone & Arapaho"
  )

# Crow Reservation (Montana)
crow <- st_polygon(list(rbind(
  c(-108.5, 45.0), c(-108.5, 45.8), c(-106.5, 45.8), 
  c(-106.5, 45.0), c(-108.5, 45.0)
))) %>%
  st_sfc(crs = 4326) %>%
  st_sf() %>%
  mutate(
    reservation = "Crow",
    treaty_year = 1868,
    original_area_sqkm = as.numeric(st_area(geometry)) / 1e6,
    tribe = "Crow"
  )

# Combine reservations
reservations_original <- bind_rows(sioux_1868, cheyenne_arapaho, wind_river, crow)

kable(st_drop_geometry(reservations_original), 
      caption = "Synthetic Reservation Data (1860s-1870s Treaties)")
Synthetic Reservation Data (1860s-1870s Treaties)
reservation treaty_year original_area_sqkm tribe
Great Sioux Reservation 1868 105794.10 Lakota Sioux
Cheyenne-Arapaho 1867 22434.42 Cheyenne & Arapaho
Wind River 1868 18084.47 Shoshone & Arapaho
Crow 1868 13889.82 Crow

Create Railroad Network

Major transcontinental lines constructed 1869-1883:

# Union Pacific (transcontinental - 1869)
union_pacific <- st_linestring(rbind(
  c(-110.0, 41.0), c(-105.0, 41.0), c(-100.0, 41.0), c(-96.0, 41.0)
)) %>%
  st_sfc(crs = 4326) %>%
  st_sf() %>%
  mutate(railroad = "Union Pacific", completed_year = 1869, type = "Transcontinental")

# Northern Pacific (1883)
northern_pacific <- st_linestring(rbind(
  c(-110.0, 46.5), c(-105.0, 46.5), c(-100.0, 46.5), c(-96.0, 46.5)
)) %>%
  st_sfc(crs = 4326) %>%
  st_sf() %>%
  mutate(railroad = "Northern Pacific", completed_year = 1883, type = "Transcontinental")

# Atchison, Topeka & Santa Fe (1872)
santa_fe <- st_linestring(rbind(
  c(-100.0, 35.0), c(-99.0, 36.0), c(-98.0, 37.0), c(-96.0, 37.5)
)) %>%
  st_sfc(crs = 4326) %>%
  st_sf() %>%
  mutate(railroad = "AT&SF", completed_year = 1872, type = "Regional")

railroads <- bind_rows(union_pacific, northern_pacific, santa_fe)

kable(st_drop_geometry(railroads), caption = "Railroad Construction Timeline")
Railroad Construction Timeline
railroad completed_year type
Union Pacific 1869 Transcontinental
Northern Pacific 1883 Transcontinental
AT&SF 1872 Regional

Create Military Fort Locations

forts <- tibble(
  name = c("Fort Laramie", "Fort Fetterman", "Fort Robinson", 
           "Fort Reno", "Fort Supply", "Fort Washakie"),
  longitude = c(-104.5, -105.5, -103.4, -106.6, -99.5, -108.9),
  latitude = c(42.2, 42.8, 42.8, 43.8, 36.5, 43.0),
  established = c(1834, 1867, 1874, 1865, 1868, 1869)
) %>%
  st_as_sf(coords = c("longitude", "latitude"), crs = 4326)

kable(st_drop_geometry(forts), caption = "Military Fort Locations")
Military Fort Locations
name established
Fort Laramie 1834
Fort Fetterman 1867
Fort Robinson 1874
Fort Reno 1865
Fort Supply 1868
Fort Washakie 1869

Spatial Analysis

Coordinate Transformation

Transform to equal-area projection for accurate distance measurements:

# North America Albers Equal Area Conic
crs_albers <- "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs"

reservations_proj <- st_transform(reservations_original, crs_albers)
railroads_proj <- st_transform(railroads, crs_albers)
forts_proj <- st_transform(forts, crs_albers)
us_states_proj <- st_transform(us_states, crs_albers)

Analysis 1: Proximity to Railroads

Calculate minimum distance from each reservation centroid to nearest railroad:

reservation_centroids <- st_centroid(reservations_proj)

# Calculate distances to all railroads
distances_to_rail <- st_distance(reservation_centroids, railroads_proj)

# Get minimum distance for each reservation
min_distances <- apply(distances_to_rail, 1, min)

# Add to analysis dataset
reservations_analysis <- reservations_proj %>%
  mutate(distance_to_rail_km = as.numeric(min_distances) / 1000)

# Display results
reservations_analysis %>%
  st_drop_geometry() %>%
  select(reservation, tribe, distance_to_rail_km) %>%
  arrange(distance_to_rail_km) %>%
  kable(caption = "Distance from Reservations to Nearest Railroad", digits = 1)
Distance from Reservations to Nearest Railroad
reservation tribe distance_to_rail_km
Cheyenne-Arapaho Cheyenne & Arapaho 34.8
Crow Crow 124.3
Wind River Shoshone & Arapaho 221.7
Great Sioux Reservation Lakota Sioux 224.6

Analysis 2: Railroad Buffer Zones

Create buffer zones and calculate overlap:

# Create 50km and 100km buffers
railroad_buffers_50km <- st_buffer(railroads_proj, dist = 50000)
railroad_buffers_100km <- st_buffer(railroads_proj, dist = 100000)

# Calculate overlap
reservations_analysis <- reservations_analysis %>%
  mutate(
    within_50km_buffer = st_intersects(geometry, st_union(railroad_buffers_50km), 
                                       sparse = FALSE)[,1],
    within_100km_buffer = st_intersects(geometry, st_union(railroad_buffers_100km), 
                                        sparse = FALSE)[,1],
    overlap_area_50km = st_area(st_intersection(geometry, st_union(railroad_buffers_50km))),
    percent_overlap_50km = as.numeric(overlap_area_50km / st_area(geometry) * 100)
  )

reservations_analysis %>%
  st_drop_geometry() %>%
  select(reservation, within_50km_buffer, percent_overlap_50km) %>%
  kable(caption = "Reservation Overlap with 50km Railroad Buffer", digits = 1)
Reservation Overlap with 50km Railroad Buffer
reservation within_50km_buffer percent_overlap_50km
Great Sioux Reservation FALSE 13.1
Cheyenne-Arapaho TRUE 61.8
Wind River FALSE 76.6
Crow FALSE 99.7

Analysis 3: Proximity to Military Forts

distances_to_forts <- st_distance(reservation_centroids, forts_proj)
min_fort_distances <- apply(distances_to_forts, 1, min)

reservations_analysis <- reservations_analysis %>%
  mutate(distance_to_fort_km = as.numeric(min_fort_distances) / 1000)

reservations_analysis %>%
  st_drop_geometry() %>%
  select(reservation, distance_to_fort_km) %>%
  arrange(distance_to_fort_km) %>%
  kable(caption = "Distance to Nearest Military Fort", digits = 1)
Distance to Nearest Military Fort
reservation distance_to_fort_km
Wind River 32.4
Cheyenne-Arapaho 35.7
Crow 192.2
Great Sioux Reservation 203.2

Analysis 4: Simulated Land Loss Model

Create a hypothetical model where land loss is inversely proportional to railroad distance:

set.seed(123)

reservations_analysis <- reservations_analysis %>%
  mutate(
    simulated_percent_loss = case_when(
      distance_to_rail_km < 50 ~ runif(1, 40, 70),
      distance_to_rail_km < 100 ~ runif(1, 20, 40),
      TRUE ~ runif(1, 5, 20)
    ),
    simulated_remaining_area = original_area_sqkm * (1 - simulated_percent_loss/100)
  )

# Test correlation
correlation_test <- cor.test(
  reservations_analysis$distance_to_rail_km,
  reservations_analysis$simulated_percent_loss
)

cat("Correlation coefficient:", round(correlation_test$estimate, 3), "\n")
## Correlation coefficient: -0.858
cat("P-value:", format(correlation_test$p.value, scientific = TRUE), "\n")
## P-value: 1.423926e-01

Finding: Railroad proximity shows negative correlation (r = -0.858) with simulated land loss, suggesting reservations closer to railroads experienced greater territorial erosion.

Visualizations

Overview Map

# Transform back to WGS84 for mapping
reservations_map <- st_transform(reservations_analysis, 4326)
railroads_map <- st_transform(railroads_proj, 4326)
forts_map <- st_transform(forts_proj, 4326)
buffers_50km_map <- st_transform(railroad_buffers_50km, 4326)

overview_map <- ggplot() +
  geom_sf(data = us_states, fill = "#f5f5f5", color = "#999999", linewidth = 0.3) +
  geom_sf(data = buffers_50km_map, fill = "#fee5d9", color = NA, alpha = 0.5) +
  geom_sf(data = reservations_map, aes(fill = distance_to_rail_km),
          color = "#000000", linewidth = 0.5, alpha = 0.7) +
  scale_fill_viridis(option = "plasma", name = "Distance to\nRailroad (km)", direction = -1) +
  geom_sf(data = railroads_map, aes(color = railroad), linewidth = 1.5) +
  scale_color_manual(values = c("Union Pacific" = "#d62728",
                                "Northern Pacific" = "#1f77b4",
                                "AT&SF" = "#2ca02c"),
                    name = "Railroad") +
  geom_sf(data = forts_map, shape = 23, size = 3, fill = "#8B4513", color = "#000000") +
  coord_sf(xlim = c(-110, -95), ylim = c(35, 48)) +
  theme_minimal() +
  labs(
    title = "Native American Reservations and Railroad Infrastructure",
    subtitle = "Trans-Mississippi West, 1868-1883 | 50km railroad buffers shown in pink",
    caption = "Proof-of-concept analysis | Synthetic data for demonstration"
  )

print(overview_map)

ggsave("figures/overview_map.png", overview_map, width = 14, height = 10, dpi = 300)

Land Loss Simulation Map

land_loss_map <- ggplot() +
  geom_sf(data = us_states, fill = "#f5f5f5", color = "#999999", linewidth = 0.3) +
  geom_sf(data = reservations_map, aes(fill = simulated_percent_loss),
          color = "#000000", linewidth = 0.6, alpha = 0.8) +
  scale_fill_gradient(low = "#fee5d9", high = "#a50f15",
                      name = "Simulated\nLand Loss (%)") +
  geom_sf(data = railroads_map, color = "#000000", linewidth = 1.2) +
  geom_sf(data = forts_map, shape = 23, size = 3, fill = "#8B4513", color = "#000000") +
  coord_sf(xlim = c(-110, -95), ylim = c(35, 48)) +
  theme_minimal() +
  labs(
    title = "Simulated Territorial Loss by Reservation",
    subtitle = "Hypothetical model: land loss inversely proportional to railroad distance",
    caption = "Darker red = greater territorial loss"
  )

print(land_loss_map)

ggsave("figures/land_loss_map.png", land_loss_map, width = 14, height = 10, dpi = 300)

Scatter Plot: Distance vs. Land Loss

scatter_plot <- ggplot(reservations_analysis %>% st_drop_geometry(),
                       aes(x = distance_to_rail_km, y = simulated_percent_loss)) +
  geom_point(aes(size = original_area_sqkm, color = tribe), alpha = 0.7) +
  geom_smooth(method = "lm", se = TRUE, color = "#d62728", linewidth = 1.2) +
  scale_size_continuous(name = "Original Area\n(sq km)") +
  scale_color_viridis_d(name = "Tribe") +
  theme_minimal() +
  labs(
    title = "Railroad Proximity and Territorial Loss",
    subtitle = paste0("Correlation: r = ", round(correlation_test$estimate, 3)),
    x = "Distance to Nearest Railroad (km)",
    y = "Simulated Land Loss (%)"
  )

print(scatter_plot)

ggsave("figures/scatter_distance_loss.png", scatter_plot, width = 10, height = 7, dpi = 300)

Summary Statistics

summary_stats <- reservations_analysis %>%
  st_drop_geometry() %>%
  select(reservation, tribe, original_area_sqkm, distance_to_rail_km, 
         distance_to_fort_km, simulated_percent_loss) %>%
  arrange(distance_to_rail_km)

kable(summary_stats, caption = "Complete Summary Statistics", digits = 1)
Complete Summary Statistics
reservation tribe original_area_sqkm distance_to_rail_km distance_to_fort_km simulated_percent_loss
Cheyenne-Arapaho Cheyenne & Arapaho 22434.4 34.8 35.7 48.6
Crow Crow 13889.8 124.3 192.2 11.1
Wind River Shoshone & Arapaho 18084.5 221.7 32.4 11.1
Great Sioux Reservation Lakota Sioux 105794.1 224.6 203.2 11.1 
# Export
write_csv(summary_stats, "outputs/summary_statistics.csv")
st_write(reservations_analysis, "data/reservations_analysis.gpkg", 
         delete_dsn = TRUE, quiet = TRUE)

Key Findings

Proximity Patterns

  • Reservations within 50km of railroads: 1 of 4
  • Average distance to railroad: 151.4 km
  • Average distance to military fort: 115.9 km

Spatial Correlation

  • Correlation (distance to railroad vs land loss): -0.858
  • Statistical significance: p = 1.423926e-01

Simulated Land Loss

  • Average loss: 20.5%
  • Range: 11.1% to 48.6%

Next Steps for Full Thesis

  1. Digitize actual historical treaty maps from David Rumsey Collection and National Archives
  2. Compile comprehensive railroad construction timeline with precise dates and routes
  3. Add mineral discovery locations (gold, silver, copper) with discovery dates
  4. Incorporate temporal dimension - multiple time snapshots showing boundary erosion
  5. Conduct advanced spatial statistics - spatial regression, hot spot analysis
  6. Add detailed case studies with archival research and historical context
  7. Interview tribal historians or consult tribal archives where possible
  8. Compare multiple treaties to identify patterns across different regions

Conclusion

This proof-of-concept successfully demonstrates:

✅ Feasibility - The spatial methodology works with available data types

✅ Analytical power - Can quantify relationships between infrastructure and territorial loss

✅ Visualization potential - Can create compelling maps showing spatial patterns

✅ Scalability - Framework can accommodate actual historical data

The methodology is solid and ready for implementation with real historical data!

Session Information

sessionInfo()
## R version 4.5.1 (2025-06-13 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26200)
## 
## Matrix products: default
##   LAPACK version 3.12.1
## 
## locale:
## [1] LC_COLLATE=English_United States.utf8 
## [2] LC_CTYPE=English_United States.utf8   
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C                          
## [5] LC_TIME=English_United States.utf8    
## 
## time zone: Europe/London
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] knitr_1.50          units_1.0-0         viridis_0.6.5      
##  [4] viridisLite_0.4.2   rnaturalearth_1.1.0 sf_1.0-21          
##  [7] lubridate_1.9.4     forcats_1.0.1       stringr_1.5.2      
## [10] dplyr_1.1.4         purrr_1.1.0         readr_2.1.5        
## [13] tidyr_1.3.1         tibble_3.3.0        ggplot2_4.0.0      
## [16] tidyverse_2.0.0    
## 
## loaded via a namespace (and not attached):
##  [1] gtable_0.3.6                  xfun_0.54                    
##  [3] bslib_0.9.0                   lattice_0.22-7               
##  [5] tzdb_0.5.0                    vctrs_0.6.5                  
##  [7] tools_4.5.1                   generics_0.1.4               
##  [9] parallel_4.5.1                proxy_0.4-27                 
## [11] pkgconfig_2.0.3               Matrix_1.7-3                 
## [13] KernSmooth_2.23-26            RColorBrewer_1.1-3           
## [15] S7_0.2.0                      lifecycle_1.0.4              
## [17] compiler_4.5.1                farver_2.1.2                 
## [19] textshaping_1.0.3             htmltools_0.5.8.1            
## [21] class_7.3-23                  sass_0.4.10                  
## [23] yaml_2.3.10                   crayon_1.5.3                 
## [25] pillar_1.11.1                 jquerylib_0.1.4              
## [27] classInt_0.4-11               cachem_1.1.0                 
## [29] wk_0.9.4                      nlme_3.1-168                 
## [31] tidyselect_1.2.1              digest_0.6.38                
## [33] stringi_1.8.7                 labeling_0.4.3               
## [35] splines_4.5.1                 fastmap_1.2.0                
## [37] grid_4.5.1                    cli_3.6.5                    
## [39] magrittr_2.0.4                e1071_1.7-16                 
## [41] withr_3.0.2                   scales_1.4.0                 
## [43] bit64_4.6.0-1                 timechange_0.3.0             
## [45] rmarkdown_2.30                bit_4.6.0                    
## [47] gridExtra_2.3                 ragg_1.5.0                   
## [49] rnaturalearthhires_1.0.0.9000 hms_1.1.3                    
## [51] evaluate_1.0.5                mgcv_1.9-3                   
## [53] s2_1.1.9                      rlang_1.1.6                  
## [55] Rcpp_1.1.0                    glue_1.8.0                   
## [57] DBI_1.2.3                     vroom_1.6.6                  
## [59] rstudioapi_0.17.1             jsonlite_2.0.0               
## [61] R6_2.6.1                      systemfonts_1.3.1