Case Study 2: ShopSmart Delivery Radius Optimization

Background

ShopSmart, the e-commerce platform from Chapters 3, 20, 23-24, and 26, operates four fulfillment centers (FCs) across the continental United States: New York City, Los Angeles, Chicago, and Dallas. Each FC ships to customers within its assigned region, and the regional assignment is based on the nearest-FC rule: every order ships from whichever FC is closest to the customer.

The problem: customer satisfaction scores have declined for three consecutive quarters, and the analytics team has traced the root cause to delivery times. Orders that arrive in 2 days or fewer have an average satisfaction score of 4.6 out of 5. Orders that take 5+ days average 3.1. The Head of Operations, Priya Sharma, wants to answer three questions:

1. What is the relationship between distance-to-FC and delivery time, and where are the "dead zones" where no FC provides fast service?
2. If ShopSmart opens a fifth fulfillment center, where should it go to maximize the reduction in average delivery time?
3. How much would the optimal fifth FC improve customer satisfaction, and is the investment justified?

The task: Build a geospatial analysis of delivery performance, identify underserved regions, evaluate candidate locations for a fifth FC, and quantify the expected impact on satisfaction.

The Data

import numpy as np
import pandas as pd
import geopandas as gpd
from shapely.geometry import Point
from math import radians, sin, cos, sqrt, atan2
import matplotlib.pyplot as plt

np.random.seed(42)

# --- Fulfillment center data ---
fcs = pd.DataFrame({
    'fc_id': ['FC_NYC', 'FC_LA', 'FC_CHI', 'FC_DAL'],
    'city': ['New York', 'Los Angeles', 'Chicago', 'Dallas'],
    'lat': [40.713, 34.052, 41.878, 32.777],
    'lon': [-74.006, -118.244, -87.630, -96.797],
    'daily_capacity': [9000, 13000, 7000, 6000],
    'avg_processing_hours': [4.2, 3.8, 5.1, 4.5],
})

# --- Haversine distance function ---
def haversine(lon1, lat1, lon2, lat2):
    """Great-circle distance in km."""
    R = 6371
    lon1, lat1, lon2, lat2 = map(radians, [lon1, lat1, lon2, lat2])
    dlat = lat2 - lat1
    dlon = lon2 - lon1
    a = sin(dlat / 2)**2 + cos(lat1) * cos(lat2) * sin(dlon / 2)**2
    return R * 2 * atan2(sqrt(a), sqrt(1 - a))

def haversine_vectorized(lon1, lat1, lon2, lat2):
    """Vectorized haversine distance in km."""
    R = 6371
    lon1, lat1, lon2, lat2 = map(np.radians, [lon1, lat1, lon2, lat2])
    dlat = lat2 - lat1
    dlon = lon2 - lon1
    a = np.sin(dlat / 2)**2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon / 2)**2
    return R * 2 * np.arctan2(np.sqrt(a), np.sqrt(1 - a))

# --- Generate order data ---
# Customer locations clustered around major metro areas
metro_areas = [
    # (lon, lat, spread, relative_weight, name)
    (-74.0, 40.7, 0.8, 0.15, 'NYC metro'),
    (-118.2, 34.0, 1.0, 0.12, 'LA metro'),
    (-87.6, 41.9, 0.6, 0.08, 'Chicago metro'),
    (-96.8, 32.8, 0.7, 0.06, 'Dallas metro'),
    (-122.4, 37.8, 0.5, 0.07, 'SF Bay'),
    (-80.2, 25.8, 0.8, 0.06, 'Miami'),
    (-84.4, 33.7, 0.7, 0.05, 'Atlanta'),
    (-95.4, 29.8, 0.7, 0.05, 'Houston'),
    (-71.1, 42.4, 0.5, 0.04, 'Boston'),
    (-77.0, 38.9, 0.5, 0.04, 'DC'),
    (-105.0, 39.7, 0.6, 0.03, 'Denver'),
    (-112.1, 33.4, 0.6, 0.04, 'Phoenix'),
    (-90.2, 38.6, 0.5, 0.03, 'St Louis'),
    (-86.8, 36.2, 0.5, 0.03, 'Nashville'),
    (-93.3, 45.0, 0.5, 0.03, 'Minneapolis'),
    (-81.7, 41.5, 0.4, 0.03, 'Cleveland'),
    (-83.0, 42.3, 0.4, 0.03, 'Detroit'),
    (-104.9, 38.8, 0.4, 0.02, 'Colorado Springs'),
    (-111.9, 40.8, 0.4, 0.02, 'Salt Lake City'),
    (-122.3, 47.6, 0.4, 0.02, 'Seattle'),
]

n_orders = 10000
orders_list = []

for i in range(n_orders):
    # Pick metro area weighted by population
    weights = [m[3] for m in metro_areas]
    weights = np.array(weights) / sum(weights)
    metro_idx = np.random.choice(len(metro_areas), p=weights)
    metro = metro_areas[metro_idx]

    # Generate customer location near metro center
    clon = metro[0] + np.random.normal(0, metro[2])
    clat = metro[1] + np.random.normal(0, metro[2])

    # Clip to continental US bounds
    clon = np.clip(clon, -125, -66)
    clat = np.clip(clat, 24, 50)

    # Find nearest FC and compute distance
    fc_distances = [
        haversine(clon, clat, row['lon'], row['lat'])
        for _, row in fcs.iterrows()
    ]
    nearest_fc_idx = np.argmin(fc_distances)
    nearest_fc_dist = fc_distances[nearest_fc_idx]
    nearest_fc_id = fcs.iloc[nearest_fc_idx]['fc_id']

    # Delivery time model:
    # base_processing (hours) + transit (distance-dependent) + variability
    processing_hours = fcs.iloc[nearest_fc_idx]['avg_processing_hours']

    # Transit: ~800 km/day for ground shipping
    transit_days = nearest_fc_dist / 800

    # Add variability (weather, carrier delays, weekends)
    variability = np.random.exponential(0.3)

    delivery_days = round((processing_hours / 24) + transit_days + variability, 1)
    delivery_days = max(1.0, delivery_days)

    # Satisfaction model:
    # Base 4.5, penalized by delivery time, with noise
    sat_penalty = max(0, (delivery_days - 2.0)) * 0.35
    satisfaction = round(np.clip(4.5 - sat_penalty + np.random.normal(0, 0.4), 1.0, 5.0), 1)

    # Order value
    order_value = round(np.random.lognormal(4.0, 0.7), 2)

    orders_list.append({
        'order_id': i,
        'customer_lon': round(clon, 4),
        'customer_lat': round(clat, 4),
        'metro_area': metro[4],
        'nearest_fc': nearest_fc_id,
        'dist_to_fc_km': round(nearest_fc_dist, 1),
        'delivery_days': delivery_days,
        'satisfaction': satisfaction,
        'order_value': order_value,
    })

orders = pd.DataFrame(orders_list)

print(f"Orders: {len(orders)}")
print(f"Average delivery time: {orders['delivery_days'].mean():.2f} days")
print(f"Average satisfaction: {orders['satisfaction'].mean():.2f}")
print(f"Orders with 5+ day delivery: {(orders['delivery_days'] >= 5).mean():.1%}")

Step 1: Current State Analysis

# Performance by fulfillment center
fc_performance = (
    orders.groupby('nearest_fc')
    .agg(
        n_orders=('order_id', 'count'),
        pct_orders=('order_id', lambda x: len(x) / len(orders)),
        avg_distance_km=('dist_to_fc_km', 'mean'),
        median_distance_km=('dist_to_fc_km', 'median'),
        avg_delivery_days=('delivery_days', 'mean'),
        avg_satisfaction=('satisfaction', 'mean'),
        pct_5plus_days=('delivery_days', lambda x: (x >= 5).mean()),
        total_revenue=('order_value', 'sum'),
    )
    .round(3)
)

print("Fulfillment Center Performance:")
print(fc_performance)
print(f"\nTotal revenue: ${orders['order_value'].sum():,.0f}")

# Distance bins: the relationship between distance and outcomes
orders['distance_bin'] = pd.cut(
    orders['dist_to_fc_km'],
    bins=[0, 200, 500, 800, 1200, 5000],
    labels=['<200km', '200-500km', '500-800km', '800-1200km', '>1200km']
)

distance_analysis = (
    orders.groupby('distance_bin', observed=True)
    .agg(
        n_orders=('order_id', 'count'),
        avg_delivery_days=('delivery_days', 'mean'),
        avg_satisfaction=('satisfaction', 'mean'),
        pct_5plus_days=('delivery_days', lambda x: (x >= 5).mean()),
    )
    .round(3)
)

print("\nPerformance by Distance to FC:")
print(distance_analysis)

The data confirms the hypothesis: delivery time increases linearly with distance, and satisfaction decreases correspondingly. Orders beyond 800 km from the nearest FC have average delivery times above 3 days and satisfaction below 4.0. Beyond 1,200 km, satisfaction drops below 3.5.

Step 2: Identifying Dead Zones

# Identify the geographic "dead zones" --- areas far from all FCs
# Create a grid of points across the continental US

lon_grid = np.arange(-125, -66, 0.5)
lat_grid = np.arange(24, 50, 0.5)
lon_mesh, lat_mesh = np.meshgrid(lon_grid, lat_grid)
grid_points = pd.DataFrame({
    'lon': lon_mesh.ravel(),
    'lat': lat_mesh.ravel(),
})

# Compute distance to nearest FC for each grid point
for _, fc_row in fcs.iterrows():
    grid_points[f'dist_{fc_row["fc_id"]}'] = haversine_vectorized(
        grid_points['lon'].values, grid_points['lat'].values,
        fc_row['lon'], fc_row['lat']
    )

dist_cols = [f'dist_{fc_id}' for fc_id in fcs['fc_id']]
grid_points['min_dist_km'] = grid_points[dist_cols].min(axis=1)
grid_points['nearest_fc'] = grid_points[dist_cols].idxmin(axis=1)

# Dead zones: grid points where minimum distance > 1000 km
dead_zone = grid_points[grid_points['min_dist_km'] > 1000]
print(f"Grid points total: {len(grid_points)}")
print(f"Dead zone points (>1000 km from any FC): {len(dead_zone)}")
print(f"Dead zone percentage: {len(dead_zone)/len(grid_points):.1%}")

# Where are the dead zones?
print("\nDead zone center of mass:")
print(f"  Lon: {dead_zone['lon'].mean():.1f}")
print(f"  Lat: {dead_zone['lat'].mean():.1f}")
print(f"  (Approximately: {''} the Mountain West / Southeast corridor)")

# Visualize dead zones with matplotlib
fig, ax = plt.subplots(figsize=(14, 8))

scatter = ax.scatter(
    grid_points['lon'], grid_points['lat'],
    c=grid_points['min_dist_km'],
    cmap='RdYlGn_r',  # Red = far, green = close
    s=4, alpha=0.6,
)

# Mark FCs
ax.scatter(
    fcs['lon'], fcs['lat'],
    c='blue', s=200, marker='*', zorder=5, edgecolors='black',
    label='Fulfillment Centers'
)
for _, fc_row in fcs.iterrows():
    ax.annotate(
        fc_row['city'], (fc_row['lon'], fc_row['lat']),
        textcoords='offset points', xytext=(10, 5),
        fontsize=9, fontweight='bold',
    )

plt.colorbar(scatter, label='Distance to Nearest FC (km)', shrink=0.7)
ax.set_xlabel('Longitude')
ax.set_ylabel('Latitude')
ax.set_title('ShopSmart: Distance to Nearest Fulfillment Center')
ax.legend(loc='lower left')
plt.tight_layout()
plt.savefig('shopsmart_dead_zones.png', dpi=150, bbox_inches='tight')
plt.show()

Two dead zones are visible: the Southeast corridor (Florida, Georgia, the Carolinas --- far from NYC, Chicago, and Dallas) and the Mountain West (Montana, Idaho, Utah, Nevada --- far from LA, Chicago, and Dallas).

Step 3: Candidate Location Evaluation

# Evaluate candidate cities for the fifth FC
candidates = pd.DataFrame({
    'city': ['Atlanta', 'Miami', 'Denver', 'Phoenix', 'Nashville',
             'Charlotte', 'Salt Lake City', 'Minneapolis', 'Seattle', 'Houston'],
    'lat': [33.749, 25.762, 39.739, 33.449, 36.163,
            35.227, 40.761, 44.978, 47.607, 29.760],
    'lon': [-84.388, -80.192, -104.990, -112.074, -86.782,
            -80.843, -111.891, -93.265, -122.332, -95.370],
    'estimated_annual_cost_M': [12.5, 14.0, 11.0, 10.5, 10.0,
                                 11.5, 10.0, 11.5, 13.0, 12.0],
})

def evaluate_candidate(candidate_row, orders_df, existing_fcs):
    """
    Evaluate a candidate FC location by computing the new average
    distance and delivery time if this FC were added.
    """
    cand_lon, cand_lat = candidate_row['lon'], candidate_row['lat']

    # Distance from every customer to the candidate
    dist_to_candidate = haversine_vectorized(
        orders_df['customer_lon'].values, orders_df['customer_lat'].values,
        cand_lon, cand_lat
    )

    # New minimum distance: min(current nearest FC, candidate)
    new_min_dist = np.minimum(orders_df['dist_to_fc_km'].values, dist_to_candidate)

    # Orders that would switch to the new FC
    switches = dist_to_candidate < orders_df['dist_to_fc_km'].values

    # New delivery time estimate
    new_delivery = (new_min_dist / 800) + (4.4 / 24) + 0.3  # avg processing + avg variability

    # New satisfaction estimate
    new_sat_penalty = np.maximum(0, (new_delivery - 2.0)) * 0.35
    new_satisfaction = np.clip(4.5 - new_sat_penalty, 1.0, 5.0)

    return {
        'city': candidate_row['city'],
        'orders_switched': switches.sum(),
        'pct_switched': switches.mean(),
        'current_avg_dist': orders_df['dist_to_fc_km'].mean(),
        'new_avg_dist': new_min_dist.mean(),
        'dist_reduction_km': orders_df['dist_to_fc_km'].mean() - new_min_dist.mean(),
        'current_avg_delivery': orders_df['delivery_days'].mean(),
        'new_avg_delivery': new_delivery.mean(),
        'delivery_improvement_days': orders_df['delivery_days'].mean() - new_delivery.mean(),
        'new_avg_satisfaction': new_satisfaction.mean(),
        'sat_improvement': new_satisfaction.mean() - orders_df['satisfaction'].mean(),
        'annual_cost_M': candidate_row['estimated_annual_cost_M'],
    }

results = []
for _, cand in candidates.iterrows():
    results.append(evaluate_candidate(cand, orders, fcs))

results_df = pd.DataFrame(results).sort_values('dist_reduction_km', ascending=False)

print("Candidate FC Evaluation (sorted by distance reduction):")
display_cols = ['city', 'orders_switched', 'pct_switched', 'dist_reduction_km',
                'delivery_improvement_days', 'new_avg_satisfaction', 'sat_improvement',
                'annual_cost_M']
print(results_df[display_cols].round(2).to_string(index=False))

Step 4: ROI Analysis

# Estimate revenue impact of the top candidate
top_candidate = results_df.iloc[0]

print(f"Top candidate: {top_candidate['city']}")
print(f"  Distance reduction: {top_candidate['dist_reduction_km']:.0f} km average")
print(f"  Delivery improvement: {top_candidate['delivery_improvement_days']:.2f} days average")
print(f"  Satisfaction improvement: {top_candidate['sat_improvement']:.3f} points")

# Revenue impact model:
# Higher satisfaction -> lower churn -> retained revenue
# Assumption: 0.1 point satisfaction increase -> 2% churn reduction
# (Based on StreamFlow's satisfaction-churn correlation from Chapter 16)

current_monthly_revenue = orders['order_value'].sum() / 12  # Monthly
churn_reduction_pct = top_candidate['sat_improvement'] * 0.02 / 0.1  # Scaled
annual_retained_revenue = current_monthly_revenue * churn_reduction_pct * 12

print(f"\nROI Estimate:")
print(f"  Current monthly revenue (sample): ${current_monthly_revenue:,.0f}")
print(f"  Estimated churn reduction: {churn_reduction_pct:.1%}")
print(f"  Annual retained revenue (sample): ${annual_retained_revenue:,.0f}")
print(f"  Annual FC cost: ${top_candidate['annual_cost_M']:.1f}M")

# Scale to full customer base (sample is ~10% of actual)
scale_factor = 10
full_annual_retained = annual_retained_revenue * scale_factor
print(f"  Scaled annual retained revenue: ${full_annual_retained:,.0f}")
print(f"  Simple ROI (year 1): {(full_annual_retained - top_candidate['annual_cost_M'] * 1e6) / (top_candidate['annual_cost_M'] * 1e6):.1%}")

Step 5: Visualization for Stakeholders

import folium
from folium.plugins import MarkerCluster

# Create the final map for Priya's presentation
m = folium.Map(location=[39.0, -98.0], zoom_start=4, tiles='cartodbpositron')

# Layer 1: Heatmap of delivery times (using customer points)
# Color code: green (<2 days), yellow (2-4 days), orange (4-6 days), red (6+ days)
def delivery_color(days):
    if days < 2:
        return 'green'
    elif days < 4:
        return 'orange'
    elif days < 6:
        return 'red'
    else:
        return 'darkred'

# Plot a sample of orders (10,000 is too many for browser)
sample_orders = orders.sample(1000, random_state=42)

for _, row in sample_orders.iterrows():
    folium.CircleMarker(
        location=[row['customer_lat'], row['customer_lon']],
        radius=2,
        color=delivery_color(row['delivery_days']),
        fill=True,
        fill_opacity=0.5,
        weight=0.5,
    ).add_to(m)

# Layer 2: Existing FCs (blue stars)
for _, fc_row in fcs.iterrows():
    folium.Marker(
        location=[fc_row['lat'], fc_row['lon']],
        popup=(
            f"<b>{fc_row['city']}</b><br>"
            f"Capacity: {fc_row['daily_capacity']:,}/day<br>"
            f"Processing: {fc_row['avg_processing_hours']}h"
        ),
        tooltip=f"FC: {fc_row['city']}",
        icon=folium.Icon(color='blue', icon='warehouse', prefix='fa'),
    ).add_to(m)

# Layer 3: Top 3 candidate locations (red triangles)
top_3 = results_df.head(3)
for _, cand_row in top_3.iterrows():
    cand_info = candidates[candidates['city'] == cand_row['city']].iloc[0]
    folium.Marker(
        location=[cand_info['lat'], cand_info['lon']],
        popup=(
            f"<b>Candidate: {cand_row['city']}</b><br>"
            f"Orders switched: {cand_row['orders_switched']:,.0f}<br>"
            f"Dist reduction: {cand_row['dist_reduction_km']:.0f} km<br>"
            f"Sat improvement: +{cand_row['sat_improvement']:.3f}<br>"
            f"Annual cost: ${cand_row['annual_cost_M']}M"
        ),
        tooltip=f"Candidate: {cand_row['city']}",
        icon=folium.Icon(color='red', icon='flag', prefix='fa'),
    ).add_to(m)

# Legend (manual HTML)
legend_html = """
<div style="position: fixed; bottom: 50px; left: 50px; z-index: 1000;
     background-color: white; padding: 10px; border: 2px solid grey;
     border-radius: 5px; font-size: 12px;">
<b>Delivery Time</b><br>
<span style="color: green;">&#9679;</span> &lt; 2 days<br>
<span style="color: orange;">&#9679;</span> 2-4 days<br>
<span style="color: red;">&#9679;</span> 4-6 days<br>
<span style="color: darkred;">&#9679;</span> 6+ days<br>
<br>
<span style="color: blue;">&#9733;</span> Existing FC<br>
<span style="color: red;">&#9873;</span> Candidate FC
</div>
"""
m.get_root().html.add_child(folium.Element(legend_html))

m.save('shopsmart_fc_optimization.html')
print("Final map saved to shopsmart_fc_optimization.html")

Findings Summary

Finding 1: Distance drives satisfaction. The relationship is near-linear: every 200 km of additional distance adds roughly 0.25 days of delivery time and reduces satisfaction by 0.09 points. Customers beyond 1,000 km from the nearest FC are the primary source of low satisfaction scores.

Finding 2: Two dead zones exist. The Southeast (centered on Georgia/Florida) and the Mountain West (centered on Utah/Montana) are more than 1,000 km from any FC. These regions account for a disproportionate share of 5+ day deliveries.

Finding 3: The optimal fifth FC location depends on the objective. If maximizing the number of customers served faster, Atlanta or Charlotte wins (large population near the Southeast dead zone). If maximizing average distance reduction across all customers, the result may differ. If minimizing cost per satisfaction-point-gained, Nashville or Denver may be more efficient due to lower operating costs.

Finding 4: The ROI is positive under reasonable assumptions. Even conservative estimates of the satisfaction-to-retention relationship suggest the fifth FC pays for itself within 18-24 months through reduced churn alone, not counting the revenue from faster delivery attracting new customers.

Domain Knowledge --- This is the kind of analysis where geospatial tools provide genuine business value. A non-spatial analysis might identify that "satisfaction is declining" and recommend "improve delivery times." The spatial analysis identifies where the problem is worst, what the root cause is (distance, not processing speed), and exactly where to invest to fix it. That specificity is what turns a data science finding into a capital allocation decision.

Discussion Questions

1. The analysis uses haversine (straight-line) distance as a proxy for shipping distance. Actual shipping routes follow road and rail networks. How would you obtain actual shipping distances, and how much would it change the analysis?

2. The candidate evaluation assumes the new FC would have the same processing time and capacity as existing FCs. What additional data would you need to make the evaluation more realistic?

3. ShopSmart's nearest-FC assignment rule is simple but may not be optimal. If FC_LA is at capacity while FC_DAL has spare capacity, some orders near LA should route to Dallas despite the greater distance. How would you model this capacity-constrained assignment problem?

4. The ROI analysis uses a simple satisfaction-to-churn-reduction multiplier. What are the weaknesses of this approach? How would you build a more rigorous ROI model?

5. A board member asks: "Instead of building a fifth FC, could we just offer free two-day shipping from the existing FCs to high-distance customers?" Outline the analysis you would run to compare these two strategies.

Return to the chapter | Previous: Case Study 1 --- StreamFlow Regional Churn Choropleth