diff --git a/postreise/plot/demo/plot_lmp.ipynb b/postreise/plot/demo/plot_lmp.ipynb
index 5071e76c..d7da9ae3 100644
--- a/postreise/plot/demo/plot_lmp.ipynb
+++ b/postreise/plot/demo/plot_lmp.ipynb
@@ -19,7 +19,7 @@
    "outputs": [],
    "source": [
     "from powersimdata.scenario.scenario import Scenario  \n",
-    "from postreise.plot.plot_lmp import plot_lmp"
+    "from postreise.plot.plot_lmp import map_lmp"
    ]
   },
   {
@@ -446,7 +446,7 @@
     }
    ],
    "source": [
-    "plot_lmp(s_grid,  lmp)"
+    "map_lmp(s_grid,  lmp)"
    ]
   }
  ],
diff --git a/postreise/plot/plot_carbon_map.py b/postreise/plot/plot_carbon_map.py
index fab18a5c..d92dc14c 100644
--- a/postreise/plot/plot_carbon_map.py
+++ b/postreise/plot/plot_carbon_map.py
@@ -49,7 +49,9 @@ def get_borders(us_states_dat, state_list=None):
     return a, b
 
 
-def plot_states(state_list, col_list, labels_list, us_states_dat=us_states.data):
+def plot_states(
+    state_list, col_list, labels_list, font_size, us_states_dat=us_states.data
+):
     """Plots US state borders and allows color coding by state,
         for example to represent different emissions goals.
 
@@ -97,7 +99,7 @@ def plot_states(state_list, col_list, labels_list, us_states_dat=us_states.data)
             y=(max(b1[0]) + min(b1[0])) / 2,
             x_units="data",
             y_units="data",
-            text_font_size="20pt",
+            text_font_size=font_size,
             text=labels_list[n],
             render_mode="css",
             border_line_color="black",
@@ -106,8 +108,7 @@ def plot_states(state_list, col_list, labels_list, us_states_dat=us_states.data)
             background_fill_alpha=0.8,
         )
         p.add_layout(citation)
-    output_notebook()
-    show(p)
+    return p
 
 
 def group_lat_lon(bus_map):
@@ -323,8 +324,7 @@ def map_carbon_emission_bar(
     )
     p_legend.add_layout(labels)
 
-    output_notebook()
-    show(row(p_legend, p))
+    return row(p_legend, p)
 
 
 def map_carbon_emission(
@@ -335,6 +335,7 @@ def map_carbon_emission(
     label_coal="Coal: tons",
     label_ng="NG: tons",
     us_states_dat=us_states.data,
+    size_factor=1,
 ):
     """Makes map of carbon emissions, color code by fuel type. Size/area
         indicates emissions.
@@ -357,8 +358,8 @@ def map_carbon_emission(
         {
             "x": bus_map["x"],
             "y": bus_map["y"],
-            "coal": (bus_map["coal"] / 1000) ** 0.5,
-            "ng": (bus_map["ng"] / 1000) ** 0.5,
+            "coal": (bus_map["coal"] / 1000 * size_factor) ** 0.5,
+            "ng": (bus_map["ng"] / 1000 * size_factor) ** 0.5,
         }
     )
 
@@ -371,6 +372,7 @@ def map_carbon_emission(
         y_axis_location=None,
         plot_width=800,
         plot_height=800,
+        output_backend="webgl",
     )
     p_legend = figure(
         x_axis_location=None,
@@ -380,6 +382,7 @@ def map_carbon_emission(
         plot_height=400,
         y_range=(0, 4),
         x_range=(0, 2),
+        output_backend="webgl",
     )
     p.add_tile(get_provider(Vendors.CARTODBPOSITRON_RETINA))
     # state borders
@@ -416,9 +419,9 @@ def map_carbon_emission(
         color=np.repeat([color_coal, color_ng], 3),
         alpha=0.25,
         size=[
-            (10000000 / 1000) ** 0.5,
-            (5000000 / 1000) ** 0.5,
-            (1000000 / 1000) ** 0.5,
+            (10000000 / 1000 * size_factor) ** 0.5,
+            (5000000 / 1000 * size_factor) ** 0.5,
+            (1000000 / 1000 * size_factor) ** 0.5,
         ]
         * 2,
     )
@@ -441,8 +444,7 @@ def map_carbon_emission(
     )
     p_legend.add_layout(labels)
 
-    output_notebook()
-    show(row(p_legend, p))
+    return row(p_legend, p)
 
 
 def map_carbon_emission_comparison(bus_info_and_emission_1, bus_info_and_emission_2):
diff --git a/postreise/plot/plot_interconnection_map.py b/postreise/plot/plot_interconnection_map.py
new file mode 100644
index 00000000..22fadb24
--- /dev/null
+++ b/postreise/plot/plot_interconnection_map.py
@@ -0,0 +1,130 @@
+from postreise.plot.projection_helpers import project_branch
+from bokeh.models import ColumnDataSource
+from bokeh.plotting import figure
+from bokeh.sampledata import us_states
+from postreise.plot import plot_carbon_map
+from bokeh.tile_providers import get_provider, Vendors
+
+
+def map_interconnections(grid, us_states_dat=us_states.data):
+    """Makes map of transmission color coded by interconnection
+
+    :param grid: grid object
+    :param dict us_states_dat: us_states data file, imported from bokeh.
+    :return:  -- map of transmission
+    """
+    # projection steps for mapping
+    branch = grid.branch
+    branch_bus = grid.bus
+    branch_map = project_branch(branch)
+    branch_west = branch_map.loc[branch_map.interconnect == "Western"]
+    branch_east = branch_map.loc[branch_map.interconnect == "Eastern"]
+    branch_tx = branch_map.loc[branch_map.interconnect == "Texas"]
+    branch_mdc = grid.dcline
+
+    branch_mdc["from_lon"] = branch_bus.loc[branch_mdc.from_bus_id, "lon"].values
+    branch_mdc["from_lat"] = branch_bus.loc[branch_mdc.from_bus_id, "lat"].values
+    branch_mdc["to_lon"] = branch_bus.loc[branch_mdc.to_bus_id, "lon"].values
+    branch_mdc["to_lat"] = branch_bus.loc[branch_mdc.to_bus_id, "lat"].values
+    branch_mdc = project_branch(branch_mdc)
+    # state borders
+    a, b = plot_carbon_map.get_borders(us_states_dat.copy())
+
+    # transmission data sources
+    multi_line_source = ColumnDataSource(
+        {
+            "xs": branch_west[["from_x", "to_x"]].values.tolist(),
+            "ys": branch_west[["from_y", "to_y"]].values.tolist(),
+            "capacity": (branch_west.rateA) * 0.000225 + 0.1,
+        }
+    )
+
+    multi_line_source2 = ColumnDataSource(
+        {
+            "xs": branch_east[["from_x", "to_x"]].values.tolist(),
+            "ys": branch_east[["from_y", "to_y"]].values.tolist(),
+            "capacity": (branch_east.rateA) * 0.000225 + 0.1,
+        }
+    )
+
+    multi_line_source3 = ColumnDataSource(
+        {
+            "xs": branch_tx[["from_x", "to_x"]].values.tolist(),
+            "ys": branch_tx[["from_y", "to_y"]].values.tolist(),
+            "capacity": (branch_tx.rateA) * 0.000225 + 0.1,
+        }
+    )
+
+    multi_line_source4 = ColumnDataSource(
+        {
+            "xs": branch_mdc[["from_x", "to_x"]].values.tolist(),
+            "ys": branch_mdc[["from_y", "to_y"]].values.tolist(),
+            "capacity": (branch_mdc.Pmax.astype(float)) * 0.000225 + 0.1,
+        }
+    )
+
+    # Set up figure
+    tools: str = "pan,wheel_zoom,reset,hover,save"
+
+    p = figure(
+        title="Interconnections",
+        tools=tools,
+        x_axis_location=None,
+        y_axis_location=None,
+        plot_width=800,
+        plot_height=800,
+        output_backend="webgl",
+    )
+
+    # for legend
+    p.multi_line(
+        [-1.084288e07, -1.084288e07],
+        [4.639031e06, 4.639031e06],
+        color="blue",
+        line_width=5,
+        legend="Western",
+    )
+    p.multi_line(
+        [-1.084288e07, -1.084288e07],
+        [4.639031e06, 4.639031e06],
+        color="red",
+        line_width=5,
+        legend="Eastern",
+    )
+    p.multi_line(
+        [-1.084288e07, -1.084288e07],
+        [4.639031e06, 4.639031e06],
+        color="purple",
+        line_width=5,
+        legend="ERCOT",
+    )
+    p.multi_line(
+        [-1.084288e07, -1.084288e07],
+        [4.639031e06, 4.639031e06],
+        color="green",
+        line_width=5,
+        legend="HVDC",
+    )
+    p.square(
+        x=[-1.084288e07], y=[4.639031e06], size=170, fill_color="white", color="white"
+    )
+    p.add_tile(get_provider(Vendors.CARTODBPOSITRON))
+
+    # state borders
+    p.patches(a, b, fill_alpha=0.0, line_color="black", line_width=2)
+    # branches
+    p.multi_line(
+        "xs", "ys", color="blue", line_width="capacity", source=multi_line_source
+    )
+    p.multi_line(
+        "xs", "ys", color="red", line_width="capacity", source=multi_line_source2
+    )
+    p.multi_line(
+        "xs", "ys", color="purple", line_width="capacity", source=multi_line_source3
+    )
+    p.multi_line(
+        "xs", "ys", color="green", line_width="capacity", source=multi_line_source4
+    )
+    p.legend.location = "bottom_right"
+
+    return p
diff --git a/postreise/plot/plot_lmp.py b/postreise/plot/plot_lmp.py
index 533edd9c..f68fa814 100644
--- a/postreise/plot/plot_lmp.py
+++ b/postreise/plot/plot_lmp.py
@@ -2,7 +2,6 @@
 #   PostREISE/postreise/plot/demo: plot_lmp.ipynb
 
 import pandas as pd
-from bokeh.io import output_notebook, show
 from bokeh.layouts import row
 from bokeh.models import ColumnDataSource
 from bokeh.plotting import figure
@@ -19,21 +18,23 @@
 default_states_list = list(default_states_dict.keys())
 
 
-def plot_lmp(s_grid, lmp, us_states_dat=us_states.data):
+def map_lmp(s_grid, lmp, us_states_dat=us_states.data):
     """ Plot average LMP by color coding buses
 
     :param s_grid: scenario grid
     :type s_grid: pandas.DataFrame
     :param lmp: lmps (locational marginal prices) calculated for the scenario
     :type lmp: pandas.DataFrame
+    :param file_name: name for output png file
+    :type file_name: str
     :param us_states_dat: us_states data file, imported from bokeh
     :type us_states_dat: dict
     """
 
     bus = project_bus(s_grid.bus)
-    lmp_split_points = list(range(0, 257, 1))
+    lmp_split_points = list(range(0, 256, 1))
     bus_segments = _construct_bus_data(bus, lmp, lmp_split_points)
-    _construct_shadowprice_visuals(lmp_split_points, bus_segments, us_states_dat)
+    return _construct_shadowprice_visuals(lmp_split_points, bus_segments, us_states_dat)
 
 
 def _construct_bus_data(bus_map, lmp, lmp_split_points):
@@ -54,28 +55,20 @@ def _construct_bus_data(bus_map, lmp, lmp_split_points):
 
     # min and max values for 'continuous' color scale, in $MW/h
     min_lmp_clamp = 20
-    max_lmp_clamp = 50
+    max_lmp_clamp = 45
 
     lmp_mean["lmp_norm"] = (
         (lmp_mean.lmp - min_lmp_clamp) / (max_lmp_clamp - min_lmp_clamp) * 256
     )
     bus_map = pd.concat([bus_map, lmp_mean], axis=1)
-    bus_map = group_lat_lon(bus_map)
+    bus_map_agg = group_lat_lon(bus_map)
 
     # set a min and a max so all values are assigned a color
-    bus_map.loc[(bus_map.lmp_norm <= 0), "lmp_norm"] = 0.01
-    bus_map.loc[(bus_map.lmp_norm > 256), "lmp_norm"] = 256
+    bus_map_agg.loc[(bus_map_agg.lmp_norm <= 0), "lmp_norm"] = 0.01
+    bus_map_agg.loc[(bus_map_agg.lmp_norm > 255), "lmp_norm"] = 255
+    bus_map_agg = pd.DataFrame(bus_map_agg)
 
-    bus_segments = []
-
-    for i in range(len(lmp_split_points) - 1):
-        bus_segment = bus_map[
-            (bus_map["lmp_norm"] > lmp_split_points[i])
-            & (bus_map["lmp_norm"] <= lmp_split_points[i + 1])
-        ]
-        bus_segments.append(bus_segment)
-
-    return bus_segments
+    return bus_map_agg
 
 
 def group_lat_lon(bus_map):
@@ -92,11 +85,11 @@ def group_lat_lon(bus_map):
     bus_map1.lat = bus_map1.lat.round(1)
     bus_map1.lon = bus_map1.lon.round(1)
 
-    bus_map1 = bus_map1.groupby(["lat", "lon"]).agg(
+    bus_map2 = bus_map1.groupby(["lat", "lon"]).agg(
         {"lmp_norm": "mean", "lmp": "mean", "x": "mean", "y": "mean"}
     )
 
-    return bus_map1
+    return bus_map2
 
 
 def _construct_shadowprice_visuals(lmp_split_points, bus_segments, us_states_dat):
@@ -112,6 +105,8 @@ def _construct_shadowprice_visuals(lmp_split_points, bus_segments, us_states_dat
     :type bus_segments: list(pandas.DataFrame)
     :param us_states_dat: us_states data file, imported from bokeh
     :type us_states_dat: dict
+    :param file_name: name for output png file
+    :type file_name: str
     """
 
     tools = "pan,wheel_zoom,reset,hover,save"
@@ -121,6 +116,7 @@ def _construct_shadowprice_visuals(lmp_split_points, bus_segments, us_states_dat
         y_axis_location=None,
         plot_width=800,
         plot_height=800,
+        output_backend="webgl",
     )
 
     # Add USA map
@@ -129,22 +125,17 @@ def _construct_shadowprice_visuals(lmp_split_points, bus_segments, us_states_dat
     a, b = get_borders(us_states_dat.copy())
     p.patches(a, b, fill_alpha=0.0, line_color="black", line_width=2)
     # Add colored circles for bus locations
-    indices = list(range(len(bus_segments)))
-    indices.reverse()  # We want the lowest prices on top
-    for i in indices:
-        bus_cds = ColumnDataSource(
-            {
-                "x": bus_segments[i]["x"],
-                "y": bus_segments[i]["y"],
-                "lmp": bus_segments[i]["lmp_norm"],
-            }
-        )
-        p.circle("x", "y", color=Turbo256[i], alpha=0.7, size=2, source=bus_cds)
+
+    bus_segments.lmp_norm = bus_segments.lmp_norm.astype(int)
+    bus_segments["col"] = bus_segments["lmp_norm"].apply(lambda x: Turbo256[x])
+    bus_cds = ColumnDataSource(
+        {"x": bus_segments["x"], "y": bus_segments["y"], "col": bus_segments["col"]}
+    )
+    p.circle("x", "y", color="col", size=2, source=bus_cds)
 
     # Add legend
     bus_legend = _construct_bus_legend(lmp_split_points)
-    output_notebook()
-    show(row(bus_legend, p))
+    return row(bus_legend, p)
 
 
 def _construct_bus_legend(lmp_split_points):
@@ -168,6 +159,7 @@ def _construct_bus_legend(lmp_split_points):
         plot_width=110,
         toolbar_location=None,
         tools="",
+        output_backend="webgl",
     )
     p.vbar_stack(
         list(bars.keys())[1:],
@@ -185,7 +177,7 @@ def _construct_bus_legend(lmp_split_points):
     p.yaxis.ticker = list(labels.keys())
     p.yaxis.major_label_overrides = labels
 
-    p.text(x=[-0.05], y=[bar_len_sum * 1.01], text=[f"$/MWh"], text_font_size="12pt")
+    p.text(x=[-0.05], y=[bar_len_sum * 1.01], text=[f" $/MWh"], text_font_size="12pt")
 
     return p
 
@@ -205,9 +197,9 @@ def _get_bus_legend_bars_and_labels(lmp_split_points, x_range):
     labels = {}  # { y-position: label_text, ... }
     for i in range(len(lmp_split_points) - 2):
         bar_length = 0
-        if i == 0 or i == len(lmp_split_points) - 2:
-            # For first and last bars, clamp bar length between 1 and 5
-            # to prevent extreme min and max vals from overwhelming the legend
+        if i == 0 or i == len(lmp_split_points) - 1:
+            # For first and last bars, clamp bar length
+            # to prevent extreme min, max vals from overwhelming the legend
             lmp_diff = lmp_split_points[i + 1] - lmp_split_points[i]
             bar_length = 1 if lmp_diff < 1 else 5 if lmp_diff > 5 else lmp_diff
         elif i != len(lmp_split_points) - 1:
@@ -217,7 +209,7 @@ def _get_bus_legend_bars_and_labels(lmp_split_points, x_range):
             bar_length = max(1, round(lmp_diff, 1))
 
         labels[round(bar_length_sum, -1)] = str(
-            round((lmp_split_points[i] * 30 / 256 + 20), 0)
+            round(((lmp_split_points[i] - 1) * (45 - 20) / 256 + 20), 0)
         )
         if i != len(lmp_split_points) - 1:
             bars[str(i)] = [bar_length]
diff --git a/postreise/plot/plot_utilization_map.py b/postreise/plot/plot_utilization_map.py
index 8fdc8091..8eca3b8f 100644
--- a/postreise/plot/plot_utilization_map.py
+++ b/postreise/plot/plot_utilization_map.py
@@ -1,9 +1,10 @@
 import pandas as pd
 from bokeh.models import ColumnDataSource, ColorBar
-from bokeh.plotting import figure, show
+from bokeh.plotting import figure
 from bokeh.sampledata import us_states
 from bokeh.tile_providers import get_provider, Vendors
 from bokeh.transform import linear_cmap
+import numpy as np
 
 from postreise.plot import plot_carbon_map
 from postreise.plot.projection_helpers import project_branch
@@ -24,9 +25,10 @@
 ]
 
 
-def map_risk_bind(congestion_stats, branch, us_states_dat=us_states.data):
-    """Makes map showing risk and binding incidents on US states map.
+def map_risk_bind(risk_or_bind, congestion_stats, branch, us_states_dat=us_states.data):
+    """Makes map showing risk or binding incidents on US states map.
 
+    :param str risk_or_bind: specify plotting "risk" or "bind"
     :param pandas.DataFrame congestion_stats: data frame as returned by
         :func:`postreise.analyze.transmission.generate_cong_stats`.
     :param pandas.DataFrame branch: branch data frame.
@@ -42,58 +44,55 @@ def map_risk_bind(congestion_stats, branch, us_states_dat=us_states.data):
         {
             "xs": branch_map_all[["from_x", "to_x"]].values.tolist(),
             "ys": branch_map_all[["from_y", "to_y"]].values.tolist(),
-            "cap": branch_map_all["rateA"] / 1000 + 1,
+            "cap": branch_map_all["rateA"] / 1000 + 2,
         }
     )
     a, b = plot_carbon_map.get_borders(us_states_dat.copy())
     tools: str = "pan,wheel_zoom,reset,hover,save"
 
-    # make two plots for risk and bind respectively
-    for param in ["risk", "bind"]:
-        branch_congestion1 = branch_congestion[branch_congestion[param] > 0]
-        branch_congestion1.sort_values(by=[param])
-        branch_map = project_branch(branch_congestion1)
-        min_val = branch_congestion1[param].min()
-        max_val = branch_congestion1[param].max()
-        mapper1 = linear_cmap(
-            field_name=param, palette=traffic_palette, low=min_val, high=max_val
-        )
-        color_bar = ColorBar(
-            color_mapper=mapper1["transform"], width=8, location=(0, 0), title=param
-        )
-        multi_line_source = ColumnDataSource(
-            {
-                "xs": branch_map[["from_x", "to_x"]].values.tolist(),
-                "ys": branch_map[["from_y", "to_y"]].values.tolist(),
-                param: branch_map[param],
-                "cap": branch_map["capacity"] / 1000 + 2,
-            }
-        )
-
-        # Set up figure
-        p = figure(
-            title=param,
-            tools=tools,
-            x_axis_location=None,
-            y_axis_location=None,
-            plot_width=800,
-            plot_height=800,
-        )
-        p.add_layout(color_bar, "right")
-        p.add_tile(get_provider(Vendors.CARTODBPOSITRON))
-        p.patches(a, b, fill_alpha=0.0, line_color="black", line_width=2)
-        p.multi_line(
-            "xs",
-            "ys",
-            color="gray",
-            line_width="cap",
-            source=multi_line_source_all,
-            alpha=0.5,
-        )
-        p.multi_line(
-            "xs", "ys", color=mapper1, line_width="cap", source=multi_line_source
-        )
-        show(p)
+    branch_congestion = branch_congestion[branch_congestion[risk_or_bind] > 0]
+    branch_congestion.sort_values(by=[risk_or_bind])
+    branch_map = project_branch(branch_congestion)
+    min_val = branch_congestion[risk_or_bind].min()
+    max_val = branch_congestion[risk_or_bind].max()
+    mapper = linear_cmap(
+        field_name=risk_or_bind, palette=traffic_palette, low=min_val, high=max_val
+    )
+    color_bar = ColorBar(
+        color_mapper=mapper["transform"], width=8, location=(0, 0), title=risk_or_bind
+    )
+    multi_line_source = ColumnDataSource(
+        {
+            "xs": branch_map[["from_x", "to_x"]].values.tolist(),
+            "ys": branch_map[["from_y", "to_y"]].values.tolist(),
+            risk_or_bind: branch_map[risk_or_bind],
+            "cap": branch_map["capacity"] / 1000 + 2,
+        }
+    )
+
+    # Set up figure
+    p = figure(
+        title=risk_or_bind,
+        tools=tools,
+        x_axis_location=None,
+        y_axis_location=None,
+        plot_width=800,
+        plot_height=800,
+        output_backend="webgl",
+    )
+    p.add_layout(color_bar, "right")
+    p.add_tile(get_provider(Vendors.CARTODBPOSITRON))
+    p.patches(a, b, fill_alpha=0.0, line_color="black", line_width=2)
+    p.multi_line(
+        "xs",
+        "ys",
+        color="gray",
+        line_width="cap",
+        source=multi_line_source_all,
+        alpha=0.5,
+    )
+    p.multi_line("xs", "ys", color=mapper1, line_width="cap", source=multi_line_source)
+    return p
 
 
 def map_utilization(utilization_df, branch, us_states_dat=us_states.data):
@@ -137,6 +136,7 @@ def map_utilization(utilization_df, branch, us_states_dat=us_states.data):
 
     branch_map = project_branch(branch_utilization)
     branch_map = branch_map.sort_values(by=["median_utilization"])
+    branch_map = branch_map[~branch_map.isin([np.nan, np.inf, -np.inf]).any(1)]
 
     # state borders
     a, b = plot_carbon_map.get_borders(us_states_dat.copy())
@@ -160,9 +160,10 @@ def map_utilization(utilization_df, branch, us_states_dat=us_states.data):
         y_axis_location=None,
         plot_width=800,
         plot_height=800,
+        output_backend="webgl",
     )
     p.add_layout(color_bar, "right")
     p.add_tile(get_provider(Vendors.CARTODBPOSITRON))
     p.patches(a, b, fill_alpha=0.0, line_color="black", line_width=2)
     p.multi_line("xs", "ys", color=mapper1, line_width="cap", source=multi_line_source)
-    show(p)
+    return p