External Station Traffic Forecasts

Single-document version of the external-station forecasting pipeline, originally migrated from 5 sequential notebooks (1-Get-Historic-AADT through 5-Export-Data) as a straight format port. Step 5.5 (CCS-derived truck split) was added afterward and is a genuine methodology change layered on top of that port — see below.

How the final AADT for each external station is assembled, by year:

How the truck/passenger split is assembled, by year — note this is a two-layer process: Steps 5 and 5.5 both run for every row, with 5.5 always having the final say.

Setup

import io
import os
import re
import zipfile
from pathlib import Path

import geopandas as gpd
import numpy as np
import pandas as pd
import pyproj
import requests
from dotenv import load_dotenv
from scipy import stats
from sklearn.linear_model import LinearRegression

load_dotenv()

PROJECT_CRS = "EPSG:3566"

# Year boundaries that recur as control-flow thresholds throughout the pipeline.
# (Distinct from the per-station judgment-call years inside the hardcoded
# reference tables below, e.g. selected_forecast_df — those are literal data,
# not thresholds, so they stay as plain literals.)
YEARS = {
    "HISTORIC_START": 1981,  # first year with any HPMS/UDOT AADT data
    "HISTORIC_END": 2023,  # last year treated as confirmed historic AADT
    "HISTORIC_FILL_CUTOFF": 2024,  # archive fallback is allowed through this year
    "TRUCK_PROJECTION_SPLIT": 2025,  # truck volume switches from %-split to AADT regression
    "PIPELINE_END": 2060,  # last year in the exported forecast grid
    "EXPORT_START": 2010,  # results/external_year_vol.csv is filtered to Year >= this
}
def fetch_github(
    url: str, mode: str = "private", token_env_var: str = "GITHUB_TOKEN"
) -> requests.Response:
    """Fetch content from GitHub repositories."""
    if mode not in ["public", "private"]:
        raise ValueError(f"mode must be 'public' or 'private', got '{mode}'")

    if mode == "public":
        response = requests.get(url, timeout=30)
    else:
        token = os.getenv(token_env_var)
        if not token:
            raise ValueError(
                f"GitHub token not found in environment variable '{token_env_var}'. "
                f"Check your .env file has: {token_env_var}=your_token_here"
            )

        headers = {
            "Authorization": f"token {token}",
            "Accept": "application/vnd.github.v3.raw",
        }
        response = requests.get(url, headers=headers, timeout=30)

    response.raise_for_status()
    return response


def format_external_label(df: pd.DataFrame) -> pd.Series:
    """Builds the 'Ext # <id> - <name>' label used across plots/exports."""
    return (
        "Ext # "
        + df["externalid"].fillna(0).astype(int).astype(str)
        + " - "
        + df["name"].fillna("")
    )

Step 1: Historic AADT (was 1-Get-Historic-AADT)

externals_df = pd.read_csv("params/externals.csv")
aadt_df = pd.read_csv(
    "data/AADTHistory 2023.xlsx - UnroundedAADT2023.csv", low_memory=False
)
external_segment_link = pd.read_csv("params/externals-segments-link.csv")

# Manual Override: make it closer to the CCS station #-319
external_segment_link.loc[external_segment_link["externalid"] == 3620, "segid"] = (
    "0189_013.1"  # Alternative: 0189_009.9 # Original: 0189_014.3
)
# This is the single source of truth for the link + override — reused as-is in Step 5.
# Melt AADT columns
aadt_year_cols = [col for col in aadt_df.columns if re.match(r"AADT\d{4}$", col)]
aadt_melted = aadt_df.melt(
    id_vars=["STATION", "RouteID", "BeginPoint", "EndPoint", "Section_Length", "DESC"],
    value_vars=aadt_year_cols,
    var_name="year",
    value_name="AADT",
)
aadt_melted["year"] = aadt_melted["year"].str.extract(r"(\d{4})").astype(int)

# Extract route number
aadt_melted["route"] = (
    aadt_melted["RouteID"].str.extract(r"^(\d{4})")[0].dropna().astype(int)
)
aadt_melted = aadt_melted[["route", "BeginPoint", "EndPoint", "year", "AADT"]]
aadt_melted.columns = ["route", "mp_begin", "mp_end", "year", "AADT"]
aadt_melted
route mp_begin mp_end year AADT
0 6.0 0.000 46.0380 2023 457.0
1 6.0 46.038 77.5560 2023 409.0
2 6.0 77.556 82.8970 2023 586.0
3 6.0 82.897 83.9110 2023 2189.0
4 6.0 83.911 87.6940 2023 4012.0
... ... ... ... ... ...
196376 3468.0 0.000 5.1248 1981 0.0
196377 3469.0 0.000 6.9300 1981 0.0
196378 3470.0 0.000 1.0394 1981 0.0
196379 3478.0 0.000 2.0404 1981 0.0
196380 3483.0 0.000 1.7030 1981 0.0

196381 rows × 5 columns

# get external route, mp
external_segment_link[["route", "milepost"]] = external_segment_link[
    "segid"
].str.split("_", expand=True)
external_segment_link["route"] = external_segment_link["route"].astype(int)
external_segment_link["milepost"] = external_segment_link["milepost"].astype(float)
external_segment_link
externalid segid route milepost
0 3601 1082_000.0 1082 0.00
1 3602 0013_006.5 13 6.50
2 3603 1112_000.0 1112 0.00
3 3604 0015_368.1 15 368.10
4 3605 0038_003.2 38 3.20
5 3606 0091_010.1 91 10.10
6 3607 3462_002.8 3462 2.80
7 3608 0039_008.7 39 8.70
8 3609 0084_087.8 84 87.80
9 3610 2688_005.5 2688 5.50
10 3611 0201_000.0 201 0.00
11 3612 0080_101.2 80 101.20
12 3613 0065_008.4 65 8.40
13 3614 0080_138.9 80 138.90
14 3615 0190_015.75 190 15.75
15 3616 1828_004.5 1828 4.50
16 3617 0006_140.3 6 140.30
17 3618 0073_015.5 73 15.50
18 3619 3108_000.0 3108 0.00
19 3620 0189_013.1 189 13.10
20 3621 2865_019.4 2865 19.40
21 3622 2863_000.0 2863 0.00
22 3623 0006_216.2 6 216.20
23 3624 0096_006.2 96 6.20
24 3625 2495_000.0 2495 0.00
25 3626 0089_296.0 89 296.00
26 3627 1822_000.0 1822 0.00
27 3628 0015_233.2 15 233.20
28 3629 1826_004.9 1826 4.90 
# Define full list of years
years_full_df = pd.DataFrame(
    {"year": list(range(YEARS["HISTORIC_START"], YEARS["HISTORIC_END"] + 1))}
)

# Create full cartesian product: segments × years
full_index_df = external_segment_link[["externalid", "segid"]].merge(
    years_full_df, how="cross"
)
full_index_df
externalid segid year
0 3601 1082_000.0 1981
1 3601 1082_000.0 1982
2 3601 1082_000.0 1983
3 3601 1082_000.0 1984
4 3601 1082_000.0 1985
... ... ... ...
1242 3629 1826_004.9 2019
1243 3629 1826_004.9 2020
1244 3629 1826_004.9 2021
1245 3629 1826_004.9 2022
1246 3629 1826_004.9 2023

1247 rows × 3 columns

# Merge external_segment_link and aadt_melted with left join
matched_df = external_segment_link.merge(aadt_melted, on=["route"], how="left").query(
    "mp_begin < milepost + .01 <= mp_end"
)
matched_df
externalid segid route milepost mp_begin mp_end year AADT
0 3601 1082_000.0 1082 0.0 0.0 12.348 2023.0 817.0
4 3601 1082_000.0 1082 0.0 0.0 12.348 2022.0 783.0
8 3601 1082_000.0 1082 0.0 0.0 12.348 2021.0 792.0
12 3601 1082_000.0 1082 0.0 0.0 12.348 2020.0 739.0
16 3601 1082_000.0 1082 0.0 0.0 12.348 2019.0 740.0
... ... ... ... ... ... ... ... ...
33529 3629 1826_004.9 1826 4.9 0.0 4.944 1985.0 0.0
33532 3629 1826_004.9 1826 4.9 0.0 4.944 1984.0 0.0
33535 3629 1826_004.9 1826 4.9 0.0 4.944 1983.0 0.0
33538 3629 1826_004.9 1826 4.9 0.0 4.944 1982.0 0.0
33541 3629 1826_004.9 1826 4.9 0.0 4.944 1981.0 0.0

989 rows × 8 columns

full_set_df = pd.merge(
    full_index_df, matched_df, how="left", on=["externalid", "segid", "year"]
)[["externalid", "segid", "year", "AADT"]].assign(
    AADT=lambda d: d["AADT"].replace(0, np.nan)
)
full_set_df
externalid segid year AADT
0 3601 1082_000.0 1981 NaN
1 3601 1082_000.0 1982 NaN
2 3601 1082_000.0 1983 NaN
3 3601 1082_000.0 1984 NaN
4 3601 1082_000.0 1985 NaN
... ... ... ... ...
1242 3629 1826_004.9 2019 3158.0
1243 3629 1826_004.9 2020 3155.0
1244 3629 1826_004.9 2021 3382.0
1245 3629 1826_004.9 2022 3345.0
1246 3629 1826_004.9 2023 3489.0

1247 rows × 4 columns

full_set_df.to_csv("intermediate/external-historic-aadt.csv", index=False)

Step 2: Previous Forecasts v900 (was 2-Prepare-Previous-Forecasts)

previous_forecasts_df = pd.read_csv("data/external-forecasts-v900.csv")
previous_forecasts_df = previous_forecasts_df.melt(
    id_vars=["TAZ"], var_name="year", value_name="AADT"
)
previous_forecasts_df["forecast"] = "v900"
previous_forecasts_df.rename(columns={"TAZ": "externalid"}, inplace=True)
previous_forecasts_df
externalid year AADT forecast
0 3601 2020 750.0 v900
1 3602 2020 9570.0 v900
2 3603 2020 1480.0 v900
3 3604 2020 28390.0 v900
4 3605 2020 3840.0 v900
... ... ... ... ...
140 3625 2060 200.0 v900
141 3626 2060 5250.0 v900
142 3627 2060 450.0 v900
143 3628 2060 54260.0 v900
144 3629 2060 3660.0 v900

145 rows × 4 columns

previous_forecasts_df.to_csv("intermediate/previous-forecasts.csv", index=False)

Step 3: Linear Forecasts (was 3-Prepare-Linear-Forecasts)

# list of years for which to produce linear forecasts
future_years = [2027, 2032, 2036, 2046, 2055, 2060]
# One projection group per "Since <year>" starting point, all ending at
# HISTORIC_END with no excluded years -- a pure derived sequence, unlike
# selected_forecast_df/manual_adj_df below, which are genuine per-station
# judgment calls and stay as literal tables.
df_proj_groups_linear = pd.DataFrame(
    [
        [f"Since {year}", year, YEARS["HISTORIC_END"], {}]
        for year in range(YEARS["HISTORIC_START"], YEARS["HISTORIC_END"])
    ],
    columns=("pg_name", "pg_year_from", "pg_year_to", "pg_years_exclude"),
)
df_proj_groups_linear
pg_name pg_year_from pg_year_to pg_years_exclude
0 Since 1981 1981 2023 {}
1 Since 1982 1982 2023 {}
2 Since 1983 1983 2023 {}
3 Since 1984 1984 2023 {}
4 Since 1985 1985 2023 {}
5 Since 1986 1986 2023 {}
6 Since 1987 1987 2023 {}
7 Since 1988 1988 2023 {}
8 Since 1989 1989 2023 {}
9 Since 1990 1990 2023 {}
10 Since 1991 1991 2023 {}
11 Since 1992 1992 2023 {}
12 Since 1993 1993 2023 {}
13 Since 1994 1994 2023 {}
14 Since 1995 1995 2023 {}
15 Since 1996 1996 2023 {}
16 Since 1997 1997 2023 {}
17 Since 1998 1998 2023 {}
18 Since 1999 1999 2023 {}
19 Since 2000 2000 2023 {}
20 Since 2001 2001 2023 {}
21 Since 2002 2002 2023 {}
22 Since 2003 2003 2023 {}
23 Since 2004 2004 2023 {}
24 Since 2005 2005 2023 {}
25 Since 2006 2006 2023 {}
26 Since 2007 2007 2023 {}
27 Since 2008 2008 2023 {}
28 Since 2009 2009 2023 {}
29 Since 2010 2010 2023 {}
30 Since 2011 2011 2023 {}
31 Since 2012 2012 2023 {}
32 Since 2013 2013 2023 {}
33 Since 2014 2014 2023 {}
34 Since 2015 2015 2023 {}
35 Since 2016 2016 2023 {}
36 Since 2017 2017 2023 {}
37 Since 2018 2018 2023 {}
38 Since 2019 2019 2023 {}
39 Since 2020 2020 2023 {}
40 Since 2021 2021 2023 {}
41 Since 2022 2022 2023 {} 
# import historic AADT (created in step 1)
df_historic_aadt = pd.read_csv("intermediate/external-historic-aadt.csv")
df_historic_aadt
externalid segid year AADT
0 3601 1082_000.0 1981 NaN
1 3601 1082_000.0 1982 NaN
2 3601 1082_000.0 1983 NaN
3 3601 1082_000.0 1984 NaN
4 3601 1082_000.0 1985 NaN
... ... ... ... ...
1242 3629 1826_004.9 2019 3158.0
1243 3629 1826_004.9 2020 3155.0
1244 3629 1826_004.9 2021 3382.0
1245 3629 1826_004.9 2022 3345.0
1246 3629 1826_004.9 2023 3489.0

1247 rows × 4 columns

Linear forecasts with assist from ChatGPT: https://chat.openai.com/share/d127492a-ad78-4f45-afd0-50e29069db1a

def fit_linear_forecast(
    filtered_group: pd.DataFrame,
    pg_name: str,
    pg_year_from: int,
    future_years: list[int],
) -> tuple[np.ndarray, str] | None:
    """
    Fits an OLS trend (AADT ~ year) on `filtered_group` and predicts
    `pg_year_from` plus each of `future_years`.

    Returns (forecast_array, proj_grp_used) if there are at least 2 valid
    (non-NaN) historic points to fit a line through; otherwise None, leaving
    the "no data" fallback to the caller.
    """
    X = filtered_group["year"].values.reshape(-1, 1)
    y = filtered_group["AADT"].values

    valid_mask = ~np.isnan(X.flatten()) & ~np.isnan(y)
    X_valid = X[valid_mask]
    y_valid = y[valid_mask]

    if len(X_valid) < 2:
        return None

    model = LinearRegression()
    model.fit(X_valid, y_valid)
    aadt = model.predict(np.array([pg_year_from] + future_years).reshape(-1, 1))
    aadt = np.rint(aadt).astype(int)
    return aadt, pg_name
# Initialize a list to store the individual result DataFrames
forecast_results_list = []

# Initialize a set to track which externalids have already been assigned 'No Data'
no_data_externalids = set()

# Open the error file
with open("intermediate/linear-forecasts-errors.txt", "w") as err_file:
    # Loop through the projection groups
    for _, row in df_proj_groups_linear.iterrows():
        pg_name = row["pg_name"]
        pg_year_from = row["pg_year_from"]
        pg_year_to = row["pg_year_to"]
        pg_years_exclude = set(row["pg_years_exclude"])

        # Group by externalid and segid and iterate through the groups
        for (externalid, segid), group in df_historic_aadt.groupby(
            ["externalid", "segid"]
        ):
            # Filter the data according to the projection group criteria
            filtered_group = group[
                (group["year"] >= pg_year_from) & (group["year"] <= pg_year_to)
            ]
            filtered_group = filtered_group[
                ~filtered_group["year"].isin(pg_years_exclude)
            ]

            fit_result = fit_linear_forecast(
                filtered_group, pg_name, pg_year_from, future_years
            )

            if fit_result is not None:
                aadt, proj_grp_used = fit_result
            else:
                # "No Data" is recorded once per externalid: the first
                # projection group (in "Since <year>" order, oldest first)
                # that lacks 2+ historic points wins; every later group for
                # the same externalid is skipped entirely via this memo set,
                # since a shorter window can't fix a lack of data.
                if externalid in no_data_externalids:
                    continue  # Skip this segid
                error_msg = f"No valid data for externalid: {externalid}, segid: {segid}, Projection Group: {pg_name}. Filling zeros."
                err_file.write(error_msg + "\n")
                aadt = np.zeros(len([pg_year_from] + future_years), dtype=int)
                proj_grp_used = "No Data"
                no_data_externalids.add(externalid)

            # Create a dictionary to store results for this group
            result_dict = {
                "externalid": externalid,
                "segid": segid,
                "PROJ_GRP": proj_grp_used,
            }
            result_dict.update(
                {
                    year: forecast
                    for year, forecast in zip([pg_year_from] + future_years, aadt)
                }
            )

            # Convert the dictionary to a DataFrame and add to the list
            result_df = pd.DataFrame([result_dict])
            result_df_melt = result_df.melt(
                id_vars=["externalid", "segid", "PROJ_GRP"],
                var_name="year",
                value_name="linear_forecast",
            )
            forecast_results_list.append(result_df_melt)

# Concatenate all the individual result DataFrames
forecast_results = pd.concat(forecast_results_list, ignore_index=True)
forecast_results
externalid segid PROJ_GRP year linear_forecast
0 3601 1082_000.0 Since 1981 1981 -26
1 3601 1082_000.0 Since 1981 2027 890
2 3601 1082_000.0 Since 1981 2032 989
3 3601 1082_000.0 Since 1981 2036 1069
4 3601 1082_000.0 Since 1981 2046 1268
... ... ... ... ... ...
6799 3629 1826_004.9 Since 2022 2032 4785
6800 3629 1826_004.9 Since 2022 2036 5361
6801 3629 1826_004.9 Since 2022 2046 6801
6802 3629 1826_004.9 Since 2022 2055 8097
6803 3629 1826_004.9 Since 2022 2060 8817

6804 rows × 5 columns

forecast_results.to_csv("results/linear-forecasts.csv", index=False)

Step 4: Finalize Forecasts (was 4-Finalize-Forecasts)

linear_forecast_df = pd.read_csv("results/linear-forecasts.csv")
externals_df = pd.read_csv("params/externals.csv")
previous_forecasts_df = pd.read_csv("intermediate/previous-forecasts.csv")
external_xy_df = pd.read_csv("data/external-x-y.csv")

# Create a lookup dictionary
name_lookup = externals_df.set_index("externalid")["name"].to_dict()

# Map names to the forecast df
linear_forecast_df["name"] = linear_forecast_df["externalid"].map(name_lookup)
linear_forecast_df["external"] = format_external_label(linear_forecast_df)

# Map names to the previous-forecasts df (used later for plotting/reference only)
previous_forecasts_df["name"] = previous_forecasts_df["externalid"].map(name_lookup)
previous_forecasts_df["external"] = format_external_label(previous_forecasts_df)
previous_forecasts_df = previous_forecasts_df[
    previous_forecasts_df["year"] > YEARS["HISTORIC_END"]
]
selected_forecast_df = pd.DataFrame(
    [
        [3601, "Since 2003", "as early as there is data"],
        [
            3602,
            "Since 1981",
            "as early as there is data, and trying not to get too high when using 2002 or higher",
        ],
        [
            3603,
            "Since 2016",
            "trying to split difference between forecast using earlier data and forecast using 2020+ data",
        ],
        [
            3604,
            "Since 1981",
            "as early as there is data, trend looks pretty consistent over time",
        ],
        [3605, "Since 2002", "earlier data jumps areound too much"],
        [
            3606,
            "Since 1981",
            "as early as there is data, trend looks pretty consistent over time",
        ],
        [3607, "Since 2007", "data does not look consistent till 2007"],
        [
            3608,
            "Since 1981",
            "as early as there is data, trend looks pretty consistent over time",
        ],
        [
            3609,
            "Since 1981",
            "as early as there is data, trend looks pretty consistent over time",
        ],
        [3610, "No Data", ""],
        [
            3611,
            "Since 1981",
            "as early as there is data, trend looks pretty consistent over time",
        ],
        [
            3612,
            "Since 1981",
            "as early as there is data, trend looks pretty consistent over time",
        ],
        [3613, "Since 2009", "pre 2009 data is much higher than post 2009 data"],
        [
            3614,
            "Since 1981",
            "as early as there is data, trend looks pretty consistent over time",
        ],
        [
            3615,
            "Since 2007",
            "erratic data. 2007 tries to split difference between 2004-2012 data and 2013-2023 data",
        ],
        [3616, "Since 2006", "data pretty consistent since 2006"],
        [3617, "Since 2011", "data pretty consistent since 2011"],
        [
            3618,
            "Since 2007",
            "data pretty consistent since 2007, but forecast from 1981 is much higher than 2007-2023 data",
        ],
        [3619, "No Data", ""],
        [
            3620,
            "Since 1981",
            "as early as there is data, trend looks pretty consistent over time",
        ],
        [3621, "No Data", ""],
        [3622, "No Data", ""],
        [
            3623,
            "Since 1981",
            "as early as there is data, trend looks pretty consistent over time. outlier in 1983 does not change the trend",
        ],
        [
            3624,
            "Since 1990",
            "data pretty consistent since 2017, but puts forecast too high",
        ],
        [3625, "No Data", ""],
        [
            3626,
            "Since 1981",
            "as early as there is data, despite some variation early on, trend looks pretty consistent over time",
        ],
        [3627, "No Data", ""],
        [
            3628,
            "Since 1984",
            "as early as there is data, trend looks pretty consistent over time except for past two years, keep eye on it",
        ],
        [3629, "Since 2011", "data pretty consistent since 2011"],
    ],
    columns=["externalid", "PROJ_GRP", "linear_forecast_notes"],
)

linear_forecast_selected_df = pd.merge(
    selected_forecast_df, linear_forecast_df, on=["externalid", "PROJ_GRP"]
)
linear_forecast_selected_df
externalid PROJ_GRP linear_forecast_notes segid year linear_forecast name external
0 3601 Since 2003 as early as there is data 1082_000.0 2003 412 FAR-1082 Bird Refuge Ext # 3601 - FAR-1082 Bird Refuge
1 3601 Since 2003 as early as there is data 1082_000.0 2027 890 FAR-1082 Bird Refuge Ext # 3601 - FAR-1082 Bird Refuge
2 3601 Since 2003 as early as there is data 1082_000.0 2032 989 FAR-1082 Bird Refuge Ext # 3601 - FAR-1082 Bird Refuge
3 3601 Since 2003 as early as there is data 1082_000.0 2036 1069 FAR-1082 Bird Refuge Ext # 3601 - FAR-1082 Bird Refuge
4 3601 Since 2003 as early as there is data 1082_000.0 2046 1268 FAR-1082 Bird Refuge Ext # 3601 - FAR-1082 Bird Refuge
... ... ... ... ... ... ... ... ...
198 3629 Since 2011 data pretty consistent since 2011 1826_004.9 2032 4248 FAR-1826 South Ridge Farms Ext # 3629 - FAR-1826 South Ridge Farms
199 3629 Since 2011 data pretty consistent since 2011 1826_004.9 2036 4588 FAR-1826 South Ridge Farms Ext # 3629 - FAR-1826 South Ridge Farms
200 3629 Since 2011 data pretty consistent since 2011 1826_004.9 2046 5437 FAR-1826 South Ridge Farms Ext # 3629 - FAR-1826 South Ridge Farms
201 3629 Since 2011 data pretty consistent since 2011 1826_004.9 2055 6202 FAR-1826 South Ridge Farms Ext # 3629 - FAR-1826 South Ridge Farms
202 3629 Since 2011 data pretty consistent since 2011 1826_004.9 2060 6627 FAR-1826 South Ridge Farms Ext # 3629 - FAR-1826 South Ridge Farms

203 rows × 8 columns

def custom_round(x: float) -> float:
    """Rounds to the nearest 10/50/100/500/1000, coarser at larger magnitudes."""
    if x < 100:
        return round(x / 10) * 10
    elif x < 1000:
        return round(x / 50) * 50
    elif x < 10000:
        return round(x / 100) * 100
    elif x < 100000:
        return round(x / 500) * 500
    else:
        return round(x / 1000) * 1000
# Example manual adjustment data
manual_adj_df = pd.DataFrame(
    [
        [3602, 2027, 500],
        [3602, 2032, 300],
        [3602, 2036, 200],
        [3602, 2046, 100],
        [3605, 2027, 100],
        [3608, 2027, 100],
        [3612, 2027, 4000],
        [3612, 2032, 2000],
        [3612, 2036, 1000],
        [3613, 2027, 50],
        [3613, 2032, 0],
        [3613, 2036, 0],
        [3613, 2046, 0],
        [3613, 2055, 0],
        [3615, 2027, -2050],  # adjust down to prior forecast
        [3615, 2032, -2200],  # adjust down to prior forecast
        [3615, 2036, -2300],  # adjust down to prior forecast
        [3615, 2046, -2600],  # adjust down to prior forecast
        [3615, 2055, -2750],  # adjust down to prior forecast
        [3615, 2060, -2900],  # adjust down to prior forecast
        [3620, 2027, 500],
        [3623, 2027, 500],
        [3623, 2032, 500],
        [3623, 2036, 500],
        [3626, 2027, 300],
        [3626, 2032, 200],
        [3626, 2036, 100],
        [3610, 2027, 100],
        [3610, 2032, 100],
        [3610, 2036, 100],
        [3610, 2046, 150],
        [3610, 2055, 150],
        [3610, 2060, 200],
        [3619, 2027, 450],
        [3619, 2032, 450],
        [3619, 2036, 500],
        [3619, 2046, 550],
        [3619, 2055, 550],
        [3619, 2060, 600],
        [3621, 2027, 350],
        [3621, 2032, 350],
        [3621, 2036, 400],
        [3621, 2046, 400],
        [3621, 2055, 450],
        [3621, 2060, 450],
        [3622, 2027, 150],
        [3622, 2032, 150],
        [3622, 2036, 150],
        [3622, 2046, 200],
        [3622, 2055, 200],
        [3622, 2060, 200],
        [3624, 2027, 600],
        [3624, 2032, 400],
        [3624, 2036, 400],
        [3624, 2046, 300],
        [3624, 2055, 100],
        [3624, 2060, 0],
        [3625, 2027, 150],
        [3625, 2032, 150],
        [3625, 2036, 150],
        [3625, 2046, 200],
        [3625, 2055, 200],
        [3625, 2060, 200],
        [3627, 2027, 350],
        [3627, 2032, 350],
        [3627, 2036, 400],
        [3627, 2046, 400],
        [3627, 2055, 450],
        [3627, 2060, 450],
    ],
    columns=["externalid", "year", "manual_adj"],
)

# Filter only needed years
forecast_df = linear_forecast_selected_df[
    linear_forecast_selected_df["year"] >= YEARS["HISTORIC_END"]
]

# Merge manual adjustments
final_forecast_df = pd.merge(
    forecast_df, manual_adj_df, on=["externalid", "year"], how="left"
)

# Compute final_forecast using manual adjustments where available
final_forecast_df["final_forecast"] = final_forecast_df[
    "linear_forecast"
] + final_forecast_df["manual_adj"].fillna(0)

# Apply custom rounding
final_forecast_df["final_forecast"] = final_forecast_df["final_forecast"].apply(
    custom_round
)
final_forecast_df
externalid PROJ_GRP linear_forecast_notes segid year linear_forecast name external manual_adj final_forecast
0 3601 Since 2003 as early as there is data 1082_000.0 2027 890 FAR-1082 Bird Refuge Ext # 3601 - FAR-1082 Bird Refuge NaN 900
1 3601 Since 2003 as early as there is data 1082_000.0 2032 989 FAR-1082 Bird Refuge Ext # 3601 - FAR-1082 Bird Refuge NaN 1000
2 3601 Since 2003 as early as there is data 1082_000.0 2036 1069 FAR-1082 Bird Refuge Ext # 3601 - FAR-1082 Bird Refuge NaN 1100
3 3601 Since 2003 as early as there is data 1082_000.0 2046 1268 FAR-1082 Bird Refuge Ext # 3601 - FAR-1082 Bird Refuge NaN 1300
4 3601 Since 2003 as early as there is data 1082_000.0 2055 1447 FAR-1082 Bird Refuge Ext # 3601 - FAR-1082 Bird Refuge NaN 1400
... ... ... ... ... ... ... ... ... ... ...
169 3629 Since 2011 data pretty consistent since 2011 1826_004.9 2032 4248 FAR-1826 South Ridge Farms Ext # 3629 - FAR-1826 South Ridge Farms NaN 4200
170 3629 Since 2011 data pretty consistent since 2011 1826_004.9 2036 4588 FAR-1826 South Ridge Farms Ext # 3629 - FAR-1826 South Ridge Farms NaN 4600
171 3629 Since 2011 data pretty consistent since 2011 1826_004.9 2046 5437 FAR-1826 South Ridge Farms Ext # 3629 - FAR-1826 South Ridge Farms NaN 5400
172 3629 Since 2011 data pretty consistent since 2011 1826_004.9 2055 6202 FAR-1826 South Ridge Farms Ext # 3629 - FAR-1826 South Ridge Farms NaN 6200
173 3629 Since 2011 data pretty consistent since 2011 1826_004.9 2060 6627 FAR-1826 South Ridge Farms Ext # 3629 - FAR-1826 South Ridge Farms NaN 6600

174 rows × 10 columns

# Export data for the externals-app Shiny app.

# Define transformer from UTM Zone 12N (EPSG:32612) to WGS84 Lat/Lon (EPSG:4326)
transformer = pyproj.Transformer.from_crs("EPSG:32612", "EPSG:4326", always_xy=True)

# Apply transformation to create 'lon' and 'lat' columns
external_xy_df["lon"], external_xy_df["lat"] = transformer.transform(
    external_xy_df["x_utm12n"].values,
    external_xy_df["y_utm12n"].values,
)

aadt_df = pd.read_csv("intermediate/external-historic-aadt.csv")
aadt_df["name"] = aadt_df["externalid"].map(name_lookup)
aadt_df["external"] = format_external_label(aadt_df)

# historic.csv: actual historic AADT by station/year (used as the "observed" series).
aadt_df.drop(columns=["external", "segid", "name"]).dropna().astype(int).to_csv(
    "externals-app/data/historic.csv", index=False
)
# forecasts-previous.csv: the prior (v9.0) forecast, kept for before/after comparison.
previous_forecasts_df.drop(columns=["external", "forecast", "name"]).rename(
    columns={"AADT": "previous_forecast"}
).sort_values(by=["externalid", "year"]).astype(int).to_csv(
    "externals-app/data/forecasts-previous.csv", index=False
)
# externals.csv: station locations (lat/lon) for the app's map view.
external_xy_df.to_csv("externals-app/data/externals.csv", index=False)
# forecasts.csv: the current final forecast (linear trend + manual adjustments + rounding).
final_forecast_df.drop(
    columns=["external", "name", "linear_forecast_notes", "segid", "PROJ_GRP"]
).fillna(0).astype(int).to_csv("externals-app/data/forecasts.csv", index=False)
# linear-forecasts.csv: the selected linear-trend line (one projection group per
# station, via selected_forecast_df above) that the app draws as the forecast trend.
linear_forecast_selected_df[linear_forecast_selected_df["linear_forecast"] > 0].drop(
    columns=["PROJ_GRP", "linear_forecast_notes", "segid", "name", "external"]
).fillna(0).astype(int).to_csv("externals-app/data/linear-forecasts.csv", index=False)
final_forecast_df.to_csv("results/final_forecast_df.csv", index=False)

Step 5: Export Data (was 5-Export-Data)

Helper functions

def sophisticated_fill(group: pd.DataFrame) -> pd.Series:
    """Interpolates AADT_Forecast gaps and back-extrapolates from the earliest trend."""
    # A. PREPARE DATA
    s = group.set_index("Year")["AADT_Forecast"].astype(float)

    # B. INTERPOLATION (Fill gaps between known years)
    s_interp = s.interpolate(method="linear", limit_area="inside")

    # C. BACKWARD EXTRAPOLATION (Fill 2023-2026)
    valid_years = s.dropna().index.sort_values()

    if len(valid_years) >= 2:
        y1, y2 = valid_years[0], valid_years[1]
        val1, val2 = s.loc[y1], s.loc[y2]

        if y2 != y1:
            slope = (val2 - val1) / (y2 - y1)

            # Find years before the first known data point
            pre_years = s_interp.index[s_interp.index < y1]
            for yr in pre_years:
                if pd.isna(s_interp.loc[yr]):
                    extrapolated_val = val1 + slope * (yr - y1)
                    s_interp.loc[yr] = max(0, extrapolated_val)

    # Fill any remaining edges
    s_final = s_interp.bfill().ffill()

    # D. RE-ALIGNMENT
    return pd.Series(s_final.loc[group["Year"]].values, index=group.index)
def project_future_volumes(
    group: pd.DataFrame,
    value_col: str,
    year_col: str,
    independent_col: str,
    split_year: int = YEARS["TRUCK_PROJECTION_SPLIT"],
) -> pd.Series:
    """
    Projects future volumes using Linear Regression (y = mx + b).

    Logic:
    1. If data varies: Calculate proper Slope (m) and Intercept (b) using OLS regression.
    2. If AADT is constant (Edge Case): Fallback to a Ratio Slope (b=0) to allow
       growth, as a unique m & b cannot be calculated from a single point.
    """
    # 1. Identify Historic Data for Training
    hist_mask = (group[year_col] < split_year) & (group[value_col].notna())

    # 2. If insufficient historic data, forward fill the last known value
    if hist_mask.sum() < 2:
        last_known = (
            group.loc[group[year_col] < split_year, value_col].iloc[-1]
            if hist_mask.any()
            else 0
        )
        group.loc[group[year_col] >= split_year, value_col] = last_known
        return group[value_col]

    # 3. Prepare Data
    x_hist = group.loc[hist_mask, independent_col].values
    y_hist = group.loc[hist_mask, value_col].values

    # 4. Check for Constant AADT (Edge Case where Regression Fails)
    if np.min(x_hist) == np.max(x_hist):
        # Fallback: Assume Intercept is 0 and calculate Slope based on Average Ratio.
        mean_x = np.mean(x_hist)
        if mean_x == 0:
            slope = 0
        else:
            slope = np.mean(y_hist) / mean_x
        intercept = 0
    else:
        # 5. Standard Calculation (The "Proper" Way)
        slope, intercept, _, _, _ = stats.linregress(x_hist, y_hist)

    # 6. Predict Future
    future_mask = group[year_col] >= split_year
    x_future = group.loc[future_mask, independent_col].values

    predicted = slope * x_future + intercept
    group.loc[future_mask, value_col] = np.maximum(predicted, 0)

    return group[value_col]

Input data

forecast_results = pd.read_csv("results/final_forecast_df.csv")
forecast_results[["externalid", "year", "final_forecast"]]
externalid year final_forecast
0 3601 2027 900
1 3601 2032 1000
2 3601 2036 1100
3 3601 2046 1300
4 3601 2055 1400
... ... ... ...
169 3629 2032 4200
170 3629 2036 4600
171 3629 2046 5400
172 3629 2055 6200
173 3629 2060 6600

174 rows × 3 columns

gdf_master_segments = gpd.read_file(
    "zip://data/updated-traffic-factors/Master_Segs_withFactors_20251120.zip"
).to_crs(PROJECT_CRS)
gdf_master_segments
SEGID BMP EMP DISTANCE CO_FIPS PLANAREA AADT2023 AADT2022 AADT2021 AADT2020 ... FAC_WIN FAC_SPR FAC_SUM FAC_FAL FAC_MAXMO FAC_MAX FACMANADJ SUTRUCKS CUTRUCKS geometry
0 0006_000.0 0.000 0.665 0.666641 27 UDOT 457.0 441.0 474.0 430.0 ... 0.8769 1.0071 1.0496 1.0664 10 1.1275 0 0.2496 0.2324 LINESTRING (916692.404 6835614.744, 920182.189...
1 0006_000.7 0.665 16.022 15.369839 27 UDOT 457.0 441.0 474.0 430.0 ... 0.8769 1.0071 1.0496 1.0664 10 1.1275 0 0.2496 0.2324 LINESTRING (920182.189 6835168.248, 923166.425...
2 0006_016.0 16.022 46.017 30.001961 27 UDOT 457.0 441.0 474.0 430.0 ... 0.8769 1.0071 1.0496 1.0664 10 1.1275 0 0.2496 0.2324 LINESTRING (999430.463 6834408.584, 999447.155...
3 0006_046.0 46.017 60.218 14.194306 27 UDOT 409.0 395.0 424.0 385.0 ... 0.8769 1.0071 1.0496 1.0664 10 1.1275 0 0.1751 0.3338 LINESTRING (1143701.911 6830874.803, 1145024.9...
4 0006_060.2 60.218 77.545 17.323237 27 UDOT 409.0 395.0 424.0 385.0 ... 0.8769 1.0071 1.0496 1.0664 10 1.1275 0 0.1751 0.3338 LINESTRING (1206604.946 6871212.763, 1206701.0...
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
9252 WFRC_8489 0.000 0.000 0.505514 35 WFRC 0.0 0.0 0.0 0.0 ... 0.9411 0.9924 1.0246 1.0419 8 1.0510 0 0.1086 0.0565 LINESTRING (1518588.41 7388436.516, 1518656.54...
9253 WFRC_8490 0.000 0.000 0.737896 35 WFRC 0.0 0.0 0.0 0.0 ... 0.9411 0.9924 1.0246 1.0419 8 1.0510 0 0.1086 0.0565 LINESTRING (1521238.95 7388420.585, 1522418.15...
9254 WFRC_8491 0.000 0.000 0.265495 35 WFRC 0.0 0.0 0.0 0.0 ... 0.9411 0.9924 1.0246 1.0419 8 1.0510 0 0.1086 0.0565 LINESTRING (1525103.396 7388831.353, 1525935.1...
9255 WFRC_8492 0.000 0.000 0.444725 35 WFRC 0.0 0.0 0.0 0.0 ... 0.9519 1.0137 1.0141 1.0203 5 1.0436 0 0.1052 0.0432 LINESTRING (1527315.01 7447662.776, 1527326.78...
9256 WFRC_8493 0.000 0.000 0.583552 35 WFRC 0.0 0.0 0.0 0.0 ... 0.9060 1.0249 1.0414 1.0278 5 1.0822 0 0.1244 0.0925 LINESTRING (1507812.247 7347416.236, 1507019.7...

9257 rows × 82 columns

# Read Processed UDOT AADT Data directly from GitHub repo
response = fetch_github(
    "https://raw.githubusercontent.com/WFRCAnalytics/DATA-UDOT-AADT-Processing/refs/heads/main/_output/udot_aadt_trkpct_data.csv",
    mode="private",
)
df_aadt_udot = pd.read_csv(io.StringIO(response.text))
df_aadt_udot
Station RouteID BeginPoint EndPoint SectionLength DESC_ YEAR AADT SUTRK CUTRK
0 001-0010 0015PM 109.029 112.071 3.042 SR 160 South Beaver Milford 1981 4125.0 0.029041 0.222085
1 001-0010 0015PM 109.029 112.071 3.042 SR 160 South Beaver Milford 1982 4325.0 0.029041 0.222085
2 001-0010 0015PM 109.029 112.071 3.042 SR 160 South Beaver Milford 1983 4635.0 0.029041 0.222085
3 001-0010 0015PM 109.029 112.071 3.042 SR 160 South Beaver Milford 1984 4810.0 0.029041 0.222085
4 001-0010 0015PM 109.029 112.071 3.042 SR 160 South Beaver Milford 1985 5120.0 0.029041 0.222085
... ... ... ... ... ... ... ... ... ... ...
199535 057-1530 3424PM 0.553 1.306 0.753 9th St (Rt 3426) via Polk Ave - Sheridan Dr 2020 2219.0 0.000000 0.000000
199536 057-1530 3424PM 0.553 1.306 0.753 9th St (Rt 3426) via Polk Ave - Sheridan Dr 2021 2405.0 0.000000 0.000000
199537 057-1530 3424PM 0.553 1.306 0.753 9th St (Rt 3426) via Polk Ave - Sheridan Dr 2022 2429.0 0.000000 0.000000
199538 057-1530 3424PM 0.553 1.306 0.753 9th St (Rt 3426) via Polk Ave - Sheridan Dr 2023 2463.0 0.000000 0.000000
199539 057-1530 3424PM 0.553 1.306 0.753 9th St (Rt 3426) via Polk Ave - Sheridan Dr 2024 2517.0 0.000000 0.000000

199540 rows × 10 columns

Reusing external_segment_link from Step 1 rather than re-reading it — same file, same CCS-319 override, so there’s a single place that logic can be edited.

external_segment_link
externalid segid route milepost
0 3601 1082_000.0 1082 0.00
1 3602 0013_006.5 13 6.50
2 3603 1112_000.0 1112 0.00
3 3604 0015_368.1 15 368.10
4 3605 0038_003.2 38 3.20
5 3606 0091_010.1 91 10.10
6 3607 3462_002.8 3462 2.80
7 3608 0039_008.7 39 8.70
8 3609 0084_087.8 84 87.80
9 3610 2688_005.5 2688 5.50
10 3611 0201_000.0 201 0.00
11 3612 0080_101.2 80 101.20
12 3613 0065_008.4 65 8.40
13 3614 0080_138.9 80 138.90
14 3615 0190_015.75 190 15.75
15 3616 1828_004.5 1828 4.50
16 3617 0006_140.3 6 140.30
17 3618 0073_015.5 73 15.50
18 3619 3108_000.0 3108 0.00
19 3620 0189_013.1 189 13.10
20 3621 2865_019.4 2865 19.40
21 3622 2863_000.0 2863 0.00
22 3623 0006_216.2 6 216.20
23 3624 0096_006.2 96 6.20
24 3625 2495_000.0 2495 0.00
25 3626 0089_296.0 89 296.00
26 3627 1822_000.0 1822 0.00
27 3628 0015_233.2 15 233.20
28 3629 1826_004.9 1826 4.90 
# Prior published version of external_year_vol.csv — used below as a fallback source
# for AWDT_FAC and historic AADT/truck-% where the primary sources have gaps.
df_external_year_archive = pd.read_csv(r"archive/v920/external_year_vol.csv")
df_external_year_archive
;Idx_WF WF_Ext Year AWDT PASS_VOL TRUCK_MD TRUCK_HV AWDT_FAC AADT PASSENGER TRUCK_SU TRUCK_MU PctTrk_SU PctTrk_MU
0 36012010 3601 2010 492 321 105 66 0.956 515 336 110 69 0.214 0.135
1 36012011 3601 2011 487 317 104 66 0.956 510 332 109 69 0.214 0.135
2 36012012 3601 2012 558 364 119 75 0.956 585 381 125 79 0.214 0.135
3 36012013 3601 2013 577 377 123 77 0.956 605 395 129 81 0.214 0.135
4 36012014 3601 2014 587 382 126 79 0.956 615 400 132 83 0.214 0.135
... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
1474 36292056 3629 2056 3409 2291 637 481 0.984 3464 2328 647 489 0.187 0.141
1475 36292057 3629 2057 3457 2323 646 488 0.984 3513 2361 656 496 0.187 0.141
1476 36292058 3629 2058 3505 2356 654 495 0.984 3562 2394 665 503 0.187 0.141
1477 36292059 3629 2059 3553 2388 663 502 0.984 3611 2427 674 510 0.187 0.141
1478 36292060 3629 2060 3602 2421 672 509 0.984 3660 2460 683 517 0.187 0.141

1479 rows × 14 columns

Prepare data

Step 1: Initialize final dataframe structure

df_external_year = (
    pd.MultiIndex.from_product(
        [
            external_segment_link["externalid"].unique(),
            range(YEARS["HISTORIC_START"], YEARS["PIPELINE_END"] + 1),
        ],
        names=["WF_Ext", "Year"],
    )
    .to_frame(index=False)
    .reset_index(drop=True)
)

# Create Index String
df_external_year[";Idx_WF"] = df_external_year["WF_Ext"].astype(
    str
) + df_external_year["Year"].astype(str)

# Map Metadata
df_external_year["segid"] = df_external_year["WF_Ext"].map(
    external_segment_link.set_index("externalid")["segid"]
)
df_external_year["route"] = df_external_year["segid"].str.split("_").str[0] + "PM"
df_external_year["milepost"] = pd.to_numeric(
    df_external_year["segid"].str.split("_").str[1], errors="coerce"
)

# Map Station Name
df_external_year["Ext_Name"] = df_external_year["WF_Ext"].map(
    forecast_results[["externalid", "external"]]
    .drop_duplicates()
    .set_index("externalid")["external"]
)
df_external_year
WF_Ext Year ;Idx_WF segid route milepost Ext_Name
0 3601 1981 36011981 1082_000.0 1082PM 0.0 Ext # 3601 - FAR-1082 Bird Refuge
1 3601 1982 36011982 1082_000.0 1082PM 0.0 Ext # 3601 - FAR-1082 Bird Refuge
2 3601 1983 36011983 1082_000.0 1082PM 0.0 Ext # 3601 - FAR-1082 Bird Refuge
3 3601 1984 36011984 1082_000.0 1082PM 0.0 Ext # 3601 - FAR-1082 Bird Refuge
4 3601 1985 36011985 1082_000.0 1082PM 0.0 Ext # 3601 - FAR-1082 Bird Refuge
... ... ... ... ... ... ... ...
2315 3629 2056 36292056 1826_004.9 1826PM 4.9 Ext # 3629 - FAR-1826 South Ridge Farms
2316 3629 2057 36292057 1826_004.9 1826PM 4.9 Ext # 3629 - FAR-1826 South Ridge Farms
2317 3629 2058 36292058 1826_004.9 1826PM 4.9 Ext # 3629 - FAR-1826 South Ridge Farms
2318 3629 2059 36292059 1826_004.9 1826PM 4.9 Ext # 3629 - FAR-1826 South Ridge Farms
2319 3629 2060 36292060 1826_004.9 1826PM 4.9 Ext # 3629 - FAR-1826 South Ridge Farms

2320 rows × 7 columns

Step 2: AWDT factors (with fallback)

# Primary: From Master Segments
df_external_year["AWDT_FAC"] = df_external_year["segid"].map(
    gdf_master_segments.set_index("SEGID")["FAC_WDAVG"]
)

# Fallback: From Archive (Match WF_Ext + Year)
archive_fac_map = df_external_year_archive.set_index(["WF_Ext", "Year"])["AWDT_FAC"]
df_external_year = df_external_year.set_index(["WF_Ext", "Year"])
df_external_year["AWDT_FAC"] = df_external_year["AWDT_FAC"].fillna(archive_fac_map)
df_external_year = df_external_year.reset_index()

# Clean up: Replace 0s with 1 (or NA) to prevent math errors, assuming 1 if missing
df_external_year["AWDT_FAC"] = df_external_year["AWDT_FAC"].replace(0, 1).fillna(1)
df_external_year
WF_Ext Year ;Idx_WF segid route milepost Ext_Name AWDT_FAC
0 3601 1981 36011981 1082_000.0 1082PM 0.0 Ext # 3601 - FAR-1082 Bird Refuge 1.0000
1 3601 1982 36011982 1082_000.0 1082PM 0.0 Ext # 3601 - FAR-1082 Bird Refuge 1.0000
2 3601 1983 36011983 1082_000.0 1082PM 0.0 Ext # 3601 - FAR-1082 Bird Refuge 1.0000
3 3601 1984 36011984 1082_000.0 1082PM 0.0 Ext # 3601 - FAR-1082 Bird Refuge 1.0000
4 3601 1985 36011985 1082_000.0 1082PM 0.0 Ext # 3601 - FAR-1082 Bird Refuge 1.0000
... ... ... ... ... ... ... ... ...
2315 3629 2056 36292056 1826_004.9 1826PM 4.9 Ext # 3629 - FAR-1826 South Ridge Farms 1.0404
2316 3629 2057 36292057 1826_004.9 1826PM 4.9 Ext # 3629 - FAR-1826 South Ridge Farms 1.0404
2317 3629 2058 36292058 1826_004.9 1826PM 4.9 Ext # 3629 - FAR-1826 South Ridge Farms 1.0404
2318 3629 2059 36292059 1826_004.9 1826PM 4.9 Ext # 3629 - FAR-1826 South Ridge Farms 1.0404
2319 3629 2060 36292060 1826_004.9 1826PM 4.9 Ext # 3629 - FAR-1826 South Ridge Farms 1.0404

2320 rows × 8 columns

Step 3: Historic data integration (1981 - 2022)

# 1. Merge UDOT Data based on Year and Route
udot_subset = df_aadt_udot[
    ["YEAR", "RouteID", "BeginPoint", "EndPoint", "AADT", "SUTRK", "CUTRK"]
].copy()

df_merged = df_external_year.merge(
    udot_subset, left_on=["Year", "route"], right_on=["YEAR", "RouteID"], how="left"
)

# 2. Filter for segment matching
# Primary Check: Does the milepost fall strictly INSIDE the segment?
df_merged["is_inside"] = (df_merged["milepost"] >= df_merged["BeginPoint"]) & (
    df_merged["milepost"] < df_merged["EndPoint"]
)

# Fallback Check: Distance to the BeginPoint (in case of edge cases/missing endpoints)
df_merged["distance"] = (df_merged["BeginPoint"] - df_merged["milepost"]).abs()

# Sort so that rows where 'is_inside' is True float to the top.
df_sorted = df_merged.sort_values(by=["is_inside", "distance"], ascending=[False, True])

# 3. Aggregate: Take the best match
df_historic_matches = df_sorted.groupby([";Idx_WF"]).first().reset_index()

# Update the main df with these matches
cols_to_update = {"AADT": "AADT_Historic", "SUTRK": "PctTrk_SU", "CUTRK": "PctTrk_MU"}
for src, dest in cols_to_update.items():
    df_external_year[dest] = df_external_year[";Idx_WF"].map(
        df_historic_matches.set_index(";Idx_WF")[src]
    )

# 4. Fallback: Fill Historic Gaps from Archive
mask_historic_fill = df_external_year["Year"] <= YEARS["HISTORIC_FILL_CUTOFF"]

fallback_cols = ["AADT", "PctTrk_SU", "PctTrk_MU"]
target_cols = ["AADT_Historic", "PctTrk_SU", "PctTrk_MU"]

for archive_col, target_col in zip(fallback_cols, target_cols):
    archive_map = df_external_year_archive.set_index(["WF_Ext", "Year"])[archive_col]

    fallback_series = pd.Series(
        df_external_year.set_index(["WF_Ext", "Year"]).index.map(archive_map),
        index=df_external_year.index,
    )

    df_external_year.loc[mask_historic_fill, target_col] = df_external_year.loc[
        mask_historic_fill, target_col
    ].fillna(fallback_series)

df_external_year
WF_Ext Year ;Idx_WF segid route milepost Ext_Name AWDT_FAC AADT_Historic PctTrk_SU PctTrk_MU
0 3601 1981 36011981 1082_000.0 1082PM 0.0 Ext # 3601 - FAR-1082 Bird Refuge 1.0000 0.0 0.0 0.0
1 3601 1982 36011982 1082_000.0 1082PM 0.0 Ext # 3601 - FAR-1082 Bird Refuge 1.0000 0.0 0.0 0.0
2 3601 1983 36011983 1082_000.0 1082PM 0.0 Ext # 3601 - FAR-1082 Bird Refuge 1.0000 0.0 0.0 0.0
3 3601 1984 36011984 1082_000.0 1082PM 0.0 Ext # 3601 - FAR-1082 Bird Refuge 1.0000 0.0 0.0 0.0
4 3601 1985 36011985 1082_000.0 1082PM 0.0 Ext # 3601 - FAR-1082 Bird Refuge 1.0000 0.0 0.0 0.0
... ... ... ... ... ... ... ... ... ... ... ...
2315 3629 2056 36292056 1826_004.9 1826PM 4.9 Ext # 3629 - FAR-1826 South Ridge Farms 1.0404 NaN NaN NaN
2316 3629 2057 36292057 1826_004.9 1826PM 4.9 Ext # 3629 - FAR-1826 South Ridge Farms 1.0404 NaN NaN NaN
2317 3629 2058 36292058 1826_004.9 1826PM 4.9 Ext # 3629 - FAR-1826 South Ridge Farms 1.0404 NaN NaN NaN
2318 3629 2059 36292059 1826_004.9 1826PM 4.9 Ext # 3629 - FAR-1826 South Ridge Farms 1.0404 NaN NaN NaN
2319 3629 2060 36292060 1826_004.9 1826PM 4.9 Ext # 3629 - FAR-1826 South Ridge Farms 1.0404 NaN NaN NaN

2320 rows × 11 columns

Step 4: Forecast data & AADT hybridization

# 1. Map Forecast Results (Sparse data: 2027, 2032, ... 2060 etc.)
df_external_year["AADT_Forecast"] = (
    df_external_year.set_index(["WF_Ext", "Year"])
    .index.map(forecast_results.set_index(["externalid", "year"])["final_forecast"])
    .astype(float)
)

# Extrapolate/Back-cast Forecast (Fill 2023-2026 using trend from 2027+)
forecast_mask = df_external_year["Year"] >= YEARS["HISTORIC_END"]

df_external_year.loc[forecast_mask, "AADT_Forecast"] = (
    df_external_year[forecast_mask]
    .copy()
    .groupby("WF_Ext", group_keys=False)[["Year", "AADT_Forecast"]]
    .apply(sophisticated_fill)
)

# Create Combined AADT Column (The "Final" Column)
# Logic: Use Historic if < YEARS["HISTORIC_END"], else use Forecast
df_external_year["AADT"] = np.where(
    df_external_year["Year"] < YEARS["HISTORIC_END"],
    df_external_year["AADT_Historic"],
    df_external_year["AADT_Forecast"],
)

df_external_year["AADT"] = (
    pd.to_numeric(df_external_year["AADT"], errors="coerce").round().astype("Int64")
)
df_external_year
WF_Ext Year ;Idx_WF segid route milepost Ext_Name AWDT_FAC AADT_Historic PctTrk_SU PctTrk_MU AADT_Forecast AADT
0 3601 1981 36011981 1082_000.0 1082PM 0.0 Ext # 3601 - FAR-1082 Bird Refuge 1.0000 0.0 0.0 0.0 NaN 0
1 3601 1982 36011982 1082_000.0 1082PM 0.0 Ext # 3601 - FAR-1082 Bird Refuge 1.0000 0.0 0.0 0.0 NaN 0
2 3601 1983 36011983 1082_000.0 1082PM 0.0 Ext # 3601 - FAR-1082 Bird Refuge 1.0000 0.0 0.0 0.0 NaN 0
3 3601 1984 36011984 1082_000.0 1082PM 0.0 Ext # 3601 - FAR-1082 Bird Refuge 1.0000 0.0 0.0 0.0 NaN 0
4 3601 1985 36011985 1082_000.0 1082PM 0.0 Ext # 3601 - FAR-1082 Bird Refuge 1.0000 0.0 0.0 0.0 NaN 0
... ... ... ... ... ... ... ... ... ... ... ... ... ...
2315 3629 2056 36292056 1826_004.9 1826PM 4.9 Ext # 3629 - FAR-1826 South Ridge Farms 1.0404 NaN NaN NaN 6280.0 6280
2316 3629 2057 36292057 1826_004.9 1826PM 4.9 Ext # 3629 - FAR-1826 South Ridge Farms 1.0404 NaN NaN NaN 6360.0 6360
2317 3629 2058 36292058 1826_004.9 1826PM 4.9 Ext # 3629 - FAR-1826 South Ridge Farms 1.0404 NaN NaN NaN 6440.0 6440
2318 3629 2059 36292059 1826_004.9 1826PM 4.9 Ext # 3629 - FAR-1826 South Ridge Farms 1.0404 NaN NaN NaN 6520.0 6520
2319 3629 2060 36292060 1826_004.9 1826PM 4.9 Ext # 3629 - FAR-1826 South Ridge Farms 1.0404 NaN NaN NaN 6600.0 6600

2320 rows × 13 columns

Step 5: Truck volume projections

# A. Calculate Base Historic Volumes (Rows < 2025)
df_external_year["TRUCK_SU"] = (
    df_external_year["AADT"] * df_external_year["PctTrk_SU"].fillna(0)
).where(df_external_year["PctTrk_SU"].notna())

df_external_year["TRUCK_MU"] = (
    df_external_year["AADT"] * df_external_year["PctTrk_MU"].fillna(0)
).where(df_external_year["PctTrk_MU"].notna())

# B. Project Future Volumes (Rows >= 2025)
# split_year=2025 lines up with mask_historic_fill (<= 2024) above: historic
# truck percentages are trusted through 2024, so 2025+ switches to regressing
# truck volume against AADT instead. TRUCK_SU and TRUCK_MU are independent
# regressions (both only depend on Year/AADT, not on each other), so both are
# computed in a single per-WF_Ext groupby pass rather than two.
def project_future_truck_volumes(
    group: pd.DataFrame, split_year: int = YEARS["TRUCK_PROJECTION_SPLIT"]
) -> pd.DataFrame:
    """Projects TRUCK_SU and TRUCK_MU forward from `split_year` for one WF_Ext group."""
    group = group.copy()
    group["TRUCK_SU"] = project_future_volumes(
        group, "TRUCK_SU", "Year", independent_col="AADT", split_year=split_year
    )
    group["TRUCK_MU"] = project_future_volumes(
        group, "TRUCK_MU", "Year", independent_col="AADT", split_year=split_year
    )
    return group[["TRUCK_SU", "TRUCK_MU"]]


df_external_year[["TRUCK_SU", "TRUCK_MU"]] = df_external_year.groupby(
    "WF_Ext", group_keys=False
).apply(project_future_truck_volumes, include_groups=False)

# Round Truck Volumes to Integers
df_external_year["TRUCK_SU"] = (
    pd.to_numeric(df_external_year["TRUCK_SU"]).round().astype("Int64")
)
df_external_year["TRUCK_MU"] = (
    pd.to_numeric(df_external_year["TRUCK_MU"]).round().astype("Int64")
)

# C. Back-Calculate Future Percentages
mask_future = df_external_year["Year"] >= YEARS["TRUCK_PROJECTION_SPLIT"]

df_external_year.loc[mask_future, "PctTrk_SU"] = (
    df_external_year.loc[mask_future, "TRUCK_SU"]
    / df_external_year.loc[mask_future, "AADT"]
)
df_external_year.loc[mask_future, "PctTrk_MU"] = (
    df_external_year.loc[mask_future, "TRUCK_MU"]
    / df_external_year.loc[mask_future, "AADT"]
)

# D. Global Rounding (Applies to ALL rows: Historic & Forecast)
df_external_year["PctTrk_SU"] = df_external_year["PctTrk_SU"].astype(float).round(4)
df_external_year["PctTrk_MU"] = df_external_year["PctTrk_MU"].astype(float).round(4)
df_external_year
WF_Ext Year ;Idx_WF segid route milepost Ext_Name AWDT_FAC AADT_Historic PctTrk_SU PctTrk_MU AADT_Forecast AADT TRUCK_SU TRUCK_MU
0 3601 1981 36011981 1082_000.0 1082PM 0.0 Ext # 3601 - FAR-1082 Bird Refuge 1.0000 0.0 0.0 0.0 NaN 0 0 0
1 3601 1982 36011982 1082_000.0 1082PM 0.0 Ext # 3601 - FAR-1082 Bird Refuge 1.0000 0.0 0.0 0.0 NaN 0 0 0
2 3601 1983 36011983 1082_000.0 1082PM 0.0 Ext # 3601 - FAR-1082 Bird Refuge 1.0000 0.0 0.0 0.0 NaN 0 0 0
3 3601 1984 36011984 1082_000.0 1082PM 0.0 Ext # 3601 - FAR-1082 Bird Refuge 1.0000 0.0 0.0 0.0 NaN 0 0 0
4 3601 1985 36011985 1082_000.0 1082PM 0.0 Ext # 3601 - FAR-1082 Bird Refuge 1.0000 0.0 0.0 0.0 NaN 0 0 0
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
2315 3629 2056 36292056 1826_004.9 1826PM 4.9 Ext # 3629 - FAR-1826 South Ridge Farms 1.0404 NaN 0.0 0.0 6280.0 6280 0 0
2316 3629 2057 36292057 1826_004.9 1826PM 4.9 Ext # 3629 - FAR-1826 South Ridge Farms 1.0404 NaN 0.0 0.0 6360.0 6360 0 0
2317 3629 2058 36292058 1826_004.9 1826PM 4.9 Ext # 3629 - FAR-1826 South Ridge Farms 1.0404 NaN 0.0 0.0 6440.0 6440 0 0
2318 3629 2059 36292059 1826_004.9 1826PM 4.9 Ext # 3629 - FAR-1826 South Ridge Farms 1.0404 NaN 0.0 0.0 6520.0 6520 0 0
2319 3629 2060 36292060 1826_004.9 1826PM 4.9 Ext # 3629 - FAR-1826 South Ridge Farms 1.0404 NaN 0.0 0.0 6600.0 6600 0 0

2320 rows × 15 columns

Step 5.5: CCS-derived / regionally-adjusted truck split (all years)

Steps A-D above give each station’s “existing” truck-split percentage for every year (historic UDOT SUTRK/CUTRK where available, regression-projected for Year >= YEARS["TRUCK_PROJECTION_SPLIT"]). This step now overrides or rescales that series, for every year 1981-2060, station by station:

  • Matched stations (have a CCS site with usable Tue/Wed/Thu length-bin count data): replace PctTrk_SU/PctTrk_MU with one fixed percentage derived from that station’s own vehicle-length distribution, held constant across all years.
  • Unmatched stations: keep the existing (already year-varying) percentage series, but rescale it by a constant regional adjustment ratio, keyed by the station’s road’s Functional Class. Stations whose Functional Class is Local or Unpaved (no ratio defined for either) fall back to the statewide “All Stations” ratio.

Matched-station CCS split

# Combine both CCS extracts' length-bin count data (WFRC_260202.zip is the original
# WFRC-county extract; 20260728_UDOT_CCS is a small supplemental pull covering a few
# sites outside that coverage, e.g. 307 and 615).
with zipfile.ZipFile("data/WFRC_260202.zip") as z:
    with z.open("JKLP_Length_WFRC_2023_Sept_Nov.csv") as f:
        jklp_a = pd.read_csv(f)

# The supplemental extract is a raw DB export with a trailing "N rows selected." footer.
jklp_b = pd.read_csv(
    "data/20260728_UDOT_CCS/JKLP_Length_WFRC_2023_Sept_Nov.csv",
    skipfooter=2,
    engine="python",
)
jklp_b["SITE"] = jklp_b["SITE"].astype(int)

jklp_all = pd.concat([jklp_a, jklp_b], ignore_index=True)
jklp_all["abs_site"] = jklp_all["SITE"].abs()
jklp_all["weekday"] = pd.to_datetime(jklp_all["DATE_ONLY"]).dt.day_name()
jklp_tue_wed_thu = jklp_all[jklp_all["weekday"].isin(["Tuesday", "Wednesday", "Thursday"])]

# Length bin -> LT(=Passenger)/MD/HV %, pre-collapsed (no FHWA-13-class intermediate step).
split_factors = pd.read_csv("data/LT_MD_HV_Split_Factors_2025.csv").dropna(
    subset=["Length Bin"]
).set_index("Length Bin")

ccs_connections = pd.read_csv("data/External Stations - CCS Connections.csv")


def compute_ccs_split(
    site_candidates: list[int], jklp_df: pd.DataFrame, split_factors_df: pd.DataFrame
) -> tuple[int, float, float, float] | None:
    """
    Tries each candidate CCS SITE id in order (a station may list more than
    one candidate site); returns (site_used, LT%, MD%, HV%) from the first
    one with usable Tue/Wed/Thu length-bin count data, or None if none of
    the candidates have any.
    """
    for site in site_candidates:
        by_bin = jklp_df.loc[jklp_df["abs_site"] == site].groupby("LENGTH_DEF")["COUNT"].sum()
        if by_bin.sum() <= 0:
            continue
        lt = md = hv = 0.0
        for bin_name, count in by_bin.items():
            if bin_name not in split_factors_df.index:
                continue
            lt += count * split_factors_df.loc[bin_name, "LT"] / 100
            md += count * split_factors_df.loc[bin_name, "MD"] / 100
            hv += count * split_factors_df.loc[bin_name, "HV"] / 100
        total = lt + md + hv
        if total > 0:
            return site, lt / total, md / total, hv / total
    return None


matched_pct = {}
for _, row in ccs_connections.iterrows():
    site_str = row["SITE"]
    if pd.isna(site_str) or str(site_str).strip() == "":
        continue
    candidates = [int(float(s.strip())) for s in str(site_str).split(",")]
    result = compute_ccs_split(candidates, jklp_tue_wed_thu, split_factors)
    if result is not None:
        _, lt_pct, md_pct, hv_pct = result
        matched_pct[int(row["N"])] = {"MD": md_pct, "HV": hv_pct}

matched_pct
{3606: {'MD': np.float64(0.06361746125939362),
  'HV': np.float64(0.07374354638507385)},
 3608: {'MD': np.float64(0.07110812283272898),
  'HV': np.float64(0.02000357349906918)},
 3609: {'MD': np.float64(0.05903544234668708),
  'HV': np.float64(0.22682773235967182)},
 3612: {'MD': np.float64(0.045269547236972886),
  'HV': np.float64(0.10275848270535058)},
 3614: {'MD': np.float64(0.05967990691219469),
  'HV': np.float64(0.10012126083565437)},
 3617: {'MD': np.float64(0.06668134306281431),
  'HV': np.float64(0.0761972208526444)},
 3620: {'MD': np.float64(0.14604709705960156),
  'HV': np.float64(0.06068312833315179)},
 3626: {'MD': np.float64(0.07838720924864032),
  'HV': np.float64(0.04144919157502055)},
 3628: {'MD': np.float64(0.060398702531700106),
  'HV': np.float64(0.17646955601772002)}}

Functional Class lookup (for unmatched stations)

def ft_to_functional_class(ft: int) -> str | None:
    """Maps a WFv1000 master-net FT_2023 facility-type code to the ratio
    table's Functional Class categories, per the UDOT/WFRC FT code legend."""
    if ft == 2:
        return "Principal Arterial"
    if ft == 3:
        return "Minor Arterial"
    if ft in (4, 5):
        return "Collector"
    if ft == 6:
        return "Local"
    if ft == 7:
        return "Unpaved"
    if 12 <= ft <= 15:
        return "Expressway"
    if 20 <= ft <= 42:
        return "Freeway"
    return None


master_net_nodes = gpd.read_file(
    "data/WFv1000_MasterNet_20260430/WFv1000_MasterNet_20260430_Node.shp"
)
master_net_links = gpd.read_file(
    "data/WFv1000_MasterNet_20260430/WFv1000_MasterNet_20260430_Link.shp"
)

external_node_ids = master_net_nodes.loc[master_net_nodes["EXTERNAL"] == 1, "N"].tolist()
connector_links = master_net_links[master_net_links["EXTERNAL"] == 1]

# Each external node's own link is a dummy centroid connector (FT_2023 == 1 always) --
# the real Functional Class comes from the link one hop further, at the far-side node.
functional_class_by_external = {}
for n in external_node_ids:
    connectors = connector_links[(connector_links["A"] == n) | (connector_links["B"] == n)]
    far_nodes = set(connectors["A"]).union(connectors["B"]) - {n}
    ft_values = set()
    for far in far_nodes:
        real_links = master_net_links[
            ((master_net_links["A"] == far) | (master_net_links["B"] == far))
            & ~((master_net_links["A"] == n) | (master_net_links["B"] == n))
        ]
        ft_values.update(real_links["FT_2023"].dropna().unique().tolist())
    classes = {ft_to_functional_class(ft) for ft in ft_values} - {None}
    # A station with 0 or >1 distinct classes falls back to "All Stations" below.
    functional_class_by_external[int(n)] = classes.pop() if len(classes) == 1 else None

functional_class_by_external
{3601: 'Collector',
 3602: 'Minor Arterial',
 3603: 'Collector',
 3604: 'Freeway',
 3605: 'Minor Arterial',
 3606: 'Principal Arterial',
 3607: 'Collector',
 3608: 'Minor Arterial',
 3609: 'Freeway',
 3610: 'Local',
 3611: 'Expressway',
 3612: 'Freeway',
 3613: 'Unpaved',
 3614: 'Freeway',
 3615: 'Minor Arterial',
 3616: 'Collector',
 3617: 'Principal Arterial',
 3618: 'Principal Arterial',
 3619: 'Collector',
 3620: 'Expressway',
 3621: 'Unpaved',
 3622: 'Unpaved',
 3623: 'Expressway',
 3624: 'Expressway',
 3625: 'Unpaved',
 3626: 'Expressway',
 3627: 'Unpaved',
 3628: 'Freeway',
 3629: 'Collector'}

Regional adjustment ratio lookup

reclass_ratio = pd.read_csv("data/CCS_Reclass_Ratio_2025_2003.csv").set_index(
    "Functional Class"
)
all_stations_ratio = reclass_ratio.loc["All Stations"]


def ratio_for_class(functional_class: str | None) -> pd.Series:
    """Looks up the MD/HV ratio for a Functional Class, falling back to the
    statewide "All Stations" ratio for Local/Unpaved (neither has its own
    row) or any unrecognized/ambiguous class."""
    if functional_class is None or functional_class not in reclass_ratio.index:
        return all_stations_ratio
    row = reclass_ratio.loc[functional_class]
    return row if pd.notna(row["MD Ratio (2025/2003)"]) else all_stations_ratio

Apply: override matched stations, rescale unmatched stations

for ext_id, pct in matched_pct.items():
    mask = df_external_year["WF_Ext"] == ext_id
    df_external_year.loc[mask, "PctTrk_SU"] = pct["MD"]
    df_external_year.loc[mask, "PctTrk_MU"] = pct["HV"]

for ext_id in df_external_year["WF_Ext"].unique():
    if ext_id in matched_pct:
        continue
    ratio = ratio_for_class(functional_class_by_external.get(ext_id))
    mask = df_external_year["WF_Ext"] == ext_id
    df_external_year.loc[mask, "PctTrk_SU"] *= ratio["MD Ratio (2025/2003)"]
    df_external_year.loc[mask, "PctTrk_MU"] *= ratio["HV Ratio (2025/2003)"]

df_external_year["PctTrk_SU"] = df_external_year["PctTrk_SU"].astype(float).round(4)
df_external_year["PctTrk_MU"] = df_external_year["PctTrk_MU"].astype(float).round(4)

# Recompute TRUCK_SU/TRUCK_MU from the (possibly overridden/rescaled) percentages,
# for every year -- the percentage is now the single source of truth throughout.
df_external_year["TRUCK_SU"] = (
    (df_external_year["AADT"] * df_external_year["PctTrk_SU"]).round().astype("Int64")
)
df_external_year["TRUCK_MU"] = (
    (df_external_year["AADT"] * df_external_year["PctTrk_MU"]).round().astype("Int64")
)
df_external_year
WF_Ext Year ;Idx_WF segid route milepost Ext_Name AWDT_FAC AADT_Historic PctTrk_SU PctTrk_MU AADT_Forecast AADT TRUCK_SU TRUCK_MU
0 3601 1981 36011981 1082_000.0 1082PM 0.0 Ext # 3601 - FAR-1082 Bird Refuge 1.0000 0.0 0.0 0.0 NaN 0 0 0
1 3601 1982 36011982 1082_000.0 1082PM 0.0 Ext # 3601 - FAR-1082 Bird Refuge 1.0000 0.0 0.0 0.0 NaN 0 0 0
2 3601 1983 36011983 1082_000.0 1082PM 0.0 Ext # 3601 - FAR-1082 Bird Refuge 1.0000 0.0 0.0 0.0 NaN 0 0 0
3 3601 1984 36011984 1082_000.0 1082PM 0.0 Ext # 3601 - FAR-1082 Bird Refuge 1.0000 0.0 0.0 0.0 NaN 0 0 0
4 3601 1985 36011985 1082_000.0 1082PM 0.0 Ext # 3601 - FAR-1082 Bird Refuge 1.0000 0.0 0.0 0.0 NaN 0 0 0
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
2315 3629 2056 36292056 1826_004.9 1826PM 4.9 Ext # 3629 - FAR-1826 South Ridge Farms 1.0404 NaN 0.0 0.0 6280.0 6280 0 0
2316 3629 2057 36292057 1826_004.9 1826PM 4.9 Ext # 3629 - FAR-1826 South Ridge Farms 1.0404 NaN 0.0 0.0 6360.0 6360 0 0
2317 3629 2058 36292058 1826_004.9 1826PM 4.9 Ext # 3629 - FAR-1826 South Ridge Farms 1.0404 NaN 0.0 0.0 6440.0 6440 0 0
2318 3629 2059 36292059 1826_004.9 1826PM 4.9 Ext # 3629 - FAR-1826 South Ridge Farms 1.0404 NaN 0.0 0.0 6520.0 6520 0 0
2319 3629 2060 36292060 1826_004.9 1826PM 4.9 Ext # 3629 - FAR-1826 South Ridge Farms 1.0404 NaN 0.0 0.0 6600.0 6600 0 0

2320 rows × 15 columns

Step 6: Passenger & weekly counts

df_external_year["PASSENGER"] = (
    (
        df_external_year["AADT"]
        - (
            df_external_year["TRUCK_SU"].fillna(0)
            + df_external_year["TRUCK_MU"].fillna(0)
        )
    )
    .round()
    .astype("Int64")
)

df_external_year["AWDT"] = (
    (df_external_year["AADT"] * df_external_year["AWDT_FAC"]).round().astype("Int64")
)
df_external_year["PASS_VOL"] = (
    (df_external_year["PASSENGER"] * df_external_year["AWDT_FAC"])
    .round()
    .astype("Int64")
)
df_external_year["TRUCK_MD"] = (
    (df_external_year["TRUCK_SU"] * df_external_year["AWDT_FAC"])
    .round()
    .astype("Int64")
)
df_external_year["TRUCK_HV"] = (
    (df_external_year["TRUCK_MU"] * df_external_year["AWDT_FAC"])
    .round()
    .astype("Int64")
)
df_external_year
WF_Ext Year ;Idx_WF segid route milepost Ext_Name AWDT_FAC AADT_Historic PctTrk_SU PctTrk_MU AADT_Forecast AADT TRUCK_SU TRUCK_MU PASSENGER AWDT PASS_VOL TRUCK_MD TRUCK_HV
0 3601 1981 36011981 1082_000.0 1082PM 0.0 Ext # 3601 - FAR-1082 Bird Refuge 1.0000 0.0 0.0 0.0 NaN 0 0 0 0 0 0 0 0
1 3601 1982 36011982 1082_000.0 1082PM 0.0 Ext # 3601 - FAR-1082 Bird Refuge 1.0000 0.0 0.0 0.0 NaN 0 0 0 0 0 0 0 0
2 3601 1983 36011983 1082_000.0 1082PM 0.0 Ext # 3601 - FAR-1082 Bird Refuge 1.0000 0.0 0.0 0.0 NaN 0 0 0 0 0 0 0 0
3 3601 1984 36011984 1082_000.0 1082PM 0.0 Ext # 3601 - FAR-1082 Bird Refuge 1.0000 0.0 0.0 0.0 NaN 0 0 0 0 0 0 0 0
4 3601 1985 36011985 1082_000.0 1082PM 0.0 Ext # 3601 - FAR-1082 Bird Refuge 1.0000 0.0 0.0 0.0 NaN 0 0 0 0 0 0 0 0
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
2315 3629 2056 36292056 1826_004.9 1826PM 4.9 Ext # 3629 - FAR-1826 South Ridge Farms 1.0404 NaN 0.0 0.0 6280.0 6280 0 0 6280 6534 6534 0 0
2316 3629 2057 36292057 1826_004.9 1826PM 4.9 Ext # 3629 - FAR-1826 South Ridge Farms 1.0404 NaN 0.0 0.0 6360.0 6360 0 0 6360 6617 6617 0 0
2317 3629 2058 36292058 1826_004.9 1826PM 4.9 Ext # 3629 - FAR-1826 South Ridge Farms 1.0404 NaN 0.0 0.0 6440.0 6440 0 0 6440 6700 6700 0 0
2318 3629 2059 36292059 1826_004.9 1826PM 4.9 Ext # 3629 - FAR-1826 South Ridge Farms 1.0404 NaN 0.0 0.0 6520.0 6520 0 0 6520 6783 6783 0 0
2319 3629 2060 36292060 1826_004.9 1826PM 4.9 Ext # 3629 - FAR-1826 South Ridge Farms 1.0404 NaN 0.0 0.0 6600.0 6600 0 0 6600 6867 6867 0 0

2320 rows × 20 columns

Step 7: Vintage labeling

df_external_year["Vintage"] = np.where(
    df_external_year["Year"]
    <= df_external_year["WF_Ext"].map(
        df_external_year.dropna(subset=["AADT_Historic"])
        .groupby("WF_Ext")["Year"]
        .max()
    ),
    "Historic",
    "Forecast",
)
df_external_year
WF_Ext Year ;Idx_WF segid route milepost Ext_Name AWDT_FAC AADT_Historic PctTrk_SU ... AADT_Forecast AADT TRUCK_SU TRUCK_MU PASSENGER AWDT PASS_VOL TRUCK_MD TRUCK_HV Vintage
0 3601 1981 36011981 1082_000.0 1082PM 0.0 Ext # 3601 - FAR-1082 Bird Refuge 1.0000 0.0 0.0 ... NaN 0 0 0 0 0 0 0 0 Historic
1 3601 1982 36011982 1082_000.0 1082PM 0.0 Ext # 3601 - FAR-1082 Bird Refuge 1.0000 0.0 0.0 ... NaN 0 0 0 0 0 0 0 0 Historic
2 3601 1983 36011983 1082_000.0 1082PM 0.0 Ext # 3601 - FAR-1082 Bird Refuge 1.0000 0.0 0.0 ... NaN 0 0 0 0 0 0 0 0 Historic
3 3601 1984 36011984 1082_000.0 1082PM 0.0 Ext # 3601 - FAR-1082 Bird Refuge 1.0000 0.0 0.0 ... NaN 0 0 0 0 0 0 0 0 Historic
4 3601 1985 36011985 1082_000.0 1082PM 0.0 Ext # 3601 - FAR-1082 Bird Refuge 1.0000 0.0 0.0 ... NaN 0 0 0 0 0 0 0 0 Historic
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
2315 3629 2056 36292056 1826_004.9 1826PM 4.9 Ext # 3629 - FAR-1826 South Ridge Farms 1.0404 NaN 0.0 ... 6280.0 6280 0 0 6280 6534 6534 0 0 Forecast
2316 3629 2057 36292057 1826_004.9 1826PM 4.9 Ext # 3629 - FAR-1826 South Ridge Farms 1.0404 NaN 0.0 ... 6360.0 6360 0 0 6360 6617 6617 0 0 Forecast
2317 3629 2058 36292058 1826_004.9 1826PM 4.9 Ext # 3629 - FAR-1826 South Ridge Farms 1.0404 NaN 0.0 ... 6440.0 6440 0 0 6440 6700 6700 0 0 Forecast
2318 3629 2059 36292059 1826_004.9 1826PM 4.9 Ext # 3629 - FAR-1826 South Ridge Farms 1.0404 NaN 0.0 ... 6520.0 6520 0 0 6520 6783 6783 0 0 Forecast
2319 3629 2060 36292060 1826_004.9 1826PM 4.9 Ext # 3629 - FAR-1826 South Ridge Farms 1.0404 NaN 0.0 ... 6600.0 6600 0 0 6600 6867 6867 0 0 Forecast

2320 rows × 21 columns

Export final results

Path("results").mkdir(parents=True, exist_ok=True)

(
    df_external_year[
        [
            ";Idx_WF",
            "WF_Ext",
            "Year",
            "AWDT",
            "PASS_VOL",
            "TRUCK_MD",
            "TRUCK_HV",
            "AWDT_FAC",
            "AADT",
            "PASSENGER",
            "TRUCK_SU",
            "TRUCK_MU",
            "PctTrk_SU",
            "PctTrk_MU",
        ]
    ][
        (df_external_year["Year"] >= YEARS["EXPORT_START"])
        &
        # External 3611 link doesn't currently exist, so it shouldn't be forecasted.
        # We'll have to assert a shift from External 3612 to 3611 using USTM model output from UDOT.
        (df_external_year["WF_Ext"] != 3611)
    ].to_csv("results/external_year_vol.csv", index=False)
)