diff --git a/README.rst b/README.rst
index 0fa419501..bda26fd43 100644
--- a/README.rst
+++ b/README.rst
@@ -57,7 +57,7 @@ as your package/environment manager.
     From your home directory ``/home/{user}/`` or another directory that you have permissions in, run the command ``git clone git@github.com:NREL/reVX.git`` and then go into your cloned repository: ``cd reVX``
 
 #. Install reVX:
-    1) Follow the installation commands installation process that we use for our automated test suite `here <https://github.com/NREL/reVX/blob/main/.github/workflows/pull_request_tests.yml#L31-L36>`_. Make sure that you call ``pip install -e .`` from within the cloned repository directory e.g. ``/home/{user}/reVX/``
+    1) Follow the installation commands installation process that we use for our automated test suite `here <https://github.com/NREL/reVX/blob/main/.github/workflows/pull_request_tests.yml#L31-L34>`_. Make sure that you call ``pip install -e .`` from within the cloned repository directory e.g. ``/home/{user}/reVX/``
 
         - NOTE: If you install using pip and want to run `exclusion setbacks <https://nrel.github.io/reVX/_cli/reVX.setbacks.setbacks.html>`_ you will need to install rtree manually:
             * ``conda install rtree``
diff --git a/reVX/config/least_cost_xmission.py b/reVX/config/least_cost_xmission.py
index 800a7b16e..8f320d83c 100644
--- a/reVX/config/least_cost_xmission.py
+++ b/reVX/config/least_cost_xmission.py
@@ -222,6 +222,14 @@ def barrier_mult(self):
         """
         return self.get('barrier_mult', self._default_barrier_mult)
 
+    @property
+    def allow_connections_within_states(self):
+        """
+        Boolean flag to allow supple curve points to connect to
+        substations outside their BA but within their own state.
+        """
+        return self.get("allow_connections_within_states", False)
+
 
 class LeastCostPathsConfig(AnalysisConfig):
     """Config framework for Least Cost Paths"""
@@ -306,6 +314,14 @@ def barrier_mult(self):
         """
         return self.get('barrier_mult', self._default_barrier_mult)
 
+    @property
+    def allow_connections_within_states(self):
+        """
+        Boolean flag to allow substations to connect to endpoints
+        outside their BA but within their own state.
+        """
+        return self.get("allow_connections_within_states", False)
+
     @property
     def save_paths(self):
         """
diff --git a/reVX/least_cost_xmission/README.md b/reVX/least_cost_xmission/README.md
index 18861f839..6f95e44d4 100644
--- a/reVX/least_cost_xmission/README.md
+++ b/reVX/least_cost_xmission/README.md
@@ -107,6 +107,7 @@ In this methodology, total interconnection costs are comprised of two components
 - Each Balancing Area has exactly one "Network Node" (typically a city or highly populated area)
 - SC points **may only connect to substations**
 - Substations that a SC point connects to **must be in the same Balancing Area as the SC point**
+    - This assumption can be relaxed to allow connections within the same state.
 - Reinforcement costs are calculated based on the distance between the substation a SC point connected to and the Network Node in that Balancing Area
     - The path used to calculate reinforcement costs is traced along existing transmission lines **for as long as possible**.
     - The reinforcement cost is taken to be half (50%) of the total greenfield cost of the transmission line being traced. If a reinforcement path traces along multiple transmission lines, the corresponding greenfield costs are used for each segment. If multiple transmission lines are available in a single raster pixel, the cost for the highest-voltage line is used. Wherever there is no transmission line, a default greenfield cost assumption (specified by the user; typically 230 kV) is used.
@@ -149,7 +150,7 @@ Next, compute the reinforcement paths on multiple nodes. Use the file below as a
 }
 ```
 
-Note that we are specifying ``"capacity_class": "400"``  to use the 230 kV (400MW capacity) greenfield costs for portions of the reinforcement paths that do no have existing transmission. If you would like to save the reinforcement path geometries, simply add `"save_paths": true` to the file, but note that this may increase your data product size significantly.
+Note that we are specifying ``"capacity_class": "400"``  to use the 230 kV (400MW capacity) greenfield costs for portions of the reinforcement paths that do no have existing transmission. If you would like to save the reinforcement path geometries, simply add `"save_paths": true` to the file, but note that this may increase your data product size significantly. If you would like to allow substations to connect to endpoints within the same state, add `"allow_connections_within_states": true` to the file.
 
 After putting together your config file, simply call
 
@@ -183,6 +184,7 @@ You should now have a file containing all of the reinforcement costs for the sub
     "name": "least_cost_transmission_1000MW"
 }
 ```
+If you would like to allow supply curve points to  connect to substations within the same state, add `"allow_connections_within_states": true` to the file.
 
 Kickoff the execution using the following command:
 
diff --git a/reVX/least_cost_xmission/least_cost_paths.py b/reVX/least_cost_xmission/least_cost_paths.py
index 117c34959..116477fd9 100644
--- a/reVX/least_cost_xmission/least_cost_paths.py
+++ b/reVX/least_cost_xmission/least_cost_paths.py
@@ -249,7 +249,7 @@ def end_features(self):
 
         Returns
         -------
-        geopandas.GeoDataFrame
+        pandas.DataFrame
         """
         end_features = self._features.drop(index=self._start_feature_ind)
         end_features['start_index'] = self._start_feature_ind
@@ -325,7 +325,7 @@ def process_least_cost_paths(self, capacity_class, barrier_mult=100,
                     end_features = end_features.drop(columns=['row', 'col'],
                                                      errors="ignore")
                     lcp = future.result()
-                    lcp = pd.concat((end_features, lcp), axis=1)
+                    lcp = pd.concat((lcp, end_features), axis=1)
                     least_cost_paths.append(lcp)
                     logger.debug('Least cost path {} of {} complete!'
                                  .format(i + 1, len(futures)))
@@ -344,7 +344,7 @@ def process_least_cost_paths(self, capacity_class, barrier_mult=100,
                                        save_paths=save_paths)
                 end_features = self.end_features.drop(columns=['row', 'col'],
                                                       errors="ignore")
-                lcp = pd.concat((end_features, lcp), axis=1)
+                lcp = pd.concat((lcp, end_features), axis=1)
                 least_cost_paths.append(lcp)
 
                 logger.debug('Least cost path {} of {} complete!'
@@ -534,14 +534,15 @@ def process_least_cost_paths(self, capacity_class, barrier_mult=100,
                                          barrier_mult=barrier_mult,
                                          save_paths=save_paths)
         feats = self._features.drop(columns=['row', 'col'])
-        least_cost_paths = pd.concat((feats, lcp), axis=1)
+        least_cost_paths = pd.concat((lcp, feats), axis=1)
 
         return least_cost_paths.drop("index", axis="columns", errors="ignore")
 
     @classmethod
     def run(cls, cost_fpath, features_fpath, network_nodes_fpath,
             transmission_lines_fpath, capacity_class, xmission_config=None,
-            barrier_mult=100, indices=None, save_paths=False):
+            barrier_mult=100, indices=None,
+            allow_connections_within_states=False, save_paths=False):
         """
         Find the reinforcement line paths between the network node and
         the substations for the given tie-line capacity class
@@ -580,6 +581,10 @@ def run(cls, cost_fpath, features_fpath, network_nodes_fpath,
         max_workers : int, optional
             Number of workers to use for processing. If 1 run in serial,
             if ``None`` use all available cores. By default, ``None``.
+        allow_connections_within_states : bool, optional
+            Allow substations to connect to network nodes outside of
+            their own BA, as long as all connections stay within the
+            same state. By default, ``False``.
         save_paths : bool, optional
             Flag to save reinforcement line path as a multi-line
             geometry. By default, ``False``.
@@ -621,7 +626,15 @@ def run(cls, cost_fpath, features_fpath, network_nodes_fpath,
             network_node = (network_nodes.iloc[index:index + 1]
                             .reset_index(drop=True))
             ba_str = network_node["ba_str"].values[0]
-            node_substations = substations[substations["ba_str"] == ba_str]
+            if allow_connections_within_states:
+                state_nodes = network_nodes[network_nodes["state"]
+                                            == network_node["state"].values[0]]
+                allowed_bas = set(state_nodes["ba_str"])
+            else:
+                allowed_bas = {ba_str}
+
+            node_substations = substations[substations["ba_str"]
+                                           .isin(allowed_bas)]
             node_substations = node_substations.reset_index(drop=True)
             logger.debug('Working on {} substations in BA {}'
                          .format(len(node_substations), ba_str))
@@ -638,7 +651,8 @@ def run(cls, cost_fpath, features_fpath, network_nodes_fpath,
         logger.info('{} paths were computed in {:.4f} hours'
                     .format(len(least_cost_paths), (time.time() - ts) / 3600))
 
-        return pd.concat(least_cost_paths, ignore_index=True)
+        costs = pd.concat(least_cost_paths, ignore_index=True)
+        return min_reinforcement_costs(costs)
 
 
 def _rasterize_transmission(transmission_lines, xmission_config, cost_shape,
@@ -675,3 +689,26 @@ def _rasterize_transmission_layer(transmission_lines, cost_shape,
                                 fill=0, transform=cost_transform)
 
     return out
+
+
+def min_reinforcement_costs(table):
+    """Filter table down to cheapest reinforcement per substation.
+
+    Parameters
+    ----------
+    table : pd.DataFrame | gpd.GeoDataFrame
+        Table containing costs for reinforced transmission. Must contain
+        a `gid` column identifying each substation with its own unique
+        ID and a `reinforcement_cost_per_mw` column with the
+        reinforcement costs to minimize.
+
+    Returns
+    -------
+    pd.DataFrame | gpd.GeoDataFrame
+        Table with a single entry for each `gid` with the least
+        `reinforcement_cost_per_mw`.
+    """
+
+    grouped = table.groupby('gid')
+    table = table.loc[grouped["reinforcement_cost_per_mw"].idxmin()]
+    return table.reset_index(drop=True)
diff --git a/reVX/least_cost_xmission/least_cost_paths_cli.py b/reVX/least_cost_xmission/least_cost_paths_cli.py
index 8cecf43a5..40bb85a89 100644
--- a/reVX/least_cost_xmission/least_cost_paths_cli.py
+++ b/reVX/least_cost_xmission/least_cost_paths_cli.py
@@ -74,6 +74,7 @@ def run_local(ctx, config):
                start_index=0, step_index=1,
                barrier_mult=config.barrier_mult,
                max_workers=config.execution_control.max_workers,
+               state_connections=config.allow_connections_within_states,
                save_paths=config.save_paths,
                out_dir=config.dirout,
                log_dir=config.log_directory,
@@ -163,6 +164,9 @@ def from_config(ctx, config, verbose):
               show_default=True, default=None,
               help=("Number of workers to use for processing, if 1 run in "
                     "serial, if None use all available cores"))
+@click.option('--state_connections', '-acws', is_flag=True,
+              help='Flag to allow substations ot connect to any endpoints '
+                   'within their state. Default is not verbose.')
 @click.option('--save_paths', '-paths', is_flag=True,
               help="Flag to save least cost path as a multi-line geometry")
 @click.option('--out_dir', '-o', type=STR, default='./',
@@ -176,7 +180,8 @@ def from_config(ctx, config, verbose):
 @click.pass_context
 def local(ctx, cost_fpath, features_fpath, capacity_class, network_nodes_fpath,
           transmission_lines_fpath, xmission_config, start_index, step_index,
-          barrier_mult, max_workers, save_paths, out_dir, log_dir, verbose):
+          barrier_mult, max_workers, state_connections, save_paths, out_dir,
+          log_dir, verbose):
     """
     Run Least Cost Paths on local hardware
     """
@@ -197,14 +202,16 @@ def local(ctx, cost_fpath, features_fpath, capacity_class, network_nodes_fpath,
         features = gpd.read_file(network_nodes_fpath)
         features, *__ = LeastCostPaths._map_to_costs(cost_fpath, features)
         indices = features.index[start_index::step_index]
+        kwargs = {"xmission_config": xmission_config,
+                  "barrier_mult": barrier_mult,
+                  "indices": indices,
+                  "allow_connections_within_states": state_connections,
+                  "save_paths": save_paths}
         least_costs = ReinforcementPaths.run(cost_fpath, features_fpath,
                                              network_nodes_fpath,
                                              transmission_lines_fpath,
                                              capacity_class,
-                                             xmission_config=xmission_config,
-                                             barrier_mult=barrier_mult,
-                                             indices=indices,
-                                             save_paths=save_paths)
+                                             **kwargs)
     else:
         features = gpd.read_file(features_fpath)
         features, *__ = LeastCostPaths._map_to_costs(cost_fpath, features)
@@ -328,6 +335,9 @@ def get_node_cmd(config, start_index=0):
             '-log {}'.format(SLURM.s(config.log_directory)),
             ]
 
+    if config.allow_connections_within_states:
+        args.append('-acws')
+
     if config.save_paths:
         args.append('-paths')
 
diff --git a/reVX/least_cost_xmission/least_cost_xmission.py b/reVX/least_cost_xmission/least_cost_xmission.py
index 2ad6bd560..834ee3031 100644
--- a/reVX/least_cost_xmission/least_cost_xmission.py
+++ b/reVX/least_cost_xmission/least_cost_xmission.py
@@ -355,7 +355,8 @@ def _clip_to_sc_point(self, sc_point, tie_line_voltage, nn_sinks=2,
         clipping_buffer : float, optional
             Buffer to increase clipping radius by, by default 1.05
         radius : None | int, optional
-            Force clipping radius if set to an int
+            Force clipping radius if set to an int. Radius will be
+            expanded to include at least one connection feature.
 
         Returns
         -------
@@ -375,29 +376,13 @@ def _clip_to_sc_point(self, sc_point, tie_line_voltage, nn_sinks=2,
                 radius = np.abs(self.sink_coords[pos] - np.array([row, col])
                                 ).max()
                 radius = int(np.ceil(radius * clipping_buffer))
-
-            if radius:
-                logger.debug('Using forced radius of {}'.format(radius))
-            else:
                 logger.debug('Radius to {} nearest sink is: {}'
                              .format(nn_sinks, radius))
-            row_min = max(row - radius, 0)
-            row_max = min(row + radius, self._shape[0])
-            col_min = max(col - radius, 0)
-            col_max = min(col + radius, self._shape[1])
-            logger.debug('Extracting all transmission features in the row '
-                         'slice {}:{} and column slice {}:{}'
-                         .format(row_min, row_max, col_min, col_max))
-
-            # Clip transmission features
-            mask = self.features['row'] >= row_min
-            mask &= self.features['row'] < row_max
-            mask &= self.features['col'] >= col_min
-            mask &= self.features['col'] < col_max
-            sc_features = self.features.loc[mask].copy(deep=True)
-            logger.debug('{} transmission features found in clipped area with '
-                         'radius {}'
-                         .format(len(sc_features), radius))
+            else:
+                logger.debug('Using forced radius of {}'.format(radius))
+
+            sc_features = self._clip_to_radius(sc_point, radius, sc_features,
+                                               clipping_buffer)
         else:
             sc_features = self.features.copy(deep=True)
 
@@ -426,6 +411,35 @@ def _clip_to_sc_point(self, sc_point, tie_line_voltage, nn_sinks=2,
 
         return sc_features, radius
 
+    def _clip_to_radius(self, sc_point, radius, sc_features, clipping_buffer):
+        """Clip features to radius.
+
+        If no features are found within the initial radius, it is
+        expanded (multiplicatively by the clipping buffer) until at
+        least one connection feature is found.
+        """
+        if radius is None or len(sc_features) == 0:
+            return sc_features
+
+        # Get pixel resolution and calculate buffer
+        with ExclusionLayers(self._cost_fpath) as ds:
+            resolution = ds.profile["transform"][0]
+        radius_m = radius * resolution
+        logger.debug('Clipping features to radius {}m'.format(radius_m))
+        buffer = sc_point["geometry"].buffer(radius_m)
+        clipped_sc_features = sc_features.clip(buffer)
+
+        while len(clipped_sc_features) <= 0:
+            radius_m *= clipping_buffer
+            logger.debug('Clipping features to radius {}m'.format(radius_m))
+            buffer = sc_point["geometry"].buffer(radius_m)
+            clipped_sc_features = sc_features.clip(buffer)
+
+        logger.debug('{} transmission features found in clipped area with '
+                     'radius {}'
+                     .format(len(clipped_sc_features), radius))
+        return clipped_sc_features.copy(deep=True)
+
     def process_sc_points(self, capacity_class, sc_point_gids=None, nn_sinks=2,
                           clipping_buffer=1.05, barrier_mult=100,
                           max_workers=None, save_paths=False, radius=None,
@@ -456,7 +470,8 @@ def process_sc_points(self, capacity_class, sc_point_gids=None, nn_sinks=2,
             Flag to return least cost paths as a multi-line geometry,
             by default False
         radius : None | int, optional
-            Force clipping radius if set to an int
+            Force clipping radius if set to an int. Radius will be
+            expanded to include at least one connection feature.
         mp_delay : float, optional
             Delay in seconds between starting multi-process workers.
             Useful for reducing memory spike at working startup.
@@ -550,7 +565,8 @@ def _process_multi_core(self, capacity_class, tie_line_voltage,
             Flag to return least cost paths as a multi-line geometry,
             by default False
         radius : None | int, optional
-            Force clipping radius if set to an int
+            Force clipping radius if set to an int. Radius will be
+            expanded to include at least one connection feature.
         mp_delay : float, optional
             Delay in seconds between starting multi-process workers.
             Useful for reducing memory spike at working startup.
@@ -572,7 +588,7 @@ def _process_multi_core(self, capacity_class, tie_line_voltage,
             for _, sc_point in self.sc_points.iterrows():
                 gid = sc_point['sc_point_gid']
                 if gid in sc_point_gids:
-                    sc_features, radius = self._clip_to_sc_point(
+                    sc_features, sc_radius = self._clip_to_sc_point(
                         sc_point, tie_line_voltage, nn_sinks=nn_sinks,
                         clipping_buffer=clipping_buffer, radius=radius)
                     if sc_features.empty:
@@ -582,7 +598,7 @@ def _process_multi_core(self, capacity_class, tie_line_voltage,
                                         self._cost_fpath,
                                         sc_point.copy(deep=True),
                                         sc_features, capacity_class,
-                                        radius=radius,
+                                        radius=sc_radius,
                                         xmission_config=self._config,
                                         barrier_mult=barrier_mult,
                                         min_line_length=self._min_line_len,
@@ -638,7 +654,8 @@ def _process_single_core(self, capacity_class, tie_line_voltage,
             Flag to return least cost paths as a multi-line geometry,
             by default False
         radius : None | int, optional
-            Force clipping radius if set to an int
+            Force clipping radius if set to an int. Radius will be
+            expanded to include at least one connection feature.
         simplify_geo : float | None, optional
             If float, simplify geometries using this value
 
@@ -653,7 +670,7 @@ def _process_single_core(self, capacity_class, tie_line_voltage,
         for i, (_, sc_point) in enumerate(self.sc_points.iterrows(), start=1):
             gid = sc_point['sc_point_gid']
             if gid in sc_point_gids:
-                sc_features, radius = self._clip_to_sc_point(
+                sc_features, sc_radius = self._clip_to_sc_point(
                     sc_point, tie_line_voltage, nn_sinks=nn_sinks,
                     clipping_buffer=clipping_buffer, radius=radius)
                 if sc_features.empty:
@@ -663,7 +680,7 @@ def _process_single_core(self, capacity_class, tie_line_voltage,
                     self._cost_fpath,
                     sc_point.copy(deep=True),
                     sc_features, capacity_class,
-                    radius=radius,
+                    radius=sc_radius,
                     xmission_config=self._config,
                     barrier_mult=barrier_mult,
                     min_line_length=self._min_line_len,
@@ -722,7 +739,8 @@ def run(cls, cost_fpath, features_fpath, capacity_class, resolution=128,
             Flag to return least costs path as a multi-line geometry,
             by default False
         radius : None | int, optional
-            Force clipping radius if set to an int
+            Force clipping radius if set to an int. Radius will be
+            expanded to include at least one connection feature.
         simplify_geo : float | None, optional
             If float, simplify geometries using this value
 
@@ -764,7 +782,8 @@ class ReinforcedXmission(LeastCostXmission):
     """
 
     def __init__(self, cost_fpath, features_fpath, balancing_areas_fpath,
-                 resolution=128, xmission_config=None, min_line_length=0):
+                 resolution=128, xmission_config=None, min_line_length=0,
+                 allow_connections_within_states=False):
         """
         Parameters
         ----------
@@ -787,6 +806,10 @@ def __init__(self, cost_fpath, features_fpath, balancing_areas_fpath,
             By default, ``None``.
         min_line_length : int | float, optional
             Minimum line length in km. By default, ``0``.
+        allow_connections_within_states : bool, optional
+            Allow supply curve points to connect to substations outside
+            of their own BA, as long as all connections stay within the
+            same state. By default, ``False``.
         """
         super().__init__(cost_fpath=cost_fpath,
                          features_fpath=features_fpath,
@@ -795,6 +818,7 @@ def __init__(self, cost_fpath, features_fpath, balancing_areas_fpath,
                          min_line_length=min_line_length)
         self._ba = (gpd.read_file(balancing_areas_fpath)
                     .to_crs(self.features.crs))
+        self.allow_connections_within_states = allow_connections_within_states
 
     @staticmethod
     def _load_trans_feats(features_fpath):
@@ -818,8 +842,23 @@ def _clip_to_sc_point(self, sc_point, tie_line_voltage, nn_sinks=2,
 
         point = self.sc_points.loc[sc_point.name:sc_point.name].centroid
         ba_str = point.apply(ba_mapper(self._ba)).values[0]
-        mask = self.features["ba_str"] == ba_str
+        if self.allow_connections_within_states:
+            state = self._ba[self._ba["ba_str"] == ba_str]["state"].values[0]
+            logger.debug('  - Clipping features to {!r}'.format(state))
+            state_nodes = self._ba[self._ba["state"] == state]
+            allowed_bas = set(state_nodes["ba_str"])
+        else:
+            allowed_bas = {ba_str}
+        logger.debug("  - Clipping features to allowed ba's: {}"
+                     .format(allowed_bas))
+        mask = self.features["ba_str"].isin(allowed_bas)
         sc_features = self.features.loc[mask].copy(deep=True)
+        logger.debug('{} transmission features found in clipped area '
+                     .format(len(sc_features)))
+
+        if radius is not None:
+            sc_features = self._clip_to_radius(sc_point, radius, sc_features,
+                                               clipping_buffer)
 
         mask = self.features['max_volts'] >= tie_line_voltage
         sc_features = sc_features.loc[mask].copy(deep=True)
@@ -839,8 +878,9 @@ def _clip_to_sc_point(self, sc_point, tie_line_voltage, nn_sinks=2,
     def run(cls, cost_fpath, features_fpath, balancing_areas_fpath,
             capacity_class, resolution=128, xmission_config=None,
             min_line_length=0, sc_point_gids=None, clipping_buffer=1.05,
-            barrier_mult=100, max_workers=None, save_paths=False,
-            simplify_geo=None):
+            barrier_mult=100, max_workers=None, simplify_geo=None,
+            allow_connections_within_states=False, save_paths=False,
+            radius=None):
         """
         Find Least Cost Transmission connections between desired
         sc_points and substations in their balancing area.
@@ -880,11 +920,21 @@ def run(cls, cost_fpath, features_fpath, balancing_areas_fpath,
         max_workers : int, optional
             Number of workers to use for processing. If 1 run in serial,
             if ``None`` use all available cores. By default, ``None``.
+        simplify_geo : float | None, optional
+            If float, simplify geometries using this value.
+        allow_connections_within_states : bool, optional
+            Allow supply curve points to connect to substations outside
+            of their own BA, as long as all connections stay within the
+            same state. By default, ``False``.
         save_paths : bool, optional
             Flag to save reinforcement line path as a multi-line
             geometry. By default, ``False``.
-        simplify_geo : float | None, optional
-            If float, simplify geometries using this value.
+        radius : None | int, optional
+            Force clipping radius. Substations beyond this radius will
+            not be considered for connection with supply curve point.
+            Radius will be expanded to include at least one connection
+            feature. This value must be given in units of pixels
+            corresponding to the cost raster.
 
         Returns
         -------
@@ -895,14 +945,15 @@ def run(cls, cost_fpath, features_fpath, balancing_areas_fpath,
         """
         ts = time.time()
         lcx = cls(cost_fpath, features_fpath, balancing_areas_fpath,
-                  resolution=resolution, xmission_config=xmission_config,
-                  min_line_length=min_line_length)
+                  resolution, xmission_config, min_line_length,
+                  allow_connections_within_states)
         least_costs = lcx.process_sc_points(capacity_class,
                                             sc_point_gids=sc_point_gids,
                                             clipping_buffer=clipping_buffer,
                                             barrier_mult=barrier_mult,
                                             max_workers=max_workers,
                                             save_paths=save_paths,
+                                            radius=radius,
                                             simplify_geo=simplify_geo)
 
         logger.info('{} connections were made to {} SC points in {:.4f} '
diff --git a/reVX/least_cost_xmission/least_cost_xmission_cli.py b/reVX/least_cost_xmission/least_cost_xmission_cli.py
index 2789adecd..bc9f9540b 100644
--- a/reVX/least_cost_xmission/least_cost_xmission_cli.py
+++ b/reVX/least_cost_xmission/least_cost_xmission_cli.py
@@ -26,6 +26,7 @@
                                                           ReinforcedXmission)
 from reVX.least_cost_xmission.config import (TRANS_LINE_CAT, LOAD_CENTER_CAT,
                                              SINK_CAT, SUBSTATION_CAT)
+from reVX.least_cost_xmission.least_cost_paths import min_reinforcement_costs
 
 TRANS_CAT_TYPES = [TRANS_LINE_CAT, LOAD_CENTER_CAT, SINK_CAT, SUBSTATION_CAT]
 
@@ -144,6 +145,9 @@ def from_config(ctx, config, verbose):
               show_default=True, default=100,
               help=("Tranmission barrier multiplier, used when computing the "
                     "least cost tie-line path"))
+@click.option('--state_connections', '-acws', is_flag=True,
+              help='Flag to allow substations ot connect to any endpoints '
+                   'within their state. Default is not verbose.')
 @click.option('--max_workers', '-mw', type=INT,
               show_default=True, default=None,
               help=("Number of workers to use for processing, if 1 run in "
@@ -170,8 +174,8 @@ def from_config(ctx, config, verbose):
 def local(ctx, cost_fpath, features_fpath, balancing_areas_fpath,
           capacity_class, resolution, xmission_config, min_line_length,
           sc_point_start_index, sc_point_step_index, nn_sinks,
-          clipping_buffer, barrier_mult, max_workers, out_dir, log_dir,
-          verbose, save_paths, radius, simplify_geo):
+          clipping_buffer, barrier_mult, state_connections, max_workers,
+          out_dir, log_dir, verbose, save_paths, radius, simplify_geo):
     """
     Run Least Cost Xmission on local hardware
     """
@@ -188,33 +192,25 @@ def local(ctx, cost_fpath, features_fpath, balancing_areas_fpath,
     sce = SupplyCurveExtent(cost_fpath, resolution=resolution)
     sc_point_gids = list(sce.points.index.values)
     sc_point_gids = sc_point_gids[sc_point_start_index::sc_point_step_index]
+    kwargs = {"resolution": resolution,
+              "xmission_config": xmission_config,
+              "min_line_length": min_line_length,
+              "sc_point_gids": sc_point_gids,
+              "clipping_buffer": clipping_buffer,
+              "barrier_mult": barrier_mult,
+              "max_workers": max_workers,
+              "save_paths": save_paths,
+              "simplify_geo": simplify_geo,
+              "radius": radius}
     if balancing_areas_fpath is not None:
+        kwargs["allow_connections_within_states"] = state_connections
         least_costs = ReinforcedXmission.run(cost_fpath, features_fpath,
                                              balancing_areas_fpath,
-                                             capacity_class,
-                                             resolution=resolution,
-                                             xmission_config=xmission_config,
-                                             min_line_length=min_line_length,
-                                             sc_point_gids=sc_point_gids,
-                                             clipping_buffer=clipping_buffer,
-                                             barrier_mult=barrier_mult,
-                                             max_workers=max_workers,
-                                             save_paths=save_paths,
-                                             simplify_geo=simplify_geo)
+                                             capacity_class, **kwargs)
     else:
+        kwargs["nn_sinks"] = nn_sinks
         least_costs = LeastCostXmission.run(cost_fpath, features_fpath,
-                                            capacity_class,
-                                            resolution=resolution,
-                                            xmission_config=xmission_config,
-                                            min_line_length=min_line_length,
-                                            sc_point_gids=sc_point_gids,
-                                            nn_sinks=nn_sinks,
-                                            clipping_buffer=clipping_buffer,
-                                            barrier_mult=barrier_mult,
-                                            max_workers=max_workers,
-                                            save_paths=save_paths,
-                                            radius=radius,
-                                            simplify_geo=simplify_geo)
+                                            capacity_class, **kwargs)
     if len(least_costs) == 0:
         logger.error('No paths found.')
         return
@@ -294,6 +290,9 @@ def merge_output(ctx, split_to_geojson, out_file, out_dir, drop,  # noqa
         logger.info('Simplifying geometries by {}'.format(simplify_geo))
         df.geometry = df.geometry.simplify(simplify_geo)
 
+    if all(col in df for col in ["gid", "reinforcement_cost_per_mw"]):
+        df = min_reinforcement_costs(df)
+
     if not split_to_geojson:
         out_file = ('combo_{}'.format(files[0])
                     if out_file is None else out_file)
@@ -349,7 +348,8 @@ def merge_reinforcement_costs(ctx, cost_fpath, reinforcement_cost_fpath,
 
     logger.info("Merging reinforcement costs into transmission costs...")
 
-    r_cols = ["reinforcement_dist_km", "reinforcement_cost_per_mw"]
+    r_cols = ["ba_str", "reinforcement_poi_lat", "reinforcement_poi_lon",
+              "reinforcement_dist_km", "reinforcement_cost_per_mw"]
     costs[r_cols] = r_costs.loc[costs["trans_gid"], r_cols].values
 
     logger.info("Writing output to {!r}".format(out_file))
@@ -394,6 +394,8 @@ def get_node_cmd(config, start_index=0):
             '-log {}'.format(SLURM.s(config.log_directory)),
             ]
 
+    if config.allow_connections_within_states:
+        args.append('-acws')
     if config.save_paths:
         args.append('--save_paths')
     if config.radius:
@@ -436,6 +438,7 @@ def run_local(ctx, config):
                nn_sinks=config.nn_sinks,
                clipping_buffer=config.clipping_buffer,
                barrier_mult=config.barrier_mult,
+               state_connections=config.allow_connections_within_states,
                max_workers=config.execution_control.max_workers,
                out_dir=config.dirout,
                log_dir=config.log_directory,
diff --git a/tests/test_xmission_least_cost.py b/tests/test_xmission_least_cost.py
index c4a6e4ea2..a1f463e52 100644
--- a/tests/test_xmission_least_cost.py
+++ b/tests/test_xmission_least_cost.py
@@ -88,7 +88,8 @@ def ri_ba():
     ba_str, shapes = zip(*[("p{}".format(int(v)), shape(p))
                            for p, v in s if int(v) != 0])
 
-    return gpd.GeoDataFrame({"ba_str": ba_str}, crs=profile['crs'],
+    return gpd.GeoDataFrame({"ba_str": ba_str, "state": "Rhode Island"},
+                            crs=profile['crs'],
                             geometry=list(shapes))
 
 
@@ -183,7 +184,8 @@ def test_resolution(resolution):
     check_baseline(truth, test)
 
 
-def test_cli(runner):
+@pytest.mark.parametrize("save_paths", [False, True])
+def test_cli(runner, save_paths):
     """
     Test CostCreator CLI
     """
@@ -201,7 +203,8 @@ def test_cli(runner):
             "cost_fpath": COST_H5,
             "features_fpath": FEATURES,
             "capacity_class": f'{capacity}MW',
-            "min_line_length": 5.76
+            "min_line_length": 5.76,
+            "save_paths": save_paths,
         }
         config_path = os.path.join(td, 'config.json')
         with open(config_path, 'w') as f:
@@ -213,16 +216,24 @@ def test_cli(runner):
                .format(traceback.print_exception(*result.exc_info)))
         assert result.exit_code == 0, msg
 
-        test = '{}_{}MW_128.csv'.format(os.path.basename(td), capacity)
-        test = os.path.join(td, test)
-        test = pd.read_csv(test)
+        if save_paths:
+            test = '{}_{}MW_128.gpkg'.format(os.path.basename(td), capacity)
+            test = os.path.join(td, test)
+            test = gpd.read_file(test)
+            assert test.geometry is not None
+        else:
+            test = '{}_{}MW_128.csv'.format(os.path.basename(td), capacity)
+            test = os.path.join(td, test)
+            test = pd.read_csv(test)
         SupplyCurve._check_substation_conns(test, sc_cols='sc_point_gid')
         check_baseline(truth, test)
 
     LOGGERS.clear()
 
 
-def test_reinforcement_cli(runner, ri_ba):
+@pytest.mark.parametrize("save_paths", [False, True])
+@pytest.mark.parametrize("state_connections", [False, True])
+def test_reinforcement_cli(runner, ri_ba, save_paths, state_connections):
     """
     Test Reinforcement cost routines and CLI
     """
@@ -254,8 +265,11 @@ def test_reinforcement_cli(runner, ri_ba):
             "balancing_areas_fpath": ri_ba_path,
             "capacity_class": f'{capacity}MW',
             "barrier_mult": 100,
-            "min_line_length": 0
+            "min_line_length": 0,
+            "save_paths": save_paths,
+            "allow_connections_within_states": state_connections,
         }
+
         config_path = os.path.join(td, 'config.json')
         with open(config_path, 'w') as f:
             json.dump(config, f)
@@ -266,16 +280,23 @@ def test_reinforcement_cli(runner, ri_ba):
                .format(traceback.print_exception(*result.exc_info)))
         assert result.exit_code == 0, msg
 
-        test = '{}_{}MW_128.csv'.format(os.path.basename(td), capacity)
-        test = os.path.join(td, test)
-        test = pd.read_csv(test)
-
-        assert len(test) == 13
+        if save_paths:
+            test = '{}_{}MW_128.gpkg'.format(os.path.basename(td), capacity)
+            test = os.path.join(td, test)
+            test = gpd.read_file(test)
+            assert test.geometry is not None
+        else:
+            test = '{}_{}MW_128.csv'.format(os.path.basename(td), capacity)
+            test = os.path.join(td, test)
+            test = pd.read_csv(test)
+
+        assert len(test) == 71 if state_connections else 13
         assert set(test.trans_gid.unique()) == {69130}
         assert set(test.ba_str.unique()) == {"p4"}
 
         assert "poi_lat" in test
         assert "poi_lon" in test
+        assert "ba_str" in test
 
         assert len(test.poi_lat.unique()) == 1
         assert len(test.poi_lon.unique()) == 1
diff --git a/tests/test_xmission_least_cost_paths.py b/tests/test_xmission_least_cost_paths.py
index 87f568f78..bac44e844 100644
--- a/tests/test_xmission_least_cost_paths.py
+++ b/tests/test_xmission_least_cost_paths.py
@@ -59,7 +59,9 @@ def ba_regions_and_network_nodes():
     ba_str, shapes = zip(*[("p{}".format(int(v)), shape(p))
                            for p, v in s if int(v) != 0])
 
-    ri_ba = gpd.GeoDataFrame({"ba_str": ba_str}, crs=profile['crs'],
+    state = ["Rhode Island"] * len(ba_str)
+    ri_ba = gpd.GeoDataFrame({"ba_str": ba_str, "state": state},
+                             crs=profile['crs'],
                              geometry=list(shapes))
 
     ri_network_nodes = ri_ba.copy()
@@ -111,7 +113,8 @@ def test_parallel(max_workers):
     check(truth, test)
 
 
-def test_cli(runner):
+@pytest.mark.parametrize("save_paths", [False, True])
+def test_cli(runner, save_paths):
     """
     Test CostCreator CLI
     """
@@ -129,6 +132,7 @@ def test_cli(runner):
             "cost_fpath": COST_H5,
             "features_fpath": FEATURES,
             "capacity_class": f'{capacity}MW',
+            "save_paths": save_paths,
         }
         config_path = os.path.join(td, 'config.json')
         with open(config_path, 'w') as f:
@@ -145,15 +149,24 @@ def test_cli(runner):
         capacity_class = xmission_config._parse_cap_class(capacity)
         cap = xmission_config['power_classes'][capacity_class]
         kv = xmission_config.capacity_to_kv(capacity_class)
-        test = '{}_{}MW_{}kV.csv'.format(os.path.basename(td), cap, kv)
-        test = os.path.join(td, test)
-        test = pd.read_csv(test)
+        if save_paths:
+            test = '{}_{}MW_{}kV.gpkg'.format(os.path.basename(td), cap, kv)
+            test = os.path.join(td, test)
+            test = gpd.read_file(test)
+            assert test.geometry is not None
+        else:
+            test = '{}_{}MW_{}kV.csv'.format(os.path.basename(td), cap, kv)
+            test = os.path.join(td, test)
+            test = pd.read_csv(test)
         check(truth, test)
 
     LOGGERS.clear()
 
 
-def test_reinforcement_cli(runner, ba_regions_and_network_nodes):
+@pytest.mark.parametrize("save_paths", [False, True])
+@pytest.mark.parametrize("state_conns", [False, True])
+def test_reinforcement_cli(runner, ba_regions_and_network_nodes, save_paths,
+                           state_conns):
     """
     Test Reinforcement cost routines and CLI
     """
@@ -200,6 +213,8 @@ def test_reinforcement_cli(runner, ba_regions_and_network_nodes):
             "transmission_lines_fpath": ALLCONNS_FEATURES,
             "capacity_class": f"{capacity}MW",
             "barrier_mult": 100,
+            "save_paths": save_paths,
+            "allow_connections_within_states": state_conns
         }
         config_path = os.path.join(td, 'config.json')
         with open(config_path, 'w') as f:
@@ -215,24 +230,37 @@ def test_reinforcement_cli(runner, ba_regions_and_network_nodes):
         capacity_class = xmission_config._parse_cap_class(capacity)
         cap = xmission_config['power_classes'][capacity_class]
         kv = xmission_config.capacity_to_kv(capacity_class)
-        test = '{}_{}MW_{}kV.csv'.format(os.path.basename(td), cap, kv)
-        test = os.path.join(td, test)
-        test = pd.read_csv(test)
+        if save_paths:
+            test = '{}_{}MW_{}kV.gpkg'.format(os.path.basename(td), cap, kv)
+            test = os.path.join(td, test)
+            test = gpd.read_file(test)
+            assert test.geometry is not None
+        else:
+            test = '{}_{}MW_{}kV.csv'.format(os.path.basename(td), cap, kv)
+            test = os.path.join(td, test)
+            test = pd.read_csv(test)
 
         assert "reinforcement_poi_lat" in test
         assert "reinforcement_poi_lon" in test
         assert "poi_lat" not in test
         assert "poi_lon" not in test
-        assert len(test["reinforcement_poi_lat"].unique()) == 4
-        assert len(test["reinforcement_poi_lon"].unique()) == 4
+        assert "ba_str" in test
 
         assert len(test) == 69
         assert np.isclose(test.reinforcement_cost_per_mw.min(), 3332.695,
                           atol=0.001)
-        assert np.isclose(test.reinforcement_cost_per_mw.max(), 569757.740,
-                          atol=0.001)
         assert np.isclose(test.reinforcement_dist_km.min(), 1.918, atol=0.001)
         assert np.isclose(test.reinforcement_dist_km.max(), 80.353, atol=0.001)
+        if state_conns:
+            assert len(test["reinforcement_poi_lat"].unique()) == 3
+            assert len(test["reinforcement_poi_lon"].unique()) == 3
+            assert np.isclose(test.reinforcement_cost_per_mw.max(), 225129.798,
+                              atol=0.001)
+        else:
+            assert len(test["reinforcement_poi_lat"].unique()) == 4
+            assert len(test["reinforcement_poi_lon"].unique()) == 4
+            assert np.isclose(test.reinforcement_cost_per_mw.max(), 569757.740,
+                              atol=0.001)
 
     LOGGERS.clear()