diff --git a/documentation/apidoc/wntr.stormwater.gis.rst b/documentation/apidoc/wntr.stormwater.gis.rst
new file mode 100644
index 000000000..4aa561e02
--- /dev/null
+++ b/documentation/apidoc/wntr.stormwater.gis.rst
@@ -0,0 +1,8 @@
+wntr.stormwater.gis module
+=====================================
+
+.. automodule:: wntr.stormwater.gis
+    :members:
+    :inherited-members:
+    :no-undoc-members:
+    :show-inheritance:
diff --git a/documentation/apidoc/wntr.stormwater.graphics.rst b/documentation/apidoc/wntr.stormwater.graphics.rst
new file mode 100644
index 000000000..1bc7a0eec
--- /dev/null
+++ b/documentation/apidoc/wntr.stormwater.graphics.rst
@@ -0,0 +1,8 @@
+wntr.stormwater.graphics module
+=====================================
+
+.. automodule:: wntr.stormwater.graphics
+    :members:
+    :inherited-members:
+    :no-undoc-members:
+    :show-inheritance:
diff --git a/documentation/apidoc/wntr.stormwater.io.rst b/documentation/apidoc/wntr.stormwater.io.rst
index fa0ac4d22..d41687e41 100644
--- a/documentation/apidoc/wntr.stormwater.io.rst
+++ b/documentation/apidoc/wntr.stormwater.io.rst
@@ -3,5 +3,6 @@ wntr.stormwater.io module
 
 .. automodule:: wntr.stormwater.io
     :members:
+    :inherited-members:
     :no-undoc-members:
     :show-inheritance:
diff --git a/documentation/apidoc/wntr.stormwater.metrics.rst b/documentation/apidoc/wntr.stormwater.metrics.rst
new file mode 100644
index 000000000..6e122af12
--- /dev/null
+++ b/documentation/apidoc/wntr.stormwater.metrics.rst
@@ -0,0 +1,8 @@
+wntr.stormwater.metrics module
+=====================================
+
+.. automodule:: wntr.stormwater.metrics
+    :members:
+    :inherited-members:
+    :no-undoc-members:
+    :show-inheritance:
diff --git a/documentation/apidoc/wntr.stormwater.network.rst b/documentation/apidoc/wntr.stormwater.network.rst
index 4e9f615ff..e286fb891 100644
--- a/documentation/apidoc/wntr.stormwater.network.rst
+++ b/documentation/apidoc/wntr.stormwater.network.rst
@@ -3,5 +3,6 @@ wntr.stormwater.network module
 
 .. automodule:: wntr.stormwater.network
     :members:
+    :inherited-members:
     :no-undoc-members:
     :show-inheritance:
diff --git a/documentation/apidoc/wntr.stormwater.scenario.rst b/documentation/apidoc/wntr.stormwater.scenario.rst
new file mode 100644
index 000000000..c51b0bb98
--- /dev/null
+++ b/documentation/apidoc/wntr.stormwater.scenario.rst
@@ -0,0 +1,8 @@
+wntr.stormwater.scenario module
+=====================================
+
+.. automodule:: wntr.stormwater.scenario
+    :members:
+    :inherited-members:
+    :no-undoc-members:
+    :show-inheritance:
diff --git a/documentation/apidoc/wntr.stormwater.sim.rst b/documentation/apidoc/wntr.stormwater.sim.rst
index a10fc24ed..bb4eae6d6 100644
--- a/documentation/apidoc/wntr.stormwater.sim.rst
+++ b/documentation/apidoc/wntr.stormwater.sim.rst
@@ -3,5 +3,6 @@ wntr.stormwater.sim module
 
 .. automodule:: wntr.stormwater.sim
     :members:
+    :inherited-members:
     :no-undoc-members:
     :show-inheritance:
diff --git a/wntr/stormwater/__init__.py b/wntr/stormwater/__init__.py
index 120890d13..324f9c5da 100644
--- a/wntr/stormwater/__init__.py
+++ b/wntr/stormwater/__init__.py
@@ -2,10 +2,10 @@
 The wntr.stormwater package contains storm water resilience functionality 
 through integration of WNTR with SWMM, pyswmm, and swmmio
 """
-from wntr.stormwater.network import *
-from wntr.stormwater.sim import *
 from wntr.stormwater.gis import *
-from wntr.stormwater.fragility import *
-from wntr.stormwater.metrics import *
 from wntr.stormwater.graphics import *
-from wntr.stormwater.io import *
\ No newline at end of file
+from wntr.stormwater.io import *
+from wntr.stormwater.metrics import *
+from wntr.stormwater.network import *
+from wntr.stormwater.sim import *
+from wntr.stormwater.scenario import *
diff --git a/wntr/stormwater/gis.py b/wntr/stormwater/gis.py
index 2b579902a..0c8f27f03 100644
--- a/wntr/stormwater/gis.py
+++ b/wntr/stormwater/gis.py
@@ -1 +1,223 @@
 from wntr.gis import snap, intersect
+
+import swmmio
+
+import pandas as pd
+import numpy as np
+
+try:
+    from shapely.geometry import LineString, Point
+    has_shapely = True
+except ModuleNotFoundError:
+    has_shapely = False
+
+try:
+    import geopandas as gpd
+    has_geopandas = True
+except ModuleNotFoundError:
+    has_geopandas = False
+
+import wntr.network.elements
+
+class StormWaterNetworkGIS:
+    """
+    Storm Water network GIS class 
+    
+    Contains methods to create GeoDataFrames from StormWaterNetworkModel.
+    The ability to a create StormWaterNetworkModel from GeoDataFrames is 
+    not implemented.
+    
+    Parameters
+    ----------
+    gis_data : dict, optional
+        Dictionary of GeoDataFrames containing data to populate an instance 
+        of StormWaterNetworkGIS.  Valid dictionary keys are 
+        'junctions', 'outfalls','storage', 
+        'conduits', 'weirs', 'orifices', and 'pumps'
+    
+    Raises
+    ------
+    ModuleNotFoundError
+        if missing either shapely or geopandas
+    """
+    
+    def __init__(self, gis_data=None) -> None:
+        
+        if not has_shapely or not has_geopandas:
+            raise ModuleNotFoundError('shapley and geopandas are required')
+        
+        self.junctions = gpd.GeoDataFrame()
+        self.outfalls = gpd.GeoDataFrame()
+        self.storage = gpd.GeoDataFrame()
+        self.conduits = gpd.GeoDataFrame()
+        self.weirs = gpd.GeoDataFrame()
+        self.orifices = gpd.GeoDataFrame()
+        self.pumps = gpd.GeoDataFrame()
+        
+        if isinstance(gis_data, dict):
+            for name in ["junctions", "outfalls", "storage", 
+                         "conduits", "weirs", "orifices", "pumps"]:
+                gdf = getattr(self, name)
+                if name in gis_data.keys():
+                    assert isinstance(gis_data[name], gpd.GeoDataFrame)
+                    gdf = gis_data[name]
+
+    def _create_gis(self, swn, crs: str = None) -> None:
+        """
+        Create GIS data from a water network model.
+        
+        This method is used by wntr.network.io.to_gis
+        
+        Note: patterns, curves, rules, controls, sources, and options are not 
+        saved to the GIS data
+
+        Parameters
+        ----------
+        wn : WaterNetworkModel
+            Water network model
+        crs : str, optional
+            Coordinate reference system, by default None
+        """
+  
+        # create an updated swmmio model from inp file
+        filename = 'temp.inp'
+        swn._swmmio_model.inp.save(filename)
+        m = swmmio.Model(filename)
+        
+        self.junctions = m.nodes.geodataframe.loc[swn.junction_name_list,:]
+        self.outfalls = m.nodes.geodataframe.loc[swn.outfall_name_list,:]
+        self.storage = m.nodes.geodataframe.loc[swn.storage_name_list,:]
+        
+        self.conduits = m.links.geodataframe.loc[swn.conduit_name_list,:]
+        self.weirs = m.links.geodataframe.loc[swn.weir_name_list,:]
+        self.orifices = m.links.geodataframe.loc[swn.orifice_name_list,:]
+        self.pumps = m.links.geodataframe.loc[swn.pump_name_list,:]
+        
+        self.set_crs(crs, allow_override=True)
+        
+    def _create_swn(self, append=None):
+        raise NotImplementedError
+    
+    def to_crs(self, crs):
+        """
+        Transform CRS of the junctions, tanks, reservoirs, pipes, pumps,
+        and valves GeoDataFrames.
+
+        Calls geopandas.GeoDataFrame.to_crs on each GeoDataFrame.
+
+        Parameters
+        ----------
+        crs : str
+            Coordinate reference system
+        """
+        for name in ["junctions", "outfalls", "storage", 
+                     "conduits", "weirs", "orifices", "pumps"]:
+            gdf = getattr(self, name)
+            if 'geometry' in gdf.columns:
+                gdf = gdf.to_crs(crs, inplace=True)
+
+    def set_crs(self, crs, allow_override=False):
+        """
+        Set CRS of the junctions, tanks, reservoirs, pipes, pumps,
+        and valves GeoDataFrames.
+
+        Calls geopandas.GeoDataFrame.set_crs on each GeoDataFrame.
+
+        Parameters
+        ----------
+        crs : str
+            Coordinate reference system
+        allow_override : bool (optional)
+            Allow override of existing coordinate reference system
+        """
+        
+        for name in ["junctions", "outfalls", "storage", 
+                     "conduits", "weirs", "orifices", "pumps"]:
+            gdf = getattr(self, name)
+            if 'geometry' in gdf.columns:
+                gdf = gdf.set_crs(crs, inplace=True,
+                                   allow_override=allow_override)
+
+    def add_node_attributes(self, values, name):
+        raise NotImplementedError
+
+    def add_link_attributes(self, values, name):
+        raise NotImplementedError
+    
+    def _read(self, files, index_col='index'):
+        
+        for name in ["junctions", "outfalls", "storage", 
+                     "conduits", "weirs", "orifices", "pumps"]:
+            gdf = getattr(self, name)
+            if name in files.keys():
+                data = gpd.read_file(files[name]).set_index(index_col)
+                gdf = pd.concat([gdf, data])
+
+    def read_geojson(self, files, index_col='index'):
+        """
+        Append information from GeoJSON files to a WaterNetworkGIS object
+
+        Parameters
+        ----------
+        files : dictionary
+            Dictionary of GeoJSON filenames, where the keys are in the set 
+            ('junction', 'tanks', 'reservoirs', 'pipes', 'pumps', 'valves') and 
+            values are the corresponding GeoJSON filename
+        index_col : str, optional
+            Column that contains the element name
+        """
+        self._read(files, index_col)
+
+    def read_shapefile(self, files, index_col='index'):
+        raise NotImplementedError
+
+    def _write(self, prefix: str, driver="GeoJSON") -> None:
+        """
+        Write the WaterNetworkGIS object to GIS files
+
+        One file will be created for each type of network element (junctions, 
+        pipes, etc.) if those elements exists in the network
+        
+        Parameters
+        ----------
+        prefix : str
+            Filename prefix, will have the element type (junctions, 
+            pipes, etc.) appended
+        driver : str, optional
+            GeoPandas driver. Use "GeoJSON" for GeoJSON files, use :code:`None` 
+            for Esri Shapefile folders, by default "GeoJSON"
+
+        """
+        
+        if driver is None or driver == "":
+            extension = ""
+        else:
+            extension = "." + driver.lower()
+        
+        for name in ["junctions", "outfalls", "storage", 
+                     "conduits", "weirs", "orifices", "pumps"]:
+            gdf = getattr(self, name)
+            if len(gdf) > 0:
+                filename = prefix + "_" + name + extension
+                gdf.to_file(filename, driver=driver)
+ 
+    def write_geojson(self, prefix: str):
+        """
+        Write the WaterNetworkGIS object to a set of GeoJSON files, one file
+        for each network element.
+
+        Parameters
+        ----------
+        prefix : str
+            File prefix
+        """
+        self._write(prefix=prefix, driver="GeoJSON")
+
+    def write_shapefile(self, prefix: str):
+        raise NotImplementedError
+
+    def _valid_names(self, complete_list=True, truncate_names=None):
+        raise NotImplementedError
+    
+    def _shapefile_field_name_map(self):
+        raise NotImplementedError
diff --git a/wntr/stormwater/graphics.py b/wntr/stormwater/graphics.py
index c2e56e3c3..6d89c6f8c 100644
--- a/wntr/stormwater/graphics.py
+++ b/wntr/stormwater/graphics.py
@@ -2,187 +2,28 @@
 import networkx as nx
 import matplotlib.pylab as plt
 
-from wntr.graphics.network import plot_network as _plot_network
+from wntr.graphics.network import plot_network
+from wntr.graphics.curve import plot_fragility_curve
+from wntr.graphics.color import custom_colormap, random_colormap
 
 logger = logging.getLogger(__name__)
 
-"""
-def plot_network(swn, node_attribute=None, link_attribute=None, title=None,
-               node_size=20, node_range=[None,None], node_alpha=1, node_cmap=None, node_labels=False,
-               link_width=1, link_range=[None,None], link_alpha=1, link_cmap=None, link_labels=False,
-               add_colorbar=True, node_colorbar_label='Node', link_colorbar_label='Link', 
-               directed=False, ax=None, filename=None):
-    
-    if ax is None:  # create a new figure
-        plt.figure(facecolor="w", edgecolor="k")
-        ax = plt.gca()
-        
-    for subcatch in swn.subcatchments['geometry']:
-        ax.plot(*subcatch.boundary.xy, c='gray', linewidth=0.5)
-    
-    _plot_network(swn, node_attribute, link_attribute, title, 
-                  node_size, node_range, node_alpha, node_cmap, node_labels,
-                  link_width, link_range, link_alpha, link_cmap, link_labels,
-                  add_colorbar, node_colorbar_label, link_colorbar_label, 
-                  directed, ax=ax, filename=filename)
-"""
-def plot_network(
-    swn,
-    node_attribute=None,
-    link_attribute=None,
-    title=None,
-    node_size=20,
-    node_range=[None, None],
-    node_alpha=1,
-    node_cmap=None,
-    node_labels=False,
-    link_width=1,
-    link_range=[None, None],
-    link_alpha=1,
-    link_cmap=None,
-    link_labels=False,
-    add_colorbar=True,
-    node_colorbar_label="Node",
-    link_colorbar_label="Link",
-    directed=False,
-    ax=None,
-    filename=None,
-    inpdata=None,
-):
-
-    if ax is None:  # create a new figure
-        plt.figure(facecolor="w", edgecolor="k")
-        ax = plt.gca()
-
-    # Graph
-    G = swn.to_graph()
-    if not directed:
-        G = G.to_undirected()
-
-    # Position
-    pos = nx.get_node_attributes(G, "coords")
-    pos = dict((k, (v[0][0], v[0][1])) for k,v in pos.items())
-    if len(pos) == 0:
-        pos = None
-
-    # Define node properties
-    add_node_colorbar = add_colorbar
-    if node_attribute is not None:
-
-        if isinstance(node_attribute, list):
-            if node_cmap is None:
-                node_cmap = ["red", "red"]
-            add_node_colorbar = False
-
-        if node_cmap is None:
-            node_cmap = plt.get_cmap("Spectral_r")
-        elif isinstance(node_cmap, list):
-            if len(node_cmap) == 1:
-                node_cmap = node_cmap * 2
-            node_cmap = custom_colormap(len(node_cmap), node_cmap)
-
-        node_attribute = dict(node_attribute)
-        nodelist, nodecolor = zip(*node_attribute.items())
-
-    else:
-        nodelist = None
-        nodecolor = "k"
-
-    add_link_colorbar = add_colorbar
-    if link_attribute is not None:
 
-        if isinstance(link_attribute, list):
-            if link_cmap is None:
-                link_cmap = ["red", "red"]
-            add_link_colorbar = False
-
-        if link_cmap is None:
-            link_cmap = plt.get_cmap("Spectral_r")
-        elif isinstance(link_cmap, list):
-            if len(link_cmap) == 1:
-                link_cmap = link_cmap * 2
-            link_cmap = custom_colormap(len(link_cmap), link_cmap)
-
-        link_attribute = dict(link_attribute)
-
-        # Replace link_attribute dictionary defined as
-        # {link_name: attr} with {(start_node, end_node, link_name): attr}
-        attr = {}
-        edge_tuples = list(G.edges(keys=True))
-        for edge_tuple in list(G.edges(keys=True)):
-            edge_name = edge_tuple[2]
-            try:
-                value = link_attribute[edge_name]
-                attr[edge_tuple] = value
-            except:
-                pass
-        link_attribute = attr
-
-        linklist, linkcolor = zip(*link_attribute.items())
-    else:
-        linklist = None
-        linkcolor = "k"
-
-    if title is not None:
-        ax.set_title(title)
+# def plot_network(swn, node_attribute=None, link_attribute=None, title=None,
+#                node_size=20, node_range=[None,None], node_alpha=1, node_cmap=None, node_labels=False,
+#                link_width=1, link_range=[None,None], link_alpha=1, link_cmap=None, link_labels=False,
+#                add_colorbar=True, node_colorbar_label='Node', link_colorbar_label='Link', 
+#                directed=False, ax=None, filename=None):
     
-    for subcatch in swn.subcatchments['geometry']:
-        ax.plot(*subcatch.boundary.xy, c='gray', linewidth=0.5)
-
-    edge_background = nx.draw_networkx_edges(G, pos, edge_color="grey", width=0.5, ax=ax)
-
-    nodes = nx.draw_networkx_nodes(
-        G,
-        pos,
-        nodelist=nodelist,
-        node_color=nodecolor,
-        node_size=node_size,
-        alpha=node_alpha,
-        cmap=node_cmap,
-        vmin=node_range[0],
-        vmax=node_range[1],
-        linewidths=0,
-        ax=ax,
-    )
-    edges = nx.draw_networkx_edges(
-        G,
-        pos,
-        edgelist=linklist,
-        edge_color=linkcolor,
-        width=link_width,
-        alpha=link_alpha,
-        edge_cmap=link_cmap,
-        edge_vmin=link_range[0],
-        edge_vmax=link_range[1],
-        ax=ax,
-    )
-    if node_labels:
-        labels = dict(zip(swn.node_name_list, swn.node_name_list))
-        nx.draw_networkx_labels(G, pos, labels, font_size=7, ax=ax)
-    # if link_labels:
-    #     labels = {}
-    #     for link_name in swn.link_name_list:
-    #         link = wn.get_link(link_name)
-    #         labels[(link.start_node_name, link.end_node_name)] = link_name
-    #     nx.draw_networkx_edge_labels(G, pos, labels, font_size=7, ax=ax)
+#     if ax is None:  # create a new figure
+#         plt.figure(facecolor="w", edgecolor="k")
+#         ax = plt.gca()
         
-    if add_node_colorbar and node_attribute:
-        clb = plt.colorbar(nodes, shrink=0.5, pad=0, ax=ax)
-        clb.ax.set_title(node_colorbar_label, fontsize=10)
-    if add_link_colorbar and link_attribute:
-        if directed:
-            vmin = min(map(abs, link_attribute.values()))
-            vmax = max(map(abs, link_attribute.values()))
-            sm = plt.cm.ScalarMappable(cmap=link_cmap, norm=plt.Normalize(vmin=vmin, vmax=vmax))
-            sm.set_array([])
-            clb = plt.colorbar(sm, shrink=0.5, pad=0.05, ax=ax)
-        else:
-            clb = plt.colorbar(edges, shrink=0.5, pad=0.05, ax=ax)
-        clb.ax.set_title(link_colorbar_label, fontsize=10)
-
-    ax.axis("off")
-
-    if filename:
-        plt.savefig(filename)
-
-    return ax
+#     #for subcatch in swn.subcatchments['geometry']:
+#     #    ax.plot(*subcatch.boundary.xy, c='gray', linewidth=0.5)
+    
+#     _plot_network(swn, node_attribute, link_attribute, title, 
+#                   node_size, node_range, node_alpha, node_cmap, node_labels,
+#                   link_width, link_range, link_alpha, link_cmap, link_labels,
+#                   add_colorbar, node_colorbar_label, link_colorbar_label, 
+#                   directed, ax=ax, filename=filename)
diff --git a/wntr/stormwater/io.py b/wntr/stormwater/io.py
index f5982b4b8..e5d8f7971 100644
--- a/wntr/stormwater/io.py
+++ b/wntr/stormwater/io.py
@@ -1,20 +1,36 @@
 import logging
 import pandas as pd
+import networkx as nx
 import pyswmm
 import swmmio
 from swmm.toolkit.shared_enum import NodeAttribute, LinkAttribute, \
                                      SubcatchAttribute, SystemAttribute
 
 from wntr.sim import SimulationResults
-from wntr.gis import WaterNetworkGIS
+from wntr.stormwater.gis import StormWaterNetworkGIS
+import wntr.stormwater
 
 logger = logging.getLogger(__name__)
 
 
 def to_graph(swn):
-
+    """
+    Convert a StormWaterNetworkModel into a NetworkX graph
+    
+    Parameters
+    ----------
+    swn : StormWaterNetworkModel
+        Storm water network model
+        
+    Returns
+    -------
+    NetworkX graph
+    
+    """
     G = swn._swmmio_model.network
-
+    geom = nx.get_node_attributes(G, 'geometry')
+    pos = dict([(k,v['coordinates']) for k,v in geom.items()])
+    nx.set_node_attributes(G, pos, 'pos')
     return G
 
 
@@ -31,42 +47,58 @@ def to_gis(swn, crs=None):
 
     Returns
     -------
-    WaterNetworkGIS object that contains GeoDataFrames
+    StormWaterNetworkGIS object that contains GeoDataFrames
         
     """
-    gis_data = WaterNetworkGIS()
-    
-    # nodes
-    gis_data.junctions = swn.nodes.loc[swn.junction_name_list,:]
-    gis_data.outfalls = swn.nodes.loc[swn.outfall_name_list,:]
-    gis_data.storages = swn.nodes.loc[swn.storage_name_list,:]
-    
-    # links
-    gis_data.conduits = swn.links.loc[swn.conduit_name_list,:]
-    gis_data.weirs = swn.links.loc[swn.weir_name_list,:]
-    gis_data.orifices = swn.links.loc[swn.orifice_name_list,:]
-    gis_data.pumps = swn.links.loc[swn.pump_name_list,:]
-    
-    # subcatchments
-    gis_data.subcatchments = swn.subcatchments.loc[swn.subcatchment_name_list,:]
-
+    gis_data = StormWaterNetworkGIS()
+    gis_data._create_gis(swn, crs)
     return gis_data
 
 
 def write_inpfile(swn, filename):
-
+    """
+    Write the StormWaterNetworkModel to an EPANET INP file
+    
+    Parameters
+    ----------
+    swn : WaterNetworkModel
+        Water network model
+    filename : string
+       Name of the inp file
+    """
     swn._swmmio_model.inp.save(filename)
 
-
 def read_inpfile(filename):
+    """
+    Create a StormWaterNetworkModel from an SWMM INP file
+    
+    Parameters
+    ----------
+    filename : string
+       Name of the inp file
+       
+    Returns
+    -------
+    StormWaterNetworkModel
+    """
+    swn = wntr.stormwater.network.StormWaterNetworkModel(filename)
 
-    model = swmmio.Model(filename)
-
-    return model
+    return swn
 
 
-def read_outfile(outfile):
+def read_outfile(filename):
+    """
+    Read a SWMM binary output file
 
+    Parameters
+    ----------
+    filename : string
+       Name of the SWMM binary output file
+    
+    Returns
+    -------
+    SimulationResults
+    """
     results = SimulationResults()
     
     # Node results = INVERT_DEPTH, HYDRAULIC_HEAD, PONDED_VOLUME, 
@@ -88,7 +120,7 @@ def read_outfile(outfile):
     # EVAP_RATE
     results.system = {}
     
-    with pyswmm.Output(outfile) as out:
+    with pyswmm.Output(filename) as out:
         times = out.times
         
         for attribute in NodeAttribute:
@@ -117,4 +149,22 @@ def read_outfile(outfile):
             temp[attribute] = ts.values()
         results.system = pd.DataFrame(data=temp, index=times)
         
-    return results
\ No newline at end of file
+    return results
+
+def write_geojson(swn, prefix: str, crs=None):
+    """
+    Write the StormWaterNetworkModel to a set of GeoJSON files, one file for each
+    network element.
+
+    Parameters
+    ----------
+    swn : wntr StormWaterNetworkModel
+        Storm water network model
+    prefix : str
+        File prefix
+    crs : str, optional
+        Coordinate reference system, by default None
+    """
+
+    swn_gis = swn.to_gis(crs)
+    swn_gis.write_geojson(prefix=prefix)
diff --git a/wntr/stormwater/network.py b/wntr/stormwater/network.py
index f175ac38f..cd8cba66e 100644
--- a/wntr/stormwater/network.py
+++ b/wntr/stormwater/network.py
@@ -1,8 +1,10 @@
 import logging
+import random
+import numpy as np
 import pandas as pd
 import swmmio
 
-from wntr.stormwater.io import to_graph, to_gis, write_inpfile
+from wntr.stormwater.io import to_graph, to_gis
 
 logger = logging.getLogger(__name__)
 
@@ -13,7 +15,7 @@ class StormWaterNetworkModel(object):
 
     Unlike the WaterNetworkModel, this class has no iterator methods,
     add/remove methods, and no component registries.
-    
+
     Parameters
     -------------------
     inp_file_name: string 
@@ -25,37 +27,37 @@ def __init__(self, inp_file_name=None):
         
         if inp_file_name:
             self._swmmio_model = swmmio.Model(inp_file_name)
-            
+  
             # Attributes of StormWaterNetworkModel link to attributes 
             # in swmmio.Model.inp, which contains dataframes from an INP file.
             # The swmmio.Model.inp object also includes a .save method to 
             # write a new INP file.
-            
+ 
             # Nodes = Junctions, outfall, and storage nodes
             self.junctions = self._swmmio_model.inp.junctions
             self.outfalls = self._swmmio_model.inp.outfalls
             self.storage = self._swmmio_model.inp.junctions
-            
+
             # Links = Conduits, weirs, orifices, and pumps
             self.conduits = self._swmmio_model.inp.conduits
             self.weirs = self._swmmio_model.inp.weirs
             self.orifices = self._swmmio_model.inp.orifices
             self.pumps = self._swmmio_model.inp.pumps
-            
+
             self.subcatchments = self._swmmio_model.inp.subcatchments
             self.subareas = self._swmmio_model.inp.subareas
-            
-            self.raingages = self._swmmio_model.inp.raingages 
+
+            self.raingages = self._swmmio_model.inp.raingages
             self.infiltration = self._swmmio_model.inp.infiltration
             self.inflows = self._swmmio_model.inp.inflows
             #self.dwf = self._swmmio_model.inp.dwf
-            
+
             self.curves = self._swmmio_model.inp.curves
             self.timeseries = self._swmmio_model.inp.timeseries
-            
+
             self.options = self._swmmio_model.inp.options
             self.files = self._swmmio_model.inp.files
-            
+
             self.coordinates = self._swmmio_model.inp.coordinates
             self.vertices = self._swmmio_model.inp.vertices
             self.polygons = self._swmmio_model.inp.polygons
@@ -63,212 +65,362 @@ def __init__(self, inp_file_name=None):
 
         else:
             self._swmmio_model = None
-
+            
+    @property
+    # def junctions(self):
+    #     """Generator to get all junctions
+        
+    #     Yields
+    #     ------
+    #     name : str
+    #         The name of the junction
+    #     node : Junction
+    #         The junction object    
+            
+    #     """
+    #     for name, junc in self.junctions.iterrows():
+    #         yield name, junc
+            
     @property
     def node_name_list(self):
         """Get a list of node names
-        
+
         Returns
         -------
         list of strings
-        
+
         """
         return self.junction_name_list + self.outfall_name_list + \
             self.storage_name_list
-    
+
     @property
     def junction_name_list(self):
         """Get a list of junction names
-        
+
         Returns
         -------
         list of strings
-        
+
         """
         return list(self.junctions.index)
 
-
     @property
     def outfall_name_list(self):
         """Get a list of outfall names
-        
+
         Returns
         -------
         list of strings
-        
+
         """
         return list(self.outfalls.index)
-    
-    
+
     @property
     def storage_name_list(self):
         """Get a list of storage names
-        
+
         Returns
         -------
         list of strings
-        
+
         """
         return list(self.storage.index)
-    
+
     @property
     def link_name_list(self):
         """Get a list of link names
-        
+
         Returns
         -------
         list of strings
-        
+
         """
         return self.conduit_name_list + self.weir_name_list + \
             self.orifice_name_list + self.pump_name_list
-    
+
     @property
     def conduit_name_list(self):
         """Get a list of conduit names
-        
+
         Returns
         -------
         list of strings
-        
+
         """
         return list(self.conduits.index)
-    
+
     @property
     def weir_name_list(self):
         """Get a list of weir names
-        
+
         Returns
         -------
         list of strings
-        
+
         """
         return list(self.weirs.index)
-    
+
     @property
     def orifice_name_list(self):
         """Get a list of orifice names
-        
+
         Returns
         -------
         list of strings
-        
+
         """
         return list(self.orifices.index)
-    
+
     @property
     def pump_name_list(self):
         """Get a list of pump names
-        
+
         Returns
         -------
         list of strings
-        
+
         """
         return list(self.pumps.index)
-    
+
     @property
     def subcatchment_name_list(self):
         """Get a list of subcatchment names
-        
+
         Returns
         -------
         list of strings
-        
+
         """
         return list(self.subcatchments.index)
-    
+
     @property
     def raingage_name_list(self):
         """Get a list of raingage names
-        
+
         Returns
         -------
         list of strings
-        
+
         """
         return list(self.raingages.index)
-    
+
     @property
     def num_nodes(self):
         """The number of nodes"""
         return len(self.node_name_list)
-    
+
     @property
     def num_junctions(self):
         """The number of junctions"""
         return len(self.junction_name_list)
-    
+
     @property
     def num_outfalls(self):
         """The number of outfalls"""
         return len(self.outfall_name_list)
-    
+
     @property
     def num_storages(self):
         """The number of storages"""
         return len(self.storage_name_list)
-    
+
     @property
     def num_links(self):
         """The number of links"""
         return len(self.link_name_list)
-    
+
     @property
     def num_conduits(self):
         """The number of conduits"""
         return len(self.conduit_name_list)
-    
+
     @property
     def num_weirs(self):
         """The number of weirs"""
         return len(self.weir_name_list)
-    
+
     @property
     def num_orifices(self):
         """The number of orifices"""
         return len(self.orifice_name_list)
-    
+
     @property
     def num_pumps(self):
         """The number of pumps"""
         return len(self.pump_name_list)
-    
+
     @property
     def num_subcatchments(self):
         """The number of subcatchments"""
         return len(self.subcatchment_name_list)
-    
+
     @property
     def num_raingages(self):
         """The number of raingages"""
         return len(self.raingage_name_list)
-    
-    
-    def to_gis(self, crs=None):
-        
-        return to_gis(self)
-    
-    
-    def to_graph(self, node_weight=None, link_weight=None, 
-                 modify_direction=False):
+
+    def get_node(self, name):
+        """Get a specific node
+
+        Parameters
+        ----------
+        name : str
+            The node name
+
+        Returns
+        -------
+        Junction, Outfall, or Storage
+
         """
-        Convert a StormWaterNetworkModel into a networkx MultiDiGraph
-        
+        return Node(self, name)
+
+    def get_link(self, name):
+        """Get a specific link
+
         Parameters
         ----------
-        node_weight :  dict or pandas Series (optional)
-            Node weights
-        link_weight : dict or pandas Series (optional)
-            Link weights.  
-        modify_direction : bool (optional)
-            If True, than if the link weight is negative, the link start and 
-            end node are switched and the abs(weight) is assigned to the link
-            (this is useful when weighting graphs by flowrate). If False, link 
-            direction and weight are not changed.
+        name : str
+            The link name
+
+        Returns
+        -------
+        Pipe, Pump, or Valve
+
+        """
+        if name in self.conduit_name_list:
+            data = self.conduits
+        elif name in self.weir_name_list:
+            data = self.weirs
+        elif name in self.orifice_name_list:
+            data = self.orifices
+        elif name in self.pump_name_list:
+            data = self.pumps
             
+        start_node_name = data.loc[name, 'InletNode']
+        end_node_name = data.loc[name, 'OutletNode']
+        
+        return Link(self, name, start_node_name, end_node_name)
+
+    def to_gis(self, crs=None):
+        """
+        Convert a StormWaterNetworkModel into GeoDataFrames
+
+        Parameters
+        ----------
+        crs : str, optional
+            Coordinate reference system, by default None
+        """
+        return to_gis(self, crs)
+
+    def to_graph(self):
+        """
+        Convert a StormWaterNetworkModel into a networkx MultiDiGraph
+
         Returns
         --------
         networkx MultiDiGraph
         """
         return to_graph(self)
+
+
+class Node(object):
+    """
+    Base class for nodes.
+    """
+
+    def __init__(self, swn, name):
+
+        # Set the node name
+        self._name = name
+
+    @property
+    def name(self):
+        """str: The name of the node (read only)"""
+        return self._name
+
+
+class Link(object):
+    """
+    Base class for links.
+    """
+
+    def __init__(self, swn, link_name, start_node_name, end_node_name):
+
+        self._link_name = link_name
+        self._start_node_name = start_node_name
+        self._end_node_name = end_node_name
+
+    @property
+    def name(self):
+        """str: The link name (read-only)"""
+        return self._link_name
+
+    @property
+    def start_node_name(self):
+        """str: The name of the start node (read only)"""
+        return self._start_node_name
+
+    @property
+    def end_node_name(self):
+        """str: The name of the end node (read only)"""
+        return self._end_node_name
+
+def generate_valve_layer(swn, placement_type='strategic', n=1, seed=None):
+    """
+    Generate valve layer data, which can be used in valve segmentation analysis.
+
+    Parameters
+    -----------
+    swn : wntr StrormWaterNetworkModel
+        A StormWaterNetworkModel object
+        
+    placement_type : string
+        Options include 'strategic' and 'random'.  
+        
+        - If 'strategic', n is the number of pipes from each node that do not 
+          contain a valve. In this case, n is generally 0, 1 or 2 
+          (i.e. N, N-1, N-2 valve placement).
+        - If 'random', then n randomly placed valves are used to define the 
+          valve layer.
+        
+    n : int
+        
+        - If 'strategic', n is the number of pipes from each node that do not 
+          contain a valve.
+        - If 'random', n is the number of number of randomly placed valves.
+        
+    seed : int or None
+        Random seed
+       
+    Returns
+    ---------
+    valve_layer : pandas DataFrame
+        Valve layer, defined by node and link pairs (for example, valve 0 is 
+        on link A and protects node B). The valve_layer DataFrame is indexed by
+        valve number, with columns named 'node' and 'link'.
+    """
+    
+    if seed is not None:
+        np.random.seed(seed)
+        random.seed(seed)
+    valve_layer = []
+    if placement_type=='random':
+        s = swn.conduits['InletNode']
+        all_valves = list(tuple(zip(s.index,s)))
+        s = swn.conduits['OutletNode']
+        all_valves.extend(tuple(zip(s.index,s)))
+        for valve_tuple in random.sample(all_valves, n):
+            pipe_name, node_name = valve_tuple
+            valve_layer.append([pipe_name, node_name])
+            
+    elif placement_type == 'strategic':
+        s = pd.concat([swn.conduits['InletNode'], swn.conduits['OutletNode']])
+        s = s.to_frame('Node')
+        for node_name, group in s.groupby('Node'):
+            links = list(group.index)
+            for l in np.random.choice(links, max(len(links)-n,0), replace=False):
+                valve_layer.append([l, node_name])
+            
+    valve_layer = pd.DataFrame(valve_layer, columns=['link', 'node'])  
+    
+    return valve_layer
diff --git a/wntr/stormwater/fragility.py b/wntr/stormwater/scenario.py
similarity index 100%
rename from wntr/stormwater/fragility.py
rename to wntr/stormwater/scenario.py
diff --git a/wntr/stormwater/sim.py b/wntr/stormwater/sim.py
index 502a5a1e4..885bf7aa0 100644
--- a/wntr/stormwater/sim.py
+++ b/wntr/stormwater/sim.py
@@ -8,12 +8,14 @@
 
 
 class SWMMSimulator(object):
-
+    """
+    SWMM simulator class.
+    """
     def __init__(self, swn):
         self._swn = swn
 
     def run_sim(self, file_prefix='temp'):
-
+        """Run a SWMM simulation"""
         inpfile = file_prefix + '.inp'
         outfile = file_prefix + '.out'