1 Set up the environment

This section establishes the computational environment for processing socioeconomic data inputs for the Lower Savannah Council of Governments regional travel demand model using both R and Python platforms.

1.1 Install and load packages

The package installation process incorporates essential libraries for comprehensive geospatial data analysis. The R environment includes tidyverse for data manipulation, sf for spatial data handling, tidycensus for Census Bureau data access, and lehdr for Longitudinal Employer-Household Dynamics data retrieval. The Python environment focuses on core data science libraries including {pandas} for data manipulation, {geopandas} for spatial analysis, and {pygris} for Census data queries. These packages form the analytical backbone for processing demographic, employment, and geographic data required for travel demand modeling.

R
Python

# # Install required packages
# install.packages(c(
#   "tidyverse",
#   "vroom",
#   "sf",
#   "tidycensus",
#   "lehdr",
#   "arcgislayers",
#   "mapview",
#   "RColorBrewer",
#   "janitor",
#   "here"
# ), dependencies = TRUE)

# Load packages
library(tidyverse)    # Data manipulation and visualization
library(vroom)        # Read rectangular data
library(sf)           # Spatial analysis
library(tidycensus)   # Accessing US Census Data
library(lehdr)        # Access LODES data
library(arcgislayers) # ArcGIS REST API access
library(mapview)      # Interactive mapping
library(RColorBrewer) # Color palettes for maps
library(janitor)      # Data cleaning and preparation
library(here)         # Filepath management

# Install required packages if not available
# pip install numpy pandas geopandas shapely folium requests pygris

# Load processing libraries & modules
import os
from pathlib import Path
import zipfile
import requests
import urllib.parse
import warnings
warnings.filterwarnings('ignore')

# Load data and visualization libraries & modules
import numpy as np
import pandas as pd
import pyproj
import geopandas as gpd
import folium
from shapely.geometry import Point

# Census data query libraries & modules
from pygris import blocks, block_groups
from pygris.helpers import validate_state, validate_county
from pygris.data import get_census, get_lodes

1.2 Project folder

The centralized directory structure organizes input data, processing files, and model outputs. By using the here package in R and pathlib in Python, we ensure that file paths work correctly regardless of the computer or operating system.

R
Python

# Define Paths relative to Project Root
# We want data in: posts/socioeconomic-demo/data
# We want output in: posts/socioeconomic-demo/output

base_dir <- here("posts", "socioeconomic-demo")
data_dir <- file.path(base_dir, "data")
output_dir <- file.path(base_dir, "output")
output_tabular <- file.path(output_dir, "tabular")
output_spatial <- file.path(output_dir, "spatial")

# Create directories if they don't exist
dir.create(data_dir, showWarnings = FALSE, recursive = TRUE)
dir.create(output_tabular, showWarnings = FALSE, recursive = TRUE)
dir.create(output_spatial, showWarnings = FALSE, recursive = TRUE)

# Verify paths
print(paste("Data Directory:", data_dir))

[1] "Data Directory: C:/Users/Pukar.Bhandari/Documents/GitHub/ar-puuk.github.io/posts/socioeconomic-demo/data"

print(paste("Output Directory:", output_dir))

[1] "Output Directory: C:/Users/Pukar.Bhandari/Documents/GitHub/ar-puuk.github.io/posts/socioeconomic-demo/output"

# Define Paths using pathlib
# We retrieve the base_dir from R, but define everything else in Python

# Get base_dir from R environment
base_dir = Path(r.base_dir)

# Define subdirectories relative to base_dir
data_dir = base_dir / "data"
output_dir = base_dir / "output"
output_tabular = output_dir / "tabular"
output_spatial = output_dir / "spatial"

# Create directories if they don't exist
data_dir.mkdir(parents=True, exist_ok=True)
output_tabular.mkdir(parents=True, exist_ok=True)
output_spatial.mkdir(parents=True, exist_ok=True)

# Verify paths
print(f"Data Directory: {data_dir.resolve()}")

Data Directory: C:\Users\Pukar.Bhandari\Documents\GitHub\ar-puuk.github.io\posts\socioeconomic-demo\data

print(f"Output Directory: {output_dir.resolve()}")

Output Directory: C:\Users\Pukar.Bhandari\Documents\GitHub\ar-puuk.github.io\posts\socioeconomic-demo\output

1.3 Set global options and parameters

Configuration settings optimize performance and establish spatial consistency. We set the tigris cache to a local folder within the project to prevent redundant downloads and ensure portability. The South Carolina State Plane coordinate system (EPSG:3361) serves as the standard projection for accurate GIS operations.

R
Python

# Set Tigris Cache to local project folder
tigris_cache <- file.path(data_dir, "tigris")
dir.create(tigris_cache, showWarnings = FALSE)
options(tigris_use_cache = TRUE)
options(tigris_cache_dir = tigris_cache)

# set project CRS
project_crs <- "EPSG:3361"

# Set project CRS
project_crs = "EPSG:3361"

# Note: pygris usually respects the tigris cache environment variable or default locations.
# We can explicitly set cache location for pygris if needed, but here we rely on defaults
# or the R setting if they share environment variables in the session.

1.4 Set census API key

API authentication enables access to detailed demographic and economic datasets from the Census Bureau. The key configuration supports both R and Python environments for automated data retrieval workflows.

Need a Census API key?

Get one for free at census.gov/developers

R
Python

# Set your API key into environment
tidycensus::census_api_key("your_api_key_here", install = TRUE)

# Set your API key into environment
os.environ['CENSUS_API_KEY'] = 'your_api_key_here'

2 Define study area

This section defines the geographic extent of the Lower Savannah Council of Governments region and loads the Traffic Analysis Zone (TAZ) geometry for spatial analysis.

2.1 Define state and counties

The study area encompasses six counties within South Carolina: Aiken, Allendale, Bamberg, Barnwell, Calhoun, and Orangeburg. These counties constitute the LSCOG planning region for travel demand modeling purposes.

R
Python

# Define state abbreviation and county names
state_abb <- "SC"
county_names <- c(
  "Aiken",
  "Allendale",
  "Bamberg",
  "Barnwell",
  "Calhoun",
  "Orangeburg"
)

# Define state abbreviation and county names
state_abb = "SC"
county_names = [
    "Aiken",
    "Allendale",
    "Bamberg",
    "Barnwell",
    "Calhoun",
    "Orangeburg"
]

FIPS code conversion translates state abbreviations and county names into standardized Federal Information Processing Standard codes. These codes enable consistent data retrieval across census datasets and ensure proper geographic matching with demographic and economic data sources.

R
Python

# Converting state abbreviation code to FIPS code
state_fips <- tidycensus:::validate_state(state = state_abb)
county_fips <- vapply(
  county_names,
  function(x) tidycensus:::validate_county(state = state_abb, county = x),
  character(1)
)

# Converting County Names to FIPS code
fips_codes <- paste(state_fips, county_fips, sep = "")
fips_codes

[1] "45003" "45005" "45009" "45011" "45017" "45075"

# Converting state abbreviation code to FIPS code
state_fips = validate_state(state_abb)

Using FIPS code '45' for input 'SC'

# Converting County Names to FIPS code
county_fips = [
    validate_county(state_fips, county)
    for county in county_names
]

Using FIPS code '003' for input 'Aiken'
Using FIPS code '005' for input 'Allendale'
Using FIPS code '009' for input 'Bamberg'
Using FIPS code '011' for input 'Barnwell'
Using FIPS code '017' for input 'Calhoun'
Using FIPS code '075' for input 'Orangeburg'

# Converting County Names to FIPS code
fips_codes = [f"{state_fips}{county}" for county in county_fips]
fips_codes

['45003', '45005', '45009', '45011', '45017', '45075']

2.2 Load TAZ geometry

The TAZ shapefile provides the fundamental spatial framework for travel demand modeling. The geometry is loaded from the TDM exports geodatabase and filtered to include only zones within the six-county study area using FIPS code matching. Coordinate transformation converts the TAZ geometry to the project’s standard coordinate reference system (EPSG:3361) for accurate spatial calculations. The attribute selection retains essential fields including TAZ identifiers, area measurements, area type classifications, and county assignments.

R
Python

# Load TAZ Shapefile
# We assume the input TAZ file is in the data directory
taz_path <- file.path(data_dir, "SE_2019_AD_10_30_2023.gpkg")

lscog_taz <- sf::read_sf(
  taz_path,
  query = paste0(
    "SELECT * FROM \"SE_2019_AD_10_30_2023\" WHERE countyID IN (",
    paste0("'", fips_codes, "'", collapse = ", "),
    ")"
  )
) |>
  sf::st_transform(project_crs) |>
  dplyr::select(
    ID,
    Area,
    Acres,
    TAZ_ID = TAZ_IDs,
    AREA_TYPE,
    COUNTY,
    COUNTYID = countyID
  )

lscog_taz

Simple feature collection with 585 features and 7 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 1691490 ymin: 331894 xmax: 2237471 ymax: 744679.9
Projected CRS: NAD83(HARN) / South Carolina (ft)
# A tibble: 585 × 8
   ID     Area  Acres TAZ_ID AREA_TYPE COUNTY COUNTYID                      geom
   <chr> <dbl>  <dbl> <chr>  <chr>     <chr>  <chr>        <MULTIPOLYGON [foot]>
 1 9050… 13.3   8518. 90501… RURAL     Bambe… 45009    (((2046194 478862.1, 204…
 2 9050… 39.2  25084. 90501… RURAL     Bambe… 45009    (((2012593 500179.5, 201…
 3 7505…  9.03  5778. 75050… RURAL     Orang… 45075    (((2056266 515976, 20561…
 4 7505… 15.5   9934. 75050… RURAL     Orang… 45075    (((2061827 488917.9, 206…
 5 7505… 17.5  11216. 75050… SUBURBAN  Orang… 45075    (((2154222 534929.2, 215…
 6 7505… 12.8   8191. 75050… RURAL     Orang… 45075    (((2195053 518745.9, 219…
 7 7505…  8.10  5181. 75050… SUBURBAN  Orang… 45075    (((2179620 542034.9, 217…
 8 7505… 13.4   8568. 75050… SUBURBAN  Orang… 45075    (((2179131 542424, 21770…
 9 7505…  6.00  3843. 75050… SUBURBAN  Orang… 45075    (((2184440 568297.6, 218…
10 7505…  5.06  3239. 75050… SUBURBAN  Orang… 45075    (((2204464 570787.3, 220…
# ℹ 575 more rows

# Load TAZ Shapefile
taz_path = data_dir / "SE_2019_AD_10_30_2023.gpkg"

lscog_taz = gpd.read_file(
    taz_path,
    layer="SE_2019_AD_10_30_2023",
    where=f"countyID IN ({', '.join([f"'{fips}'" for fips in fips_codes])})"
)

lscog_taz = lscog_taz.to_crs(project_crs)

lscog_taz = lscog_taz.rename(
  columns={
    'TAZ_IDs': 'TAZ_ID',
    'countyID': 'COUNTYID'
})[['ID', 'Area', 'Acres', 'TAZ_ID', 'AREA_TYPE', 'COUNTY', 'COUNTYID', 'geometry']]

lscog_taz

           ID  ...                                           geometry
0     9050130  ...  MULTIPOLYGON (((2046194.08 478862.054, 2046134...
1     9050132  ...  MULTIPOLYGON (((2012593.472 500179.47, 2013190...
2    75050131  ...  MULTIPOLYGON (((2056266.071 515975.986, 205617...
3    75050045  ...  MULTIPOLYGON (((2061826.693 488917.873, 206173...
4    75050182  ...  MULTIPOLYGON (((2154221.857 534929.221, 215441...
..        ...  ...                                                ...
580   5050049  ...  MULTIPOLYGON (((1924248.136 433664.04, 1924004...
581   5050055  ...  MULTIPOLYGON (((1905461.23 427627.626, 1905456...
582   5050066  ...  MULTIPOLYGON (((1880812.362 414251.546, 188090...
583   5050054  ...  MULTIPOLYGON (((1871871.454 413403.458, 187167...
584   5050065  ...  MULTIPOLYGON (((1849074.015 398566.33, 1849218...

[585 rows x 8 columns]

The interactive map visualization displays the TAZ structure colored by county, providing spatial context for the analysis area and enabling quality assurance of the geometric data loading process.

R
Python

# Create interactive map
mapview::mapview(lscog_taz, zcol = "COUNTY", lwd = 1.6, map.types = "CartoDB.Voyager", col.regions = RColorBrewer::brewer.pal(6, "Dark2"))