如何基于pyplot在犹他州地图上绘制州际公路？

我编写了一段绘制美国州和县的代码，现将核心代码简化后展示如下，用于绘制犹他州的各个县。我希望在该地图上绘制州际公路系统，但一直未找到可行的实现方法及数据来源，恳请指点可用的数据及绘制示例。

#! /usr/bin/python3

# Based on this YouTube video and the data from the US Census Bureau:
# https://www.youtube.com/watch?v=X706iRqe1hA
# https://www.census.gov/geographies/mapping-files/time-series/geo/carto-boundary-file.html
#
# Projections:
# https://scitools.org.uk/cartopy/docs/latest/reference/projections.html

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import numpy as np
import cartopy.feature as cfeature
import cartopy.io.shapereader as shpreader
from metpy.cbook import get_test_data


def drawUtah():
    fig = plt.figure()
    ax = fig.add_subplot(1, 1, 1, projection=ccrs.Miller())
    ax.stock_img()

    for rec in shpreader.Reader(get_test_data('us_counties_500k.shp',
            as_file_obj=False)).records():
        fips_code = rec.attributes['GEOID']
        if(fips_code.startswith("49")):
            ax.add_geometries([rec.geometry], crs=ccrs.PlateCarree(),
                edgecolor='black', facecolor="#FFFAFA")

    ax.set_extent( (-114.75, -108.5, 36.5, 42.5) )
    ax.add_feature(cfeature.STATES)

    fig.savefig("./utah.png")


if __name__ == "__main__":
    drawUtah()

更新

提交问题两天后（问题在一天内通过审核），我继续搜索后终于找到了道路数据。我下载了该数据压缩包，解压到代码所在目录，并添加了读取和绘制州际公路的代码，更新后的代码如下。此更新旨在帮助其他遇到同类问题的开发者，我也曾从多年前的解决方案中获益。

#! /usr/bin/python3

# Based on this YouTube video and the data from the US Census Bureau:
# https://www.youtube.com/watch?v=X706iRqe1hA
# https://www.census.gov/geographies/mapping-files/time-series/geo/carto-boundary-file.html
#
# Projections:
# https://scitools.org.uk/cartopy/docs/latest/reference/projections.html

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import numpy as np
import cartopy.feature as cfeature
import cartopy.io.shapereader as shpreader
from metpy.cbook import get_test_data


def drawUtah():
    fig = plt.figure()
    ax = fig.add_subplot(1, 1, 1, projection=ccrs.Miller())
    ax.stock_img()

    for rec in shpreader.Reader(get_test_data('us_counties_500k.shp',
            as_file_obj=False)).records():
        fips_code = rec.attributes['GEOID']
        if(fips_code.startswith("49")):
            ax.add_geometries([rec.geometry], crs=ccrs.PlateCarree(),
                edgecolor='black', facecolor="#FFFAFA")

    for rec in shpreader.Reader("./tl_2017_us_primaryroads.shp").records():
        if(rec.attributes["FULLNAME"].startswith("I- ")):
            ax.add_geometries([rec.geometry], crs=ccrs.PlateCarree(),
                edgecolor='red', facecolor='#FF000000')

    ax.set_extent( (-114.75, -108.5, 36.5, 42.5) )
    ax.add_feature(cfeature.STATES)

    fig.savefig("./utah.png")


if __name__ == "__main__":
    drawUtah()

注：新增代码为第二个for循环。

内容的提问来源于stack exchange，提问作者Val