Repository: USEPA/WNTR (language: python)

Issue:
Bug fix in EPANET INP reader and minor update for pandas
This PR includes the following changes:

* Bug fix in EPANET INP reader, allow DEMANDS section to override demands in JUNCTIONS section, this update closes #68
* Update use of timedelta, based on changes in pandas



Discussion / hints:

Fix patch:
diff --git a/documentation/acronyms.rst b/documentation/acronyms.rst
index b7fce00aa..c72b00c1f 100644
--- a/documentation/acronyms.rst
+++ b/documentation/acronyms.rst
@@ -18,6 +18,8 @@ Abbreviations
 
 **EPA**: Environmental Protection Agency
 
+**GIS**: Geographic information system
+ 
 **HDF**: Hierarchical Data Format
 
 **IDE**: Integrated development environment
diff --git a/documentation/apidoc/wntr.gis.geospatial.rst b/documentation/apidoc/wntr.gis.geospatial.rst
new file mode 100644
index 000000000..b1ecdef10
--- /dev/null
+++ b/documentation/apidoc/wntr.gis.geospatial.rst
@@ -0,0 +1,7 @@
+wntr.gis.geospatial module
+=====================================
+
+.. automodule:: wntr.gis.geospatial
+    :members:
+    :no-undoc-members:
+    :show-inheritance:
diff --git a/documentation/apidoc/wntr.gis.network.rst b/documentation/apidoc/wntr.gis.network.rst
new file mode 100644
index 000000000..236c31d19
--- /dev/null
+++ b/documentation/apidoc/wntr.gis.network.rst
@@ -0,0 +1,7 @@
+wntr.gis.network module
+=====================================
+
+.. automodule:: wntr.gis.network
+    :members:
+    :no-undoc-members:
+    :show-inheritance:
diff --git a/documentation/apidoc/wntr.gis.rst b/documentation/apidoc/wntr.gis.rst
new file mode 100644
index 000000000..69a81375c
--- /dev/null
+++ b/documentation/apidoc/wntr.gis.rst
@@ -0,0 +1,15 @@
+wntr.gis package
+====================
+
+.. automodule:: wntr.gis
+    :members:
+    :no-undoc-members:
+    :show-inheritance:
+
+Submodules
+----------
+
+.. toctree::
+
+   wntr.gis.network
+   wntr.gis.geospatial
diff --git a/documentation/apidoc/wntr.network.io.rst b/documentation/apidoc/wntr.network.io.rst
new file mode 100644
index 000000000..190a5e351
--- /dev/null
+++ b/documentation/apidoc/wntr.network.io.rst
@@ -0,0 +1,8 @@
+wntr.network.io module
+===================================
+
+.. automodule:: wntr.network.io
+    :members:
+    :inherited-members:
+    :no-undoc-members:
+    :show-inheritance:
diff --git a/documentation/apidoc/wntr.network.rst b/documentation/apidoc/wntr.network.rst
index 6bce88993..05dbe11ca 100644
--- a/documentation/apidoc/wntr.network.rst
+++ b/documentation/apidoc/wntr.network.rst
@@ -14,6 +14,7 @@ Submodules
    wntr.network.base
    wntr.network.controls
    wntr.network.elements
+   wntr.network.io
    wntr.network.model
    wntr.network.layer
    wntr.network.options
diff --git a/documentation/apidoc/wntr.rst b/documentation/apidoc/wntr.rst
index aa8906a59..abab02a50 100644
--- a/documentation/apidoc/wntr.rst
+++ b/documentation/apidoc/wntr.rst
@@ -16,6 +16,7 @@ Subpackages
     
     wntr.epanet
     wntr.graphics
+    wntr.gis
     wntr.metrics
     wntr.morph
     wntr.network
diff --git a/documentation/disaster_models.rst b/documentation/disaster_models.rst
index a93fd1df8..6b946ec59 100644
--- a/documentation/disaster_models.rst
+++ b/documentation/disaster_models.rst
@@ -150,6 +150,14 @@ The method :class:`~wntr.network.elements.Pump.add_outage` adds time controls to
 When simulating power outages, consider placing check bypasses around pumps 
 and check valves next to reservoirs.
 
+.. note:: 
+   The power outage is added to the WaterNetworkModel as a Rule.  
+   If the WaterNetworkModel includes Controls that impact pumps, 
+   it is recommended that those Controls are converted to Rules using 
+   to ensure that the power outage is evaluated in a consistent manner.  
+   The method :class:`~wntr.network.model.WaterNetworkModel.convert_controls_to_rules`
+   converts all Controls to Rules, with an option input argument to set priority.  
+   
 Fires
 ----------------
 WNTR can be used to simulate damage caused to system components due to fire and/or to simulate water usage due to fighting fires. To fight fires, additional water is drawn from the system. Fire codes vary by 
diff --git a/documentation/figures/intersect_demographics.png b/documentation/figures/intersect_demographics.png
new file mode 100644
index 000000000..68bb2fc83
Binary files /dev/null and b/documentation/figures/intersect_demographics.png differ
diff --git a/documentation/figures/intersect_earthquake.png b/documentation/figures/intersect_earthquake.png
new file mode 100644
index 000000000..19c3dc3e1
Binary files /dev/null and b/documentation/figures/intersect_earthquake.png differ
diff --git a/documentation/figures/intersect_landslide.png b/documentation/figures/intersect_landslide.png
new file mode 100644
index 000000000..fe2fbcabc
Binary files /dev/null and b/documentation/figures/intersect_landslide.png differ
diff --git a/documentation/figures/snap_lines.png b/documentation/figures/snap_lines.png
new file mode 100644
index 000000000..6684bfc99
Binary files /dev/null and b/documentation/figures/snap_lines.png differ
diff --git a/documentation/figures/snap_points.png b/documentation/figures/snap_points.png
new file mode 100644
index 000000000..526ec8f48
Binary files /dev/null and b/documentation/figures/snap_points.png differ
diff --git a/documentation/framework.rst b/documentation/framework.rst
index 850151d22..8d40b27bd 100644
--- a/documentation/framework.rst
+++ b/documentation/framework.rst
@@ -35,6 +35,7 @@ run simulations are described in more detail below, followed by a list of softwa
    :class:`~wntr.sim`		                           Contains classes and methods to run hydraulic and water quality simulations using the water network model.
    :class:`~wntr.metrics`	                           Contains functions to compute resilience, including hydraulic, water quality, water security, and economic metrics. Methods to compute topographic metrics are included in the wntr.network.graph module.
    :class:`~wntr.morph`	                              Contains methods to modify water network model morphology, including network skeletonization, modifying node coordinates, and splitting or breaking pipes.
+   :class:`~wntr.gis`	                              Contains geospatial capabilities, including a function to convert the water network model to GeoDataFrames.
    :class:`~wntr.graphics`                            Contains functions to generate graphics.
    :class:`~wntr.epanet`                              Contains EPANET 2.00.12 compatibility class and methods for WNTR.
    :class:`~wntr.utils`                               Contains helper functions.
diff --git a/documentation/gis.rst b/documentation/gis.rst
new file mode 100644
index 000000000..3c0fc08e3
--- /dev/null
+++ b/documentation/gis.rst
@@ -0,0 +1,652 @@
+.. raw:: latex
+
+    \clearpage
+
+.. doctest::
+    :hide:
+
+    >>> import wntr
+    >>> import numpy as np
+    >>> import pandas as pd
+    >>> pd.options.display.expand_frame_repr = False
+    >>> pd.options.display.precision = 3
+    >>> try:
+    ...    import geopandas as gpd
+    ... except ModuleNotFoundError:
+    ...    gpd = None
+    >>> import matplotlib.pylab as plt
+    >>> from os.path import isdir
+    >>> examples_dir = '../examples'
+    >>> if not isdir(examples_dir):
+    ...    examples_dir = 'examples'
+    >>> wn = wntr.network.model.WaterNetworkModel(examples_dir+'/networks/Net1.inp')
+    
+.. doctest::
+    :hide:
+    :skipif: gpd is None
+	
+    >>> demographic_data = gpd.read_file(examples_dir+'/data/Net1_demographic_data.geojson')
+    >>> landslide_data = gpd.read_file(examples_dir+'/data/Net1_landslide_data.geojson')
+    >>> earthquake_data = gpd.read_file(examples_dir+'/data/Net1_earthquake_data.geojson')
+    >>> hydrant_data = gpd.read_file(examples_dir+'/data/Net1_hydrant_data.geojson')
+    >>> valve_data = gpd.read_file(examples_dir+'/data/Net1_valve_data.geojson')
+
+.. _geospatial:
+
+Geospatial capabilities
+======================================
+
+The junctions, tanks, reservoirs, pipes, pumps, and valves in a water network model can be converted to 
+geospatial data objects. These objects can be used 
+directly within WNTR, 
+in Python with geospatial packages such as GeoPandas and Shapely, and 
+in geographic information system (GIS) platforms.
+Open source GIS platforms include QGIS and GRASS GIS.
+
+.. note:: 
+   The GIS capabilities described below require the Python packages **geopandas** [JVFM21]_ 
+   and **rtree** [rtree]_, both are optional dependencies of WNTR.
+   Note that **shapely** is installed with geopandas.
+	
+The following capabilities use data stored in GeoPandas GeoDataFrames.  
+Each GeoDataFrame contain a `geometry` column which contains 
+geometric objects commonly used in geospatial analysis.
+
+:numref:`table-geometry-type` includes water network model components and the 
+geometry type that defines each component.  
+Geometry types include 
+``shapely.geometry.Point``, ``shapely.geometry.LineString``, and ``shapely.geometry.MultiLineString``.
+A few components can be defined using multiple types:
+
+* Pumps and valves can be stored as Points or LineStrings. While Pumps are defined as 
+  lines within WNTR (and EPANET), converting the geometry to Points can be useful for geospatial analysis and visualization.
+* Pipes that do not contain vertices are stored as a LineString while pipes that contain 
+  vertices are stored as a MultiLineString.
+
+.. _table-geometry-type:
+.. table:: Geometry types for water network model components
+   
+   ==============================  ===============================
+   Water Network model component   Shapely geometry type
+   ==============================  ===============================
+   Junction                        Point
+   Tank                            Point
+   Reservoir                       Point
+   Pipe                            LineString or MultiLineString 
+   Pump                            LineString or Point
+   Valve                           LineString or Point
+   ==============================  ===============================
+
+External datasets that are defined using geospatial data objects can be utilized within 
+WNTR to add attributes to the water network model and the analysis.  Example external datasets include:
+
+* Point geometries that contain utility billing data, hydrant locations, isolation valve locations, or the location of emergency services.
+  These geometries can be associated with points and lines in a water network model by snapping the point to the nearest component.
+* LineString or MultiLineString geometries that contain street layout or earthquake fault lines.
+  These geometries can be associated with points and lines in a water network model by finding the intersection.
+* Polygon geometries that contain elevation, building footprints, zoning, land cover, hazard maps, census data, demographics, or social vulnerability.
+  These geometries can be associated with points and lines in a water network model by finding the intersection.
+
+Convert a water network model to GeoDataFrames
+----------------------------------------------
+
+A water network model can be converted to a group of GeoDataFrames using the function
+:class:`~wntr.model.network.io.to_gis`.
+One GeoDataFrame is created for each type of model component (Junction, Tank, Reservoir, Pipe, Pump, Valve).
+Component attributes (e.g., junction elevation) are stored in the GeoDataFrame along with the 
+component's geometry.
+
+.. note:: 
+   When converting a water network model into GeoDataFrames, the user can specify the 
+   coordinate reference system (CRS) for the node coordinates.
+   By default, the CRS is not specified (and is set to None).  
+   CRSs can be geographic (e.g., latitude/longitude where the units are in degrees) or 
+   projected (e.g., Universal Transverse Mercator where units are in meters).
+   Projected CRSs are preferred for more accurate distance calculations.  
+   GeoPandas includes documentation on managing projections and includes the ``to_crs`` method to convert between CRSs 
+   Additionally, WNTR includes methods to modify coordinates, see :ref:`modify_node_coords` for more information.
+   
+The following example creates GeoDataFrames from a water network model.  Note that Net1 is used in the following example.
+The CRS EPSG:4326 assumes that the coordinates are in degrees latitude and degrees longitude (which is not realistic for this network).  
+
+.. doctest::
+    :skipif: gpd is None
+
+    >>> import wntr # doctest: +SKIP
+	
+    >>> wn = wntr.network.WaterNetworkModel('networks/Net1.inp') # doctest: +SKIP
+    >>> wn_gis = wn.to_gis(crs='EPSG:4326')
+	
+Individual GeoDataFrames are obtained as follows (Note that Net1 has no valves and the GeoDataFrame for valves is empty).
+
+.. doctest::
+    :skipif: gpd is None
+
+    >>> wn_gis.junctions # doctest: +SKIP
+    >>> wn_gis.tanks # doctest: +SKIP
+    >>> wn_gis.reservoirs # doctest: +SKIP
+    >>> wn_gis.pipes # doctest: +SKIP
+    >>> wn_gis.pumps # doctest: +SKIP
+    >>> wn_gis.valves # doctest: +SKIP
+	
+For example, the junctions GeoDataFrame contains the following information
+
+.. doctest::
+    :skipif: gpd is None
+
+    >>> print(wn_gis.junctions.head()) # doctest: +SKIP
+              type  elevation  ... base_demand                   geometry
+    name                       ...                                       
+    10    Junction    216.408  ...    0.000000  POINT (20.00000 70.00000)
+    11    Junction    216.408  ...    0.009464  POINT (30.00000 70.00000)
+    12    Junction    213.360  ...    0.009464  POINT (50.00000 70.00000)
+    13    Junction    211.836  ...    0.006309  POINT (70.00000 70.00000)
+    21    Junction    213.360  ...    0.009464  POINT (30.00000 40.00000)
+
+The GeoDataFrames can be saved to GeoJSON files using the :class:`~wntr.gis.network.WaterNetworkGIS.write` method.
+
+.. doctest::
+    :skipif: gpd is None
+
+    >>> wn_gis.write('Net1', driver='GeoJSON')
+	
+This creates the following GeoJSON files for junctions, tanks, reservoirs, pipes, pumps, and valves:
+
+* Net1_junctions.geojson
+* Net1_tanks.geojson
+* Net1_reservoirs.geojson
+* Net1_pipes.geojson
+* Net1_pumps.geojson
+* Note, Net1_valves.geojson not created since Net1 has no valves
+
+The `write` method includes an optional argument that sets the GeoPandas driver.  
+This allows the user to create Shapefiles in addition to GeoJSON files.
+These files can be loaded into GIS platforms for further analysis and visualization.
+
+The water network model also contains methods to write GeoJSON and Shapefiles directly (without converting the water network model to GeoDataFrames).
+See :ref:`giswrite` for more information.
+
+Add attributes to the GeoDataFrames
+----------------------------------------
+
+Additional attributes can be added to the GeoDataFrames using the 
+:class:`~wntr.gis.network.WaterNetworkGIS.add_node_attributes` and 
+:class:`~wntr.gis.network.WaterNetworkGIS.add_link_attributes` methods.
+Additional attributes, such as simulation results or resilience metric, can be used in further analysis and visualization.
+
+The following example adds simulated pressure at hour 1 to nodes (which includes junctions, tanks, and reservoirs).
+
+.. doctest::
+    :skipif: gpd is None
+
+    >>> sim = wntr.sim.EpanetSimulator(wn)
+    >>> results = sim.run_sim()
+    >>> wn_gis.add_node_attributes(results.node['pressure'].loc[3600,:], 'Pressure_1hr')
+	
+Attributes can also be added directly to individual GeoDataFrames, as shown below.
+
+.. doctest::
+    :skipif: gpd is None
+
+    >>> wn_gis.junctions['new attribute'] = 10
+	
+
+.. _gisread:
+
+Convert GeoDataFrames to a water network model
+----------------------------------------------
+
+A water network model can be created from a group of GeoDataFrames using the function
+:class:`~wntr.network.io.from_gis`. The water network model can be created from 
+a :class:`~wntr.gis.network.WaterNetworkGIS` object or a dictionary of GeoDataFrames, with the following keys:
+junctions, tanks, reservoirs, pipes, pumps, valves.
+
+.. note:: 
+   A water network model created from GeoDataFrames only contains 
+   junction, tank, reservoir, pipe, pump and valve
+   attributes and topographic connectivity of the network.  
+   The network will **NOT** contain patterns, curves, rules, controls, 
+   or sources.  Water network model options are set to default values. 
+   Additional functionality could be added to WNTR in a future release.
+
+
+The following example creates a water network model from the group of GeoDataFrames created above.
+
+.. doctest::
+    :skipif: gpd is None
+
+    >>> wn2 = wntr.network.from_gis(wn_gis)
+
+Snap point geometries to the nearest point or line
+----------------------------------------------------
+
+The :class:`~wntr.gis.snap` function is used to find the nearest point or line to a set of points. 
+This functionality can be used to assign hydrants to junctions or assign isolation valves to pipes.
+
+When snapping Point geometries in GeoDataFrame A to Point or Line geometries in GeoDataFrame B, 
+the function returns the following information (one entry for each point in A):
+
+* Nearest point or line in B
+* Distance between original and snapped point
+* Coordinates of the snapped point
+* If B contains Lines, the nearest endpoint along the nearest line
+* If B contains Lines, the relative distance from the line's start node (line position)
+
+The following examples used geospatial data stored in the `examples/data <https://github.com/USEPA/WNTR/blob/main/examples/data>`_ directory.
+The GeoPandas ``read_file`` method is used to read the GeoJSON files into GeoDataFrames.  
+
+Snap hydrants to junctions
+^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+The dataset used in this example defines hydrant locations.  
+The dataset is a GeoDataFrame with a `geometry` column that contains ``shapely.geometry.Point`` geometries and a 
+`demand` column that defines fire flow requirements.
+
+.. doctest::
+    :skipif: gpd is None
+    
+    >>> import geopandas as gpd # doctest: +SKIP
+	
+    >>> hydrant_data = gpd.read_file('data/Net1_hydrant_data.geojson') # doctest: +SKIP
+    >>> print(hydrant_data)
+       demand                   geometry
+    0    5000  POINT (48.20000 37.20000)
+    1    1500  POINT (71.80000 68.30000)
+    2    8000  POINT (51.20000 71.10000)
+
+The following example uses the function :class:`~wntr.gis.snap` to snap hydrant locations to the nearest junction.
+	
+.. doctest::
+    :skipif: gpd is None
+
+    >>> snapped_to_junctions = wntr.gis.snap(hydrant_data, wn_gis.junctions, tolerance=5.0)
+    >>> print(snapped_to_junctions)
+      node  snap_distance                   geometry
+    0   22          3.329  POINT (50.00000 40.00000)
+    1   13          2.476  POINT (70.00000 70.00000)
+    2   12          1.628  POINT (50.00000 70.00000)
+
+The data, water network model, and snapped points can be plotted as follows.
+
+.. doctest::
+    :skipif: gpd is None
+
+    >>> ax = hydrant_data.plot()
+    >>> ax = wntr.graphics.plot_network(wn, node_attribute=snapped_to_junctions['node'].to_list(), ax=ax)
+
+.. doctest::
+    :skipif: gpd is None
+    :hide:
+    
+    >>> bounds = ax.axis('equal')
+    >>> plt.tight_layout()
+    >>> plt.savefig('snap_points.png', dpi=300)
+
+:numref:`fig-snap-points` illustrates hydrants snapped to junctions.
+
+.. _fig-snap-points:
+.. figure:: figures/snap_points.png
+   :width: 800
+   :alt: napped points to points
+
+   Example hydrants snapped to junctions
+
+** By reversing the order of GeoDataFrames in the snap function**,
+the nearest hydrant to each junction can also be identified.
+Note that the tolerance is increased to ensure all junctions are assigned a hydrant.
+   
+.. doctest::
+    :skipif: gpd is None
+
+    >>> snapped_to_hydrants = wntr.gis.snap(wn_gis.junctions, hydrant_data, tolerance=100.0)
+    >>> print(snapped_to_hydrants)
+        node  snap_distance                   geometry
+    10     2         31.219  POINT (51.20000 71.10000)
+    11     2         21.229  POINT (51.20000 71.10000)
+    12     2          1.628  POINT (51.20000 71.10000)
+    13     1          2.476  POINT (71.80000 68.30000)
+    21     0         18.414  POINT (48.20000 37.20000)
+    22     0          3.329  POINT (48.20000 37.20000)
+    23     0         21.979  POINT (48.20000 37.20000)
+    31     0         32.727  POINT (48.20000 37.20000)
+    32     0         27.259  POINT (48.20000 37.20000)
+
+Snap valves to pipes
+^^^^^^^^^^^^^^^^^^^^^
+
+The dataset used in this example defines valve locations.  
+The dataset is a GeoDataFrame with a `geometry` column that contains ``shapely.geometry.Point`` geometries.
+
+.. doctest::
+    :skipif: gpd is None
+
+    >>> valve_data = gpd.read_file('data/Net1_valve_data.geojson') # doctest: +SKIP
+    >>> print(valve_data)
+                        geometry
+    0  POINT (56.50000 41.50000)
+    1  POINT (32.10000 67.60000)
+    2  POINT (52.70000 86.30000)
+	
+The following example uses the function :class:`~wntr.gis.snap` to snap valve locations to the nearest pipe.
+
+.. doctest::
+    :skipif: gpd is None
+
+    >>> snapped_to_pipes = wntr.gis.snap(valve_data, wn_gis.pipes, tolerance=5.0)
+    >>> print(snapped_to_pipes)
+      link node  snap_distance  line_position                   geometry
+    0   22   22            1.5          0.325  POINT (56.50000 40.00000)
+    1  111   11            2.1          0.080  POINT (30.00000 67.60000)
+    2  110    2            2.7          0.185  POINT (50.00000 86.30000)
+	
+The snapped locations can be used to define a :ref:`valvelayer` and then create network segments.
+
+.. doctest::
+    :skipif: gpd is None
+
+    >>> valve_layer = snapped_to_pipes[['link', 'node']]
+    >>> G = wn.get_graph()
+    >>> node_segments, link_segments, segment_size = wntr.metrics.valve_segments(G, valve_layer)
+
+The data, water network model, and valve layer can be plotted as follows.
+
+.. doctest::
+    :skipif: gpd is None
+
+    >>> ax = valve_data.plot()
+    >>> ax = wntr.graphics.plot_valve_layer(wn, valve_layer, add_colorbar=False, ax=ax)
+
+.. doctest::
+    :skipif: gpd is None 
+    :hide:
+    
+    >>> bounds = ax.axis('equal')
+    >>> plt.tight_layout()
+    >>> plt.savefig('snap_lines.png', dpi=300)
+
+:numref:`fig-snap-lines` illustrates the valve layer created by snapping points to lines.
+
+.. _fig-snap-lines:
+.. figure:: figures/snap_lines.png
+   :width: 600
+   :alt: Snapped points to lines
+
+   Example valve layer created by snapping points to lines
+
+Find the intersect between geometries
+--------------------------------------
+
+The :class:`~wntr.gis.intersect`  function is used to find the intersection between geometries.
+This functionality can be used to identify faults or landslides that intersect pipes,
+or assign demographic data to network components.
+
+When finding the intersection of GeoDataFrame A with GeoDataFrame B (where A and B can contain Points, Lines, or Polygons),
+the function returns the following information (one entry for each geometry in A):
+
+* List of intersecting B geometry indices
+* Number of intersecting B geometries
+
+The following additional information is returned when geometries in B are assigned a value:
+
+* List of intersecting B geometry values
+* Minimum B geometry value
+* Maximum B geometry value
+* Mean B geometry value
+* If A contains Lines and B contains Polygons, weighted mean value (weighted by intersecting length)
+
+When the B geometry contains Polygons, the user can optionally include the background in the intersection.
+This is useful when working with geometries that do not cover the entire region of interest.
+For example, while census tracts cover the entire region, hazard maps might contain gaps (regions with no hazard) 
+that the user might want to include in the intersection.
+
+Assign earthquake probability to pipes
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+The dataset used in this example contains earthquake fault lines.
+The dataset is a GeoDataFrame with a `geometry` column that contains ``shapely.geometry.LineString`` geometries and a 
+`Pr` column which contains probability of an earthquake over magnitude 7.
+
+.. doctest::
+    :skipif: gpd is None
+
+    >>> earthquake_data = gpd.read_file('data/Net1_earthquake_data.geojson') # doctest: +SKIP
+    >>> print(earthquake_data)
+         Pr                                           geometry
+    0  0.50  LINESTRING (36.00000 2.00000, 44.00000 44.0000...
+    1  0.75  LINESTRING (42.00000 2.00000, 45.00000 27.0000...
+    2  0.90  LINESTRING (40.00000 2.00000, 50.00000 50.0000...
+    3  0.25  LINESTRING (30.00000 2.00000, 35.00000 30.0000...
+	
+The following example uses the function :class:`~wntr.gis.intersect` to assign earthquake probability to pipes. 
+
+.. doctest::
+    :skipif: gpd is None
+
+    >>> pipe_Pr = wntr.gis.intersect(wn_gis.pipes, earthquake_data, 'Pr')
+    >>> print(pipe_Pr)
+        intersections                  values  n   sum   min   max  mean
+    10             []                      []  0   NaN   NaN   NaN   NaN
+    11            [1]                  [0.75]  1  0.75  0.75  0.75  0.75
+    12      [3, 2, 0]        [0.25, 0.9, 0.5]  3  1.65  0.25  0.90  0.55
+    21   [1, 3, 2, 0]  [0.75, 0.25, 0.9, 0.5]  4  2.40  0.25  0.90  0.60
+    22             []                      []  0   NaN   NaN   NaN   NaN
+    31   [1, 3, 2, 0]  [0.75, 0.25, 0.9, 0.5]  4  2.40  0.25  0.90  0.60
+    110            []                      []  0   NaN   NaN   NaN   NaN
+    111            []                      []  0   NaN   NaN   NaN   NaN
+    112     [3, 2, 0]        [0.25, 0.9, 0.5]  3  1.65  0.25  0.90  0.55
+    113           [0]                   [0.5]  1  0.50  0.50  0.50  0.50
+    121            []                      []  0   NaN   NaN   NaN   NaN
+    122            []                      []  0   NaN   NaN   NaN   NaN
+	
+The data and water network model can be plotted as follows.  
+Pipes are colored with their max probability.
+
+.. doctest::
+    :skipif: gpd is None
+
+    >>> ax = earthquake_data.plot(column='Pr', alpha=0.5, cmap='bone', vmin=0, vmax=1)
+    >>> ax = wntr.graphics.plot_network(wn, link_attribute=pipe_Pr['max'], link_width=1.5, 
+    ...     node_range=[0,1], link_range=[0,1], ax=ax)
+
+.. doctest::
+    :skipif: gpd is None
+    :hide:
+    
+    >>> bounds = ax.axis('equal')
+    >>> plt.tight_layout()
+    >>> plt.savefig('intersect_earthquake.png', dpi=300)
+	
+:numref:`fig-intersect-earthquake` illustrates the intersection of junctions and pipes with the lines. 
+
+.. _fig-intersect-earthquake:
+.. figure:: figures/intersect_earthquake.png
+   :width: 800
+   :alt: Intersection of junctions and pipes with earthquake fault lines
+
+   Example intersection of junctions and pipes with earthquake fault lines
+
+The intersect function can also be used to identify pipes that cross each fault simply by reversing 
+the order in which the geometries intersect, as shown below:
+
+.. doctest::
+    :skipif: gpd is None
+
+    >>> pipes_that_intersect_each_fault = wntr.gis.intersect(earthquake_data, wn_gis.pipes)
+    >>> print(pipes_that_intersect_each_fault)
+                intersections  n
+    0  [21, 31, 112, 12, 113]  5
+    1            [21, 31, 11]  3
+    2       [21, 31, 112, 12]  4
+    3       [21, 31, 112, 12]  4
+
+Assign landslide probability to pipes
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+The dataset used in this example contains landslide hazard zones.
+The dataset is a GeoDataFrame with a `geometry` column that contains ``shapely.geometry.LineString`` geometries and a 
+`Pr` column which contains probability of damage from a landslide in that zone.
+
+.. doctest::
+    :skipif: gpd is None
+
+    >>> landslide_data = gpd.read_file('data/Net1_landslide_data.geojson') # doctest: +SKIP
+    >>> print(landslide_data)
+         Pr                                           geometry
+    0  0.50  POLYGON ((28.84615 22.23077, 28.76040 22.05079...
+    1  0.75  POLYGON ((40.00708 1.83192, 33.00708 84.83192,...
+    2  0.90  POLYGON ((58.05971 44.48507, 58.11776 44.67615...
+	
+The following example uses the function :class:`~wntr.gis.intersect` to assign earthquake hazard map with pipes. 
+This is very similar to the earthquake example above, except that the landslide hazards are polygons.  Additionally, since the 
+hazard map does not include a "background" value that defines the probability of damage outside landslide zones, 
+the background conditions are included in the intersection function.
+
+.. doctest::
+    :skipif: gpd is None 
+
+    >>> pipe_Pr = wntr.gis.intersect(wn_gis.pipes, landslide_data, 'Pr', include_background=True, 
+    ...    background_value=0)
+    >>> print(pipe_Pr)
+              intersections            values  n   sum  min   max   mean  weighted_mean
+    10         [BACKGROUND]             [0.0]  1  0.00  0.0  0.00  0.000          0.000
+    11      [BACKGROUND, 1]       [0.0, 0.75]  2  0.75  0.0  0.75  0.375          0.201
+    12         [BACKGROUND]             [0.0]  1  0.00  0.0  0.00  0.000          0.000
+    21   [BACKGROUND, 1, 0]  [0.0, 0.75, 0.5]  3  1.25  0.0  0.75  0.417          0.394
+    22      [BACKGROUND, 2]        [0.0, 0.9]  2  0.90  0.0  0.90  0.450          0.246
+    31      [BACKGROUND, 1]       [0.0, 0.75]  2  0.75  0.0  0.75  0.375          0.212
+    110        [BACKGROUND]             [0.0]  1  0.00  0.0  0.00  0.000          0.000
+    111     [BACKGROUND, 0]        [0.0, 0.5]  2  0.50  0.0  0.50  0.250          0.352
+    112        [BACKGROUND]             [0.0]  1  0.00  0.0  0.00  0.000          0.000
+    113        [BACKGROUND]             [0.0]  1  0.00  0.0  0.00  0.000          0.000
+    121     [BACKGROUND, 0]        [0.0, 0.5]  2  0.50  0.0  0.50  0.250          0.250
+    122        [BACKGROUND]             [0.0]  1  0.00  0.0  0.00  0.000          0.000
+
+The data and water network model can be plotted as follows.  
+Pipes are colored with the weighted mean probability.
+
+.. doctest::
+    :skipif: gpd is None
+
+    >>> ax = landslide_data.plot(column='Pr', alpha=0.5, cmap='bone', vmin=0, vmax=1)
+    >>> ax = wntr.graphics.plot_network(wn, link_attribute=pipe_Pr['weighted_mean'], link_width=1.5, 
+    ...     node_range=[0,1], link_range=[0,1], ax=ax)
+
+.. doctest::
+    :skipif: gpd is None
+    :hide:
+    
+    >>> bounds = ax.axis('equal')
+    >>> plt.tight_layout()
+    >>> plt.savefig('intersect_landslide.png', dpi=300)
+	
+:numref:`fig-intersect-landslide` illustrates the intersection. 
+
+.. _fig-intersect-landslide:
+.. figure:: figures/intersect_landslide.png
+   :width: 800
+   :alt: Intersection of junctions and pipes with landslide zones
+
+   Example intersection of pipes with landslide zones
+   
+**By reversing the order of GeoDataFrames in the intersection function**, the pipes that intersect each landslide zone and information about 
+the intersecting pipe diameters can also be identified:
+
+.. doctest::
+    :skipif: gpd is None 
+
+    >>> pipes_that_intersect_each_landslide = wntr.gis.intersect(landslide_data, wn_gis.pipes, 'diameter')
+    >>> print(pipes_that_intersect_each_landslide)
+        intersections                                             values  n    sum    min    max   mean
+    0  [111, 121, 21]                             [0.254, 0.2032, 0.254]  3  0.711  0.203  0.254  0.237
+    1    [21, 11, 31]  [0.254, 0.35559999999999997, 0.15239999999999998]  3  0.762  0.152  0.356  0.254
+    2            [22]                              [0.30479999999999996]  1  0.305  0.305  0.305  0.305
+	
+Assign demographic data to pipes and junctions
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+The dataset used in this example contains demographic data.
+The dataset is a GeoDataFrame with a `geometry` column that contains ``shapely.geometry.Polygon`` geometries along with 
+columns that store mean income, mean age, and population.
+
+.. doctest::
+    :skipif: gpd is None
+
+    >>> demographic_data = gpd.read_file('data/Net1_demographic_data.geojson') # doctest: +SKIP
+    >>> print(demographic_data)
+       mean_income  mean_age  population                                           geometry
+    0      63326.0      35.0      3362.0  POLYGON ((41.67813 82.75023, 41.98596 60.85779...
+    1      78245.0      31.0      5618.0  POLYGON ((23.21084 40.19160, 22.99063 27.71777...
+    2      91452.0      40.0      5650.0  POLYGON ((22.99063 27.71777, 61.93720 16.36165...
+    3      54040.0      39.0      5546.0  POLYGON ((61.93720 16.36165, 22.99063 27.71777...
+    4      26135.0      38.0      5968.0  POLYGON ((61.93720 16.36165, 64.04456 22.10119...
+    5      57620.0      31.0      4315.0  POLYGON ((44.48497 87.21487, 79.81144 71.92669...
+    6      44871.0      54.0      4547.0  POLYGON ((64.04456 22.10119, 51.72994 45.92347...
+    7      69067.0      55.0      2541.0  POLYGON ((46.01047 99.15725, 46.40654 99.33204...
+	
+The following example uses the function :class:`~wntr.gis.intersect`
+to assign demographic data to junctions and pipes.  
+
+.. doctest::
+    :skipif: gpd is None
+
+    >>> junction_demographics = wntr.gis.intersect(wn_gis.junctions, demographic_data, 'mean_income')
+    >>> print(junction_demographics)
+       intersections     values  n      sum      min      max     mean
+    10           [0]  [63326.0]  1  63326.0  63326.0  63326.0  63326.0
+    11           [0]  [63326.0]  1  63326.0  63326.0  63326.0  63326.0
+    12           [5]  [57620.0]  1  57620.0  57620.0  57620.0  57620.0
+    13           [5]  [57620.0]  1  57620.0  57620.0  57620.0  57620.0
+    21           [3]  [54040.0]  1  54040.0  54040.0  54040.0  54040.0
+    22           [3]  [54040.0]  1  54040.0  54040.0  54040.0  54040.0
+    23           [6]  [44871.0]  1  44871.0  44871.0  44871.0  44871.0
+    31           [2]  [91452.0]  1  91452.0  91452.0  91452.0  91452.0
+    32           [2]  [91452.0]  1  91452.0  91452.0  91452.0  91452.0
+	
+.. doctest::
+    :skipif: gpd is None
+
+    >>> pipe_demographics = wntr.gis.intersect(wn_gis.pipes, demographic_data, 'mean_income')
+    >>> print(pipe_demographics)
+        intersections              values  n       sum      min      max     mean  weighted_mean
+    10            [0]           [63326.0]  1   63326.0  63326.0  63326.0  63326.0      63326.000
+    11         [0, 5]  [63326.0, 57620.0]  2  120946.0  57620.0  63326.0  60473.0      61002.920
+    12            [5]           [57620.0]  1   57620.0  57620.0  57620.0  57620.0      57620.000
+    21            [3]           [54040.0]  1   54040.0  54040.0  54040.0  54040.0      54040.000
+    22         [3, 6]  [54040.0, 44871.0]  2   98911.0  44871.0  54040.0  49455.5      47067.895
+    31            [2]           [91452.0]  1   91452.0  91452.0  91452.0  91452.0      91452.000
+    110        [5, 7]  [57620.0, 69067.0]  2  126687.0  57620.0  69067.0  63343.5      60580.117
+    111        [0, 3]  [63326.0, 54040.0]  2  117366.0  54040.0  63326.0  58683.0      60953.558
+    112        [5, 3]  [57620.0, 54040.0]  2  111660.0  54040.0  57620.0  55830.0      56596.728
+    113        [5, 6]  [57620.0, 44871.0]  2  102491.0  44871.0  57620.0  51245.5      53707.370
+    121        [3, 2]  [54040.0, 91452.0]  2  145492.0  54040.0  91452.0  72746.0      73586.482
+    122        [3, 2]  [54040.0, 91452.0]  2  145492.0  54040.0  91452.0  72746.0      66314.037
+	
+The data and water network model can be plotted as follows. 
+Junctions and pipes are colored with their average value (weighted average for pipes).
+
+.. doctest::
+    :skipif: gpd is None
+
+    >>> ax = demographic_data.plot(column='mean_income', alpha=0.5,  
+    ...     cmap='bone', vmin=10000, vmax=100000)
+    >>> ax = wntr.graphics.plot_network(wn, node_attribute=junction_demographics['mean'], 
+    ...     link_attribute=pipe_demographics['weighted_mean'], link_width=1.5,
+    ...     node_range=[40000,80000], link_range=[40000,80000], ax=ax)
+						   
+.. doctest::
+    :skipif: gpd is None
+    :hide:
+    
+    >>> bounds = ax.axis('equal')
+    >>> plt.tight_layout()
+    >>> plt.savefig('intersect_demographics.png', dpi=300)
+
+:numref:`fig-intersect-demographics` illustrates the intersection of junctions and pipes with the polygons. Note that the color scale for 
+the polygons is different than the junction and pipe attributes.  
+
+.. _fig-intersect-demographics:
+.. figure:: figures/intersect_demographics.png
+   :width: 800
+   :alt: Intersection of junctions and pipes with demographic data
+
+   Example intersection of junctions and pipes with demographic data
\ No newline at end of file
diff --git a/documentation/index.rst b/documentation/index.rst
index c63d54128..b5c0280bf 100644
--- a/documentation/index.rst
+++ b/documentation/index.rst
@@ -27,6 +27,7 @@ designed to simulate and analyze resilience of water distribution networks.
    fragility
    morph
    graphics
+   gis
    advancedsim
    license
    whatsnew
diff --git a/documentation/installation.rst b/documentation/installation.rst
index f2243fc0c..da6056f81 100644
--- a/documentation/installation.rst
+++ b/documentation/installation.rst
@@ -218,15 +218,23 @@ The following Python packages are optional:
   http://python-visualization.github.io/folium/
 * utm [Bieni19]_: used to translate node coordinates to utm and lat/long,
   https://pypi.org/project/utm/
+* geopandas [JVFM21]_: used to work with geospatial data,
+  https://geopandas.org/
+* rtree [rtree]_: used for overlay operations in geopandas,
+  https://rtree.readthedocs.io/
 * openpyxl [GaCl18]_: used to read/write to Microsoft® Excel® spreadsheets,
   https://openpyxl.readthedocs.io
 * numpydoc [VaCV11]_: used to build the user manual,
   https://github.com/numpy/numpydoc
-* nose: used to run software tests,
-  http://nose.readthedocs.io
 
-These packages are included in the Anaconda Python distribution.
+All of these packages **except geopandas** are included in the Anaconda Python distribution.
 
+.. note:: 
+   Proper installation of geopandas requires installing several geopandas dependencies, including 
+   fiona, pyproj, and shapely.  See https://geopandas.org/en/stable/getting_started/install.html for more information.
+   On Linux and Mac OS X, installing geopandas through the conda-forge channel will install the dependencies.
+   On Windows, the dependencies must be installed manually, see https://geoffboeing.com/2014/09/using-geopandas-windows/ for more information.
+   
 .. The following is not shown in the UM
    WNTR includes a beta version of a Pyomo hydraulic simulator which requires installing 
    Pyomo, Interior Point OPTimizer (Ipopt), and HSL.
diff --git a/documentation/reference.rst b/documentation/reference.rst
index 14e092a95..422332ace 100644
--- a/documentation/reference.rst
+++ b/documentation/reference.rst
@@ -22,6 +22,8 @@ References
 
 .. [ELLT12] Eto, J.H., LaCommare, K.H., Larsen, P.H., Todd, A., and Fisher, E. (2012). An Examination of Temporal Trends in Electricity Reliability Based on Reports from U.S. Electric Utilities. Lawrence Berkeley National Laboratory Report Number LBNL-5268E. Berkeley, CA: Ernest Orlando Lawrence Berkeley National Laboratory, 68p
 
+.. [JVFM21] Jordahl, K., Van den Bossche, J., Fleischmann, M,  McBride, J. and others (2021) geopandas, 10.5281/zenodo.5573592
+
 .. [Folium] python-visualization/folium. (n.d.). Retrieved on February 5, 2019 from https://github.com/python-visualization/folium
 
 .. [GaCl18] Gazoni, E. and Clark, C. (2018) openpyxl - A Python library to read/write Excel 2010 xlsx/xlsm files, Retrieved on May 4, 2018 from https://openpyxl.readthedocs.io.
@@ -48,6 +50,8 @@ References
 
 .. [Ross00] Rossman, L.A. (2000). EPANET 2 Users Manual. Cincinnati, OH: U.S. Environmental Protection Agency. U.S. Environmental Protection Agency Technical Report, EPA/600/R--00/057, 200p.
 
+.. [rtree] Toblerity/rtree. (n.d.). Retrieved on March 17, 2022 from https://github.com/Toblerity/rtree
+
 .. [RWTS20] Rossman, L., Woo, H., Tryby M., Shang, F., Janke, R., and  Haxton, T. (2020) EPANET 2.2 User Manual. U.S. Environmental Protection Agency, Washington, DC, EPA/600/R-20/133.
 
 .. [SOKZ12] Salomons, E., Ostfeld, A., Kapelan, Z., Zecchin, A., Marchi, A., and Simpson, A. (2012). The battle of the water networks II - Problem description. Water Distribution Systems Analysis Conference 2012, September 24-27, Adelaide, South Australia, Australia. Retrieved on May 23, 2017 from https://emps.exeter.ac.uk/media/universityofexeter/emps/research/cws/downloads/WDSA2012-BWNII-ProblemDescription.pdf.
diff --git a/documentation/resilience.rst b/documentation/resilience.rst
index b2d6f30a0..83b2ec724 100644
--- a/documentation/resilience.rst
+++ b/documentation/resilience.rst
@@ -520,8 +520,8 @@ Economic metrics included in WNTR are listed in  :numref:`table-economic-metrics
                                           Greenhouse gas emissions can be computed 
                                           using the :class:`~wntr.metrics.economic.annual_ghg_emissions` method.
 
-   Pump operating power, energy and       The power, energy and cost required to operate a pump can be computed using the :class:`~wntr.metrics.economic.pump_energy`, 
-   cost                                   :class:`~wntr.metrics.economic.pump_energy` and 
+   Pump operating power, energy and cost  The power, energy and cost required to operate a pump can be computed using the :class:`~wntr.metrics.economic.pump_energy`, 
+                                          :class:`~wntr.metrics.economic.pump_energy` and 
                                           :class:`~wntr.metrics.economic.pump_cost` methods. These
                                           use the flowrates and pressures from simulation results to compute pump power, energy and cost.
    =====================================  ================================================================================================================================================
diff --git a/documentation/users.rst b/documentation/users.rst
index 5fb5192c9..331bee015 100644
--- a/documentation/users.rst
+++ b/documentation/users.rst
@@ -33,6 +33,8 @@ Publications
 
 * Antonowicz, A., Bałut, A., Urbaniak, A., & Zakrzewski, P. (2019). Algorithm for Early Warning System for Contamination in Water Network. In 2019 20th International Carpathian Control Conference (ICCC) (pp. 1-5). IEEE.
 
+* Antonowicz, A., & Urbaniak, A. (2022). Optimization of the process of restoring the continuity of the WDS based on the matrix and genetic algorithm approach. Bulletin of the Polish Academy of Sciences: Technical Sciences, e141594-e141594.
+
 * Bjerke, M. (2019). Leak Detection in Water Distribution Networks using Gated Recurrent Neural Networks, Master's thesis, Norwegian University of Science and Technology (NTNU)
 
 * Bunn, B. B. (2018). An operational model of interdependent water and power distribution infrastructure systems. Naval Postgraduate School, Monterey, CA
@@ -43,12 +45,22 @@ Publications
 
 * Han, Q., Mehrotra, S., & Venkatasubramanian, N. (2019). Aquaeis: Middleware support for event identification in community water infrastructures. In Proceedings of the 20th International Middleware Conference (pp. 293-305).
 
+* Huang, H., & Burton, H. V. (2022). Dynamic seismic damage assessment of distributed infrastructure systems using graph neural networks and semi-supervised machine learning. Advances in Engineering Software, 168, 103113.
+
+* Iannacone, L., Sharma, N., Tabandeh, A., & Gardoni, P. (2022). Modeling time-varying reliability and resilience of deteriorating infrastructure. Reliability Engineering & System Safety, 217, 108074.
+
+* Kammoun, M., Kammoun, A., & Abid, M. (2022). Experiments based comparative evaluations of machine learning techniques for leak detection in water distribution systems. Water Supply, 22(1), 628-642.
+
+* Liu, J., & Kang, Y. (2022). Segment-based resilience response and intervention evaluation of water distribution systems. AQUA—Water Infrastructure, Ecosystems and Society, 71(1), 100-119.
+
 * Liu, Y., Barrows, C., Macknick, J., & Mauter, M. (2020). Optimization Framework to Assess the Demand Response Capacity of a Water Distribution System. Journal of Water Resources Planning and Management, 146(8), 04020063.
 
 * Logan, K. T., Lestakova, M., Thiessen, N., Engels, J. I., & Pelz, P. F. (2021). Water Distribution in a Socio-Technical System: Resilience Assessment for Critical Events Causing Demand Relocation. Water, 13(15), 2062.
 
 * Lorenz, I. S., & Pelz, P. F. (2020). Optimal resilience enhancement of water distribution systems. Water, 12(9), 2602.
 
+* Marlim, M. S., & Kang, D. (2022). Contaminant Flushing in Water Distribution Networks Incorporating Customer Faucet Control. Sustainability, 14(4), 2249.
+
 * Mazumder, R. K., Salman, A. M., Li, Y., & Yu, X. (2019). A Decision-making Framework for Water Distribution Systems using Fuzzy Inference and Centrality Analysis. 13th International Conference on Applications of Statistics and Probability in Civil Engineering, ICASP13, Seoul, South Korea, May 26-30, 2019
 
 * Mazumder, R. K., Salman, A. M., & Li, Y. (2020). Post-disaster sequential recovery planning for water distribution systems using topological and hydraulic metrics. Structure and Infrastructure Engineering, 1-16.
@@ -59,16 +71,26 @@ Publications
 
 * Nikolopoulos, D., Ostfeld, A., Salomons, E., & Makropoulos, C. (2021). Resilience Assessment of Water Quality Sensor Designs under Cyber-Physical Attacks. Water, 13(5), 647.
 
+* Nikolopoulos, D., Kossieris, P., Tsoukalas, I., & Makropoulos, C. (2022). Stress-testing framework for urban water systems: A source to tap approach for stochastic resilience assessment. Water, 14(2), 154.
+
+* Nikolopoulos, D., & Makropoulos, C. (2022). Stress-testing water distribution networks for cyber-physical attacks on water quality. Urban Water Journal, 19(3), 256-270.
+
 * Nyahora, P. P., Babel, M. S., Ferras, D., & Emen, A. (2020). Multi-objective optimization for improving equity and reliability in intermittent water supply systems. Water Supply, 20(5), 1592-1603.
 
 * Pagani, A., Wei, Z., Silva, R., & Guo, W. (2020). Neural Network Approximation of Graph Fourier Transforms for Sparse Sampling of Networked Flow Dynamics. arXiv preprint arXiv:2002.05508.
 
+* Rahimi-Golkhandan, A., Aslani, B., & Mohebbi, S. (2022). Predictive resilience of interdependent water and transportation infrastructures: A sociotechnical approach. Socio-Economic Planning Sciences, 80, 101166.
+
+* Randeniya, A., Radhakrishnan, M., Sirisena, T. A. J. G., Maish, I., & Pathirana, A. (2022). Equity–performance trade-off in water rationing regimes with domestic storage. Water Supply, 22(5), 4781-4797.
+
 * Sharma, N., Tabandeh, A., & Gardoni, P. (2019). Recovery optimization of interdependent infrastructure: a multi-scale approach. 13th International Conference on Applications of Statistics and Probability in Civil Engineering, ICASP13
 
 * Sharma, N., Tabandeh, A., & Gardoni, P. (2020). Regional resilience analysis: A multiscale approach to optimize the resilience of interdependent infrastructure. Computer‐Aided Civil and Infrastructure Engineering, 35(12), 1315-1330.
 
 * Tabandeh, S. (2018). Societal risk and resilience analysis: A multi-scale approach to model the dynamics of infrastructure-social systems (Doctoral dissertation, University of Illinois at Urbana-Champaign).
 
+* Tabandeh, A., Sharma, N., & Gardoni, P. (2022). Uncertainty propagation in risk and resilience analysis of hierarchical systems. Reliability Engineering & System Safety, 219, 108208.
+
 * Tomar, A., Burton, H. V., Mosleh, A., & Yun Lee, J. (2020). Hindcasting the Functional Loss and Restoration of the Napa Water System Following the 2014 Earthquake Using Discrete-Event Simulation. Journal of Infrastructure Systems, 26(4), 04020035.
 
 * Vrachimis, S. G., & Kyriakou, M. S. (2018). LeakDB: A benchmark dataset for leakage diagnosis in water distribution networks. In WDSA/CCWI Joint Conference Proceedings (Vol. 1).
diff --git a/documentation/waternetworkmodel.rst b/documentation/waternetworkmodel.rst
index fd40a4673..c1a1c0bc8 100644
--- a/documentation/waternetworkmodel.rst
+++ b/documentation/waternetworkmodel.rst
@@ -2,6 +2,15 @@
 
     \clearpage
 
+.. doctest::
+    :hide:
+	
+    >>> try:
+    ...    import geopandas as gpd
+    ... except ModuleNotFoundError:
+    ...    gpd = None
+
+
 Water network model
 ======================================
 
@@ -263,21 +272,29 @@ When using the EpanetSimulator, this step is not needed since EPANET starts at t
 
     >>> wn.reset_initial_values()
 
-Write a model to an INP file
----------------------------------
+Read and write a model to an INP file
+--------------------------------------
+
+As described above, the water network model can be created from an INP file.
+This can be accomplished by creating the `WaterNetworkModel` with an INP file (as shown above), 
+or using the function :class:`~wntr.network.io.read_inpfile`.
 
-The water network model can be written to a file in EPANET INP format.
+.. doctest::
 
+    >>> wn = wntr.network.read_inpfile('networks/Net3.inp') # doctest: +SKIP
+	
+The water network model can also be written to a file in EPANET INP format
+using the :class:`~wntr.network.io.write_inpfile` function.
 By default, files are written in the LPS (liter per second) EPANET unit convention.
 The EPANET INP file will not include features not supported by EPANET (i.e., custom element attributes).
 EPANET INP files can be saved in EPANET 2.00.12 or 2.2.0 format.
 
 .. doctest::
 
-    >>> wn.write_inpfile('filename.inp', version=2.2)
+    >>> wntr.network.write_inpfile(wn, 'filename.inp', version=2.2)
 
 Read and write a model to a dictionary or JSON file
----------------------------------------------------
+---------------------------------------------------------
 
 The water network model can be converted to a dictionary representation.
 The dictionary contains keys for each of the following water network model objects: 
@@ -293,16 +310,47 @@ To create a water network model from a dictionary, use the :class:`~wntr.network
 
 The water network model can also be converted to a JSON (JavaScript Object Notation) file.  
 The JSON file is a formatted version of the dictionary representation.
-The :class:`~wntr.network.model.WaterNetworkModel.write_json` and :class:`~wntr.network.model.WaterNetworkModel.read_json` methods wrap 
+The :class:`~wntr.network.io.write_json` and :class:`~wntr.network.io.read_json` methods wrap 
 the dictionary representation with the Python standard library JSON loader and dumper.
 
 .. doctest::
 
-    >>> wn.write_json('Net3.json')
+    >>> wntr.network.write_json(wn, 'Net3.json')
 	
 Note that these methods do not check for a valid dictionary/JSON schema prior to building a model.
 They simply ignore extraneous or invalid dictionary keys.
 
+.. _giswrite:
+
+Read and write a model to GeoDataFrames, GeoJSON files, or Shapefiles
+----------------------------------------------------------------------
+
+The water network model can also be converted to GeoDataFrames and GIS formatted files, including 
+GeoJSON files and Shapefiles. 
+See :ref:`geospatial` for more information.  
+
+To create GeoDataFrames, use the :class:`~wntr.network.model.WaterNetworkModel.to_gis` method on the water network model.
+To create a water network model from GeoDataFrames, use the :class:`~wntr.network.model.WaterNetworkModel.from_gis` method.
+
+.. doctest::
+    :skipif: gpd is None
+
+    >>> wn_gis = wn.to_gis() # doctest: +SKIP
+
+The water network model can also be saved as a series of GeoJSON files and Shapefiles and 
+created from a series of GeoJSON files and Shapefiles using the 
+:class:`~wntr.network.io.write_geojson`, 
+:class:`~wntr.network.io.read_geojson`, 
+:class:`~wntr.network.io.write_shapefile`, and 
+:class:`~wntr.network.io.read_shapefile`
+functions.  A file or directory prefix is given to name the files.
+
+.. doctest::
+
+    >>> wntr.network.write_geojson(wn, 'Net3') # doctest: +SKIP
+
+Note, the GeoPandas ``read_file`` method can also be used to read GeoJSON and Shapefiles into GeoDataFrames.
+   
 Build a model from scratch
 ---------------------------------
 
diff --git a/documentation/whatsnew/v0.5.0.rst b/documentation/whatsnew/v0.5.0.rst
index 05f8cac53..076cb0067 100644
--- a/documentation/whatsnew/v0.5.0.rst
+++ b/documentation/whatsnew/v0.5.0.rst
@@ -2,11 +2,29 @@
 
 v0.5.0 (main)
 ---------------------------------------------------
+* Added GIS capabilities in a `wntr.gis` module, which includes the following:
 
-* Added support for Python 3.10, dropped support for Python 3.6.  
+  * Convert between water network models and GIS formatted data (GeoPandas) 
+    and GIS formatted files (GeoJSON and Shapefiles)
+  
+  * Snap GeoPandas point data to elements in a water network model
+  
+  * Find the intersection between GeoPandas points, line or polygons and 
+    elements in a water network model
+  
+  * See :ref:`geospatial` section of the user manual for more information.
 
-* The build workflow now creates wheel artifacts.  
-  The `_evaluator` and `_network_isolation` binaries are also distributed with the code.
+* Removed the `read_inpfile`, `write_inpfile`, `read_json`, and `write_json` methods from the WaterNetworkModel.  
+  This functionality is still available using `wntr.network.read_inpfile`, `wntr.network.write_inpfile`, 
+  `wntr.network.read_json` and `wntr.network.write_json`
+  and is located with other read/write functions.
 
-* Updated documentation
+* Added support for Python 3.10, dropped support for Python 3.6.
+  The build workflow now creates wheel artifacts.  
+  The `_evaluator` and `_network_isolation` binaries are also 
+  distributed with the code. `#287 <https://github.com/USEPA/WNTR/pull/287>`_, 
+  `#294 <https://github.com/USEPA/WNTR/pull/294>`_
+
+* Updated documentation `#291 <https://github.com/USEPA/WNTR/pull/291>`_, 
+  `#292 <https://github.com/USEPA/WNTR/pull/292>`_
  
\ No newline at end of file
diff --git a/examples/data/Net1_demographic_data.geojson b/examples/data/Net1_demographic_data.geojson
new file mode 100644
index 000000000..933a6dcc1
--- /dev/null
+++ b/examples/data/Net1_demographic_data.geojson
@@ -0,0 +1,13 @@
+{
+"type": "FeatureCollection",
+"features": [
+{ "type": "Feature", "properties": { "mean_income": 63326.0, "mean_age": 35.0, "population": 3362.0 }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 41.678134161727016, 82.750228129522114 ], [ 41.985961673378931, 60.857792546547728 ], [ 23.210840621475665, 40.191596652399582 ], [ 9.577878109518823, 39.643589069759734 ], [ 0.513167019494862, 66.837722339831615 ], [ 0.255395387480345, 67.754408553262493 ], [ 0.085983915505732, 68.691456887837106 ], [ 0.006468758743196, 69.640370565687036 ], [ 0.017570928691562, 70.592545218777573 ], [ 0.119189755333684, 71.539346909860541 ], [ 0.310403799971716, 72.47219042154876 ], [ 0.589479210463072, 73.382617103618628 ], [ 0.953885443095125, 74.262371572652214 ], [ 1.400318208539037, 75.103476568537843 ], [ 1.548522088801759, 75.328103149448637 ], [ 41.678134161727016, 82.750228129522114 ] ] ] } },
+{ "type": "Feature", "properties": { "mean_income": 78245.0, "mean_age": 31.0, "population": 5618.0 }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 23.210840621475665, 40.191596652399582 ], [ 22.990631170537007, 27.717773173337257 ], [ 18.873500743899402, 11.756721166618 ], [ 9.577878109518823, 39.643589069759734 ], [ 23.210840621475665, 40.191596652399582 ] ] ] } },
+{ "type": "Feature", "properties": { "mean_income": 91452.0, "mean_age": 40.0, "population": 5650.0 }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 22.990631170537007, 27.717773173337257 ], [ 61.937201479083768, 16.361652197633106 ], [ 59.380301104236082, 6.042695279034173 ], [ 58.320502943378436, 4.452998037747709 ], [ 57.736073394141762, 3.663347221091248 ], [ 57.07698251241132, 2.934851840263761 ], [ 56.349591225877091, 2.274542650155613 ], [ 55.560919633745755, 1.688792336427817 ], [ 54.718579255320243, 1.183254012320562 ], [ 53.830699570746681, 0.762806660100635 ], [ 52.905849562888463, 0.431508043695656 ], [ 51.952955017531934, 0.192555546958381 ], [ 50.981212380064534, 0.04825531551336 ], [ 50.0, 0.0 ], [ 30.0, 0.0 ], [ 29.04067314748271, 0.046121761341499 ], [ 28.090195463793474, 0.184061601992147 ], [ 27.157334489909807, 0.412547116269948 ], [ 26.250695264072156, 0.729470673303004 ], [ 25.378640945879916, 1.131908858559115 ], [ 24.549215671561164, 1.616149440451267 ], [ 23.770070352026465, 2.177725613269478 ], [ 23.048392098171398, 2.811457200570191 ], [ 22.3908379244339, 3.511498438948591 ], [ 21.803473342148003, 4.271391901419545 ], [ 21.291716409130647, 5.084128062998836 ], [ 20.860287751607721, 5.942209959030879 ], [ 20.513167019494862, 6.83772233983162 ], [ 18.873500743899402, 11.756721166618 ], [ 22.990631170537007, 27.717773173337257 ] ] ] } },
+{ "type": "Feature", "properties": { "mean_income": 54040.0, "mean_age": 39.0, "population": 5546.0 }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 61.937201479083768, 16.361652197633106 ], [ 22.990631170537007, 27.717773173337257 ], [ 23.210840621475665, 40.191596652399582 ], [ 41.985961673378931, 60.857792546547728 ], [ 51.729942246639531, 45.92347392620681 ], [ 64.044563367543333, 22.10119024305617 ], [ 61.937201479083768, 16.361652197633106 ] ] ] } },
+{ "type": "Feature", "properties": { "mean_income": 26135.0, "mean_age": 38.0, "population": 5968.0 }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 61.937201479083768, 16.361652197633106 ], [ 64.044563367543333, 22.10119024305617 ], [ 71.664299776493777, 24.468693287420724 ], [ 59.380301104236082, 6.042695279034173 ], [ 61.937201479083768, 16.361652197633106 ] ] ] } },
+{ "type": "Feature", "properties": { "mean_income": 57620.0, "mean_age": 31.0, "population": 4315.0 }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 44.484966890794844, 87.214865127957253 ], [ 79.811444522502882, 71.926689821038252 ], [ 79.95184726672197, 70.9801714032956 ], [ 80.0, 70.0 ], [ 80.0, 51.004078682461177 ], [ 51.729942246639531, 45.92347392620681 ], [ 41.985961673378931, 60.857792546547728 ], [ 41.678134161727016, 82.750228129522114 ], [ 44.484966890794844, 87.214865127957253 ] ] ] } },
+{ "type": "Feature", "properties": { "mean_income": 44871.0, "mean_age": 54.0, "population": 4547.0 }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 64.044563367543333, 22.10119024305617 ], [ 51.729942246639531, 45.92347392620681 ], [ 80.0, 51.004078682461177 ], [ 80.0, 40.0 ], [ 79.952017995138718, 39.021563581812764 ], [ 79.80853243511298, 38.052516631820303 ], [ 79.570920264890532, 37.102158513115704 ], [ 79.241461706132696, 36.179609243273276 ], [ 78.823318375273061, 35.293721975003102 ], [ 78.320502943378443, 34.452998037747712 ], [ 71.664299776493777, 24.468693287420724 ], [ 64.044563367543333, 22.10119024305617 ] ] ] } },
+{ "type": "Feature", "properties": { "mean_income": 69067.0, "mean_age": 55.0, "population": 2541.0 }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 46.010471610447148, 99.157253485762936 ], [ 46.406538345668572, 99.332043363532435 ], [ 47.318360002498913, 99.633732761697431 ], [ 48.25491805132296, 99.84655721520987 ], [ 49.207573350778652, 99.968553556339245 ], [ 50.167538272793557, 99.998596447859541 ], [ 51.125957762761786, 99.936408763556202 ], [ 52.073991021687078, 99.782564144535996 ], [ 53.002893056827034, 99.53848170775936 ], [ 53.904095348586786, 99.206412955605614 ], [ 54.769284889562144, 98.789421007221932 ], [ 55.590480866648306, 98.291352343233228 ], [ 56.360108278871309, 97.71680132444996 ], [ 57.071067811865476, 97.071067811865476 ], [ 77.071067811865476, 77.071067811865476 ], [ 77.730104533627369, 76.343932841636459 ], [ 78.314696123025456, 75.555702330196027 ], [ 78.819212643483553, 74.713967368259972 ], [ 79.238795325112875, 73.8268343236509 ], [ 79.569403357322088, 72.902846772544621 ], [ 79.807852804032308, 71.95090322016128 ], [ 79.811444522502882, 71.926689821038252 ], [ 44.484966890794844, 87.214865127957253 ], [ 46.010471610447148, 99.157253485762936 ] ] ] } }
+]
+}
diff --git a/examples/data/Net1_earthquake_data.geojson b/examples/data/Net1_earthquake_data.geojson
new file mode 100644
index 000000000..4acc0f4c2
--- /dev/null
+++ b/examples/data/Net1_earthquake_data.geojson
@@ -0,0 +1,9 @@
+{
+"type": "FeatureCollection",
+"features": [
+{ "type": "Feature", "properties": { "Pr": 0.5 }, "geometry": { "type": "LineString", "coordinates": [ [ 36.0, 2.0 ], [ 44.0, 44.0 ], [ 85.0, 85.0 ] ] } },
+{ "type": "Feature", "properties": { "Pr": 0.75 }, "geometry": { "type": "LineString", "coordinates": [ [ 42.0, 2.0 ], [ 45.0, 27.0 ], [ 38.0, 56.0 ], [ 30.0, 85.0 ] ] } },
+{ "type": "Feature", "properties": { "Pr": 0.9 }, "geometry": { "type": "LineString", "coordinates": [ [ 40.0, 2.0 ], [ 50.0, 50.0 ], [ 60.0, 85.0 ] ] } },
+{ "type": "Feature", "properties": { "Pr": 0.25 }, "geometry": { "type": "LineString", "coordinates": [ [ 30.0, 2.0 ], [ 35.0, 30.0 ], [ 40.0, 50.0 ], [ 60.0, 80.0 ] ] } }
+]
+}
diff --git a/examples/data/Net1_hydrant_data.geojson b/examples/data/Net1_hydrant_data.geojson
new file mode 100644
index 000000000..71f323f19
--- /dev/null
+++ b/examples/data/Net1_hydrant_data.geojson
@@ -0,0 +1,8 @@
+{
+"type": "FeatureCollection",
+"features": [
+{ "type": "Feature", "properties": { "demand": 5000 }, "geometry": { "type": "Point", "coordinates": [ 48.2, 37.2 ] } },
+{ "type": "Feature", "properties": { "demand": 1500 }, "geometry": { "type": "Point", "coordinates": [ 71.8, 68.3 ] } },
+{ "type": "Feature", "properties": { "demand": 8000 }, "geometry": { "type": "Point", "coordinates": [ 51.2, 71.1 ] } }
+]
+}
diff --git a/examples/data/Net1_landslide_data.geojson b/examples/data/Net1_landslide_data.geojson
new file mode 100644
index 000000000..ff96ae537
--- /dev/null
+++ b/examples/data/Net1_landslide_data.geojson
@@ -0,0 +1,8 @@
+{
+"type": "FeatureCollection",
+"features": [
+{ "type": "Feature", "properties": { "Pr": 0.5 }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 28.846153846153847, 22.23076923076923 ], [ 28.76039947938807, 22.050793134785565 ], [ 28.657153252394828, 21.880248644529846 ], [ 28.537441050828335, 21.720830341499664 ], [ 28.402452371232886, 21.574122254040251 ], [ 28.253528501846422, 21.441582117961083 ], [ 28.092149195155997, 21.324526892033145 ], [ 27.919917964630738, 21.224119672289362 ], [ 27.738546151728322, 21.14135813515156 ], [ 27.549835921490054, 21.077064624216614 ], [ 27.355662355685425, 21.031877979202434 ], [ 27.157954821434057, 21.006247188243805 ], [ 26.958677800432024, 21.000426926610867 ], [ 26.759811369269169, 21.014475026178808 ], [ 26.563331524790847, 21.048251900792902 ], [ 26.371190549997237, 21.101421933238662 ], [ 26.185297615570569, 21.173456810035688 ], [ 26.007499809779347, 21.263640770919807 ], [ 25.839563785252171, 21.371077720852917 ], [ 25.683158204984334, 21.494700133893609 ], [ 25.539837161998971, 21.63327966045691 ], [ 25.411024737409999, 21.785439332566281 ], [ 25.29800085032263, 21.949667245822777 ], [ 25.201888540171062, 22.124331582143192 ], [ 25.12364280785992, 22.30769682399675 ], [ 25.064041126587188, 22.497940999030796 ], [ 25.023675716636024, 22.693173783737659 ], [ 25.002947660895369, 22.891455286279207 ], [ 25.002062919579402, 23.090815321837301 ], [ 25.993211742958316, 44.896089436173433 ], [ 23.022127294263406, 64.703319094139516 ], [ 23.002742283780481, 64.895302268402844 ], [ 23.001948406132485, 65.088260003645388 ], [ 23.019753050994115, 65.280396185696034 ], [ 23.055990486865117, 65.469922347671471 ], [ 23.110323403753689, 65.65507431761516 ], [ 23.182246052969035, 65.834128639991064 ], [ 23.27108895479676, 66.005418618175653 ], [ 23.376025130236659, 66.167349828619137 ], [ 23.496077798795685, 66.318414962265209 ], [ 23.63012947068259, 66.45720785508027 ], [ 23.776932348771066, 66.582436577091585 ], [ 23.935119943506496, 66.692935458097082 ], [ 24.103219792640928, 66.78767593810764 ], [ 24.279667167396791, 66.865776141522304 ], [ 24.462819637477875, 66.926509085916962 ], [ 24.650972359351638, 66.969309449036245 ], [ 24.842373945494472, 66.993778830999318 ], [ 25.035242766883794, 66.999689462737251 ], [ 25.227783536987953, 66.9869863261425 ], [ 25.418204022884741, 66.955787666195647 ], [ 25.60473172795589, 66.906383890302124 ], [ 25.785630390869578, 66.839234865084421 ], [ 25.959216147272901, 66.754965635792033 ], [ 26.123873203756062, 66.654360608174372 ], [ 26.278068878189739, 66.538356246973635 ], [ 26.420367866435324, 66.408032359002434 ], [ 26.549445602628698, 66.264602041946233 ], [ 26.664100588675687, 66.109400392450453 ], [ 32.664100588675687, 57.109400392450461 ], [ 32.757143617226482, 56.955220554866898 ], [ 32.83594129468684, 56.793296642160804 ], [ 32.899854799864492, 56.624941388797254 ], [ 32.948365979142544, 56.451519668807364 ], [ 32.981081547206202, 56.274437430608636 ], [ 32.997736275448872, 56.095130298830895 ], [ 33.997736275448872, 35.095130298830895 ], [ 33.998328098481785, 34.918239307621683 ], [ 33.983287613774685, 34.741987905229372 ], [ 33.952732478425865, 34.567754852310884 ], [ 33.906901715886569, 34.396903120594331 ], [ 33.846153846153847, 34.230769230769234 ], [ 28.846153846153847, 22.23076923076923 ] ] ] } },
+{ "type": "Feature", "properties": { "Pr": 0.75 }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 40.007075068197658, 1.831921993703417 ], [ 33.007075068197658, 84.831921993703418 ], [ 33.00020764471715, 85.028819017198018 ], [ 33.012746005956245, 85.225436383836197 ], [ 33.044568480913313, 85.41986613838499 ], [ 33.095366267478646, 85.610221553980011 ], [ 33.164646429019257, 85.794655440734047 ], [ 33.251736677786532, 85.971378070682647 ], [ 33.355791898729272, 86.138674545125028 ], [ 33.475802350405438, 86.29492143582938 ], [ 33.610603463412041, 86.438602538617729 ], [ 33.758887141249915, 86.568323586459499 ], [ 33.919214453961189, 86.682825779299094 ], [ 34.090029601361451, 86.78099799932555 ], [ 34.269675010368893, 86.861887593148495 ], [ 34.456407419927515, 86.924709616251278 ], [ 34.648414797437781, 86.968854450013836 ], [ 34.843833922539488, 86.993893717390819 ], [ 35.040768467615528, 86.999584439839509 ], [ 35.237307399565132, 86.985871395159222 ], [ 35.431543524277686, 86.952887653362069 ], [ 35.621591993853379, 86.90095328537484 ], [ 35.805608596978388, 86.830572257102816 ], [ 35.981807654966808, 86.742427538995116 ], [ 36.148479349807744, 86.637374478567807 ], [ 36.304006316068026, 86.516432500197311 ], [ 36.446879335644411, 86.380775212728409 ], [ 36.575711983064856, 86.231719020891461 ], [ 36.689254079222472, 86.07071035104282 ], [ 36.786403822988341, 85.899311615188324 ], [ 36.866218482979335, 85.719186049493672 ], [ 36.927923545728994, 85.532081574407286 ], [ 36.970920231488698, 85.339813833015299 ], [ 46.970920231488698, 27.339813833015292 ], [ 46.994524452920331, 27.14789255120812 ], [ 46.999484389471924, 26.954588809109559 ], [ 46.985753676973843, 26.76170955876314 ], [ 43.985753676973843, 1.761709558763139 ], [ 43.95293900721969, 1.568688935825005 ], [ 43.901408461991302, 1.379801756961632 ], [ 43.831655881221231, 1.19685821128042 ], [ 43.744349734641951, 1.021611527431689 ], [ 43.640326715545719, 0.855741171645598 ], [ 43.520583722383876, 0.700836752666857 ], [ 43.386268305049498, 0.558382787833333 ], [ 43.238667667400875, 0.429744476291942 ], [ 43.079196331419467, 0.316154615693365 ], [ 42.90938258122214, 0.218701787748676 ], [ 42.730853816840707, 0.13831992587089 ], [ 42.54532095813002, 0.075779364879339 ], [ 42.354562048268662, 0.031679458541488 ], [ 42.160405213987367, 0.006442835701553 ], [ 41.964711145824573, 0.000311350041966 ], [ 41.7693552663087, 0.013343762292915 ], [ 41.576209756957283, 0.045415177102362 ], [ 41.3871256163361, 0.096218239963259 ], [ 41.203914921123911, 0.165266082727293 ], [ 41.028333460182772, 0.251896989476988 ], [ 40.862063908059305, 0.355280738041067 ], [ 40.706699699172709, 0.474426556379525 ], [ 40.563729757230234, 0.608192617588882 ], [ 40.434524226214918, 0.755296982532783 ], [ 40.320321339691709, 0.914329885229888 ], [ 40.222215554269383, 1.083767243262782 ], [ 40.14114706094032, 1.261985263731686 ], [ 40.077892774815759, 1.447276004777674 ], [ 40.033058889605762, 1.637863743542464 ], [ 40.007075068197658, 1.831921993703417 ] ] ] } },
+{ "type": "Feature", "properties": { "Pr": 0.9 }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 58.059714999709335, 44.485071250072664 ], [ 58.117760887496061, 44.676147856138279 ], [ 58.194572508410467, 44.860483336626999 ], [ 58.289384058379298, 45.036239885486523 ], [ 58.401250275861493, 45.201665227743177 ], [ 58.529055866009912, 45.3551100895057 ], [ 58.671526620047864, 45.495044641058776 ], [ 58.82724211900149, 45.620073749110169 ], [ 58.994649895132859, 45.728950886128246 ], [ 59.172080909885985, 45.820590558095034 ], [ 59.357766194031967, 45.894079126771587 ], [ 59.549854484112913, 45.948683918578574 ], [ 59.746430679350873, 45.983860529277585 ], [ 59.945534935008503, 45.999258251626209 ], [ 60.14518220184236, 45.994723571893658 ], [ 60.343382016842611, 45.970301700377156 ], [ 60.538158347947487, 45.92623612065978 ], [ 60.727569294882557, 45.862966162103888 ], [ 60.909726449709623, 45.781122619781897 ], [ 61.082813724062405, 45.681521465513335 ], [ 61.245105455363422, 45.56515571270856 ], [ 61.394983611504969, 45.433185516125718 ], [ 61.530953922465166, 45.286926605245426 ], [ 61.651660778028607, 45.127837166581216 ], [ 61.755900743082734, 44.957503305707874 ], [ 61.842634555743956, 44.777623233949633 ], [ 61.910997488693233, 44.589990337385416 ], [ 61.960307970419791, 44.396475296971502 ], [ 61.990074380419976, 44.199007438041995 ], [ 62.990074380419976, 34.199007438041995 ], [ 62.999997826116257, 34.002948818450129 ], [ 62.990652564211146, 33.806861788867138 ], [ 62.962128630362059, 33.612635523165181 ], [ 62.914700834538195, 33.422141267939352 ], [ 62.848826113401167, 33.23721431423634 ], [ 62.765139128010773, 33.059636315691975 ], [ 62.664446149269153, 32.891118123431042 ], [ 62.547717290013196, 32.733283303104351 ], [ 62.416077158594106, 32.587652492865857 ], [ 62.270794023990895, 32.455628752990712 ], [ 62.113267596844999, 32.338484048281522 ], [ 61.945015544138116, 32.237346993496068 ], [ 61.767658867435649, 32.153191979861631 ], [ 61.582906285567155, 32.086829787434922 ], [ 61.392537772206666, 32.03889977375173 ], [ 61.198387406957941, 32.00986371402346 ], [ 61.002325705163479, 32.000001352226583 ], [ 60.806241596678277, 32.009407705947247 ], [ 60.612024227231565, 32.037992150947211 ], [ 60.421544757708688, 32.085479294270833 ], [ 60.236638336705489, 32.151411627481053 ], [ 60.059086420038284, 32.235153934462382 ], [ 59.890599607549539, 32.335899411324313 ], [ 56.890599607549539, 34.335899411324313 ], [ 56.738612849156254, 34.447936065850925 ], [ 56.597868893308771, 34.573806338659146 ], [ 56.469622197454228, 34.712388342132932 ], [ 56.3550158303782, 34.86244688840754 ], [ 56.255071283979241, 35.022644498656632 ], [ 56.170679368663919, 35.191553324103459 ], [ 56.102592273512109, 35.367667872503652 ], [ 56.051416861980599, 35.549418426667856 ], [ 56.017609262900784, 35.735185035424308 ], [ 56.001470804981096, 35.923311952319231 ], [ 56.00314533105017, 36.112122393362249 ], [ 56.022617915979012, 36.299933482280352 ], [ 56.059714999709335, 36.485071250072664 ], [ 58.059714999709335, 44.485071250072664 ] ] ] } }
+]
+}
diff --git a/examples/data/Net1_valve_data.geojson b/examples/data/Net1_valve_data.geojson
new file mode 100644
index 000000000..d130a7663
--- /dev/null
+++ b/examples/data/Net1_valve_data.geojson
@@ -0,0 +1,8 @@
+{
+"type": "FeatureCollection",
+"features": [
+{ "type": "Feature", "properties": { }, "geometry": { "type": "Point", "coordinates": [ 56.5, 41.5 ] } },
+{ "type": "Feature", "properties": { }, "geometry": { "type": "Point", "coordinates": [ 32.1, 67.6 ] } },
+{ "type": "Feature", "properties": { }, "geometry": { "type": "Point", "coordinates": [ 52.7, 86.3 ] } }
+]
+}
diff --git a/requirements.txt b/requirements.txt
index acbb38847..f4c7bd03a 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -10,9 +10,10 @@ plotly
 folium
 utm
 openpyxl
+geopandas; sys_platform == "darwin" or sys_platform == "linux"
+rtree; sys_platform == "darwin" or sys_platform == "linux"
 
 # Testing
-# nose
 pytest
 ipykernel
 nbformat
diff --git a/wntr/__init__.py b/wntr/__init__.py
index b61fe9579..930b80d76 100644
--- a/wntr/__init__.py
+++ b/wntr/__init__.py
@@ -5,6 +5,7 @@
 from wntr import sim
 from wntr import scenario
 from wntr import graphics
+from wntr import gis
 from wntr import utils
 
 __version__ = '0.5.0'
diff --git a/wntr/gis/__init__.py b/wntr/gis/__init__.py
new file mode 100644
index 000000000..babbb2afb
--- /dev/null
+++ b/wntr/gis/__init__.py
@@ -0,0 +1,7 @@
+"""
+The wntr.gis package contains methods to convert between water network models
+and GIS formatted data and geospatial functions to snap data and find intersections.
+"""
+from wntr.gis.network import WaterNetworkGIS
+from wntr.gis.geospatial import snap, intersect
+
diff --git a/wntr/gis/geospatial.py b/wntr/gis/geospatial.py
new file mode 100644
index 000000000..b31d0ef04
--- /dev/null
+++ b/wntr/gis/geospatial.py
@@ -0,0 +1,271 @@
+"""
+The wntr.gis.geospatial module contains functions to snap data and find 
+intersects with polygons.
+"""
+
+import pandas as pd
+import numpy as np
+
+try:
+    from shapely.geometry import MultiPoint, LineString, Point, shape
+    has_shapely = True
+except ModuleNotFoundError:
+    has_shapely = False
+
+try:
+    import geopandas as gpd
+    has_geopandas = True
+except ModuleNotFoundError:
+    gpd = None
+    has_geopandas = False
+
+
+def snap(A, B, tolerance):  
+    """
+    Snap Points in A to Points or Lines in B
+
+    For each Point geometry in A, the function returns snapped Point geometry 
+    and associated element in B. Note the CRS of A must equal the CRS of B.
+    
+    Parameters
+    ----------
+    A : geopandas GeoDataFrame
+        GeoDataFrame containing Point geometries.
+    B : geopandas GeoDataFrame
+        GeoDataFrame containing Point, LineString, or MultiLineString geometries.
+    tolerance : float
+        Maximum allowable distance (in the coordinate reference system units) 
+        between Points in A and Points or Lines in B.  
+    
+    Returns
+    -------
+    GeoPandas GeoDataFrame
+        Snapped points (index = A.index, columns = defined below)
+        
+        If B contains Points, columns include:
+            - node: closest Point in B to Point in A
+            - snap_distance: distance between Point in A and snapped point
+            - geometry: GeoPandas Point object of the snapped point
+        
+        If B contains Lines or MultiLineString, columns include:
+            - link: closest Line in B to Point in A
+            - node: start or end node of Line in B that is closest to the snapped point
+            - snap_distance: distance between Point A and snapped point
+            - line_position: normalized distance of snapped point along Line in B from the start node (0.0) and end node (1.0)
+            - geometry: GeoPandas Point object of the snapped point
+    """   
+    if not has_shapely or not has_geopandas:
+        raise ModuleNotFoundError('shapley and geopandas are required')
+        
+    isinstance(A, gpd.GeoDataFrame)
+    assert(A['geometry'].geom_type).isin(['Point']).all()
+    isinstance(B, gpd.GeoDataFrame)
+    assert (B['geometry'].geom_type).isin(['Point', 'LineString', 'MultiLineString']).all()
+    assert A.crs == B.crs
+    
+    # Modify B to include "indexB" as a separate column
+    B = B.reset_index()
+    B.rename(columns={'index':'indexB'}, inplace=True)
+    
+    # Define the coordinate reference system, based on B
+    crs = B.crs
+    
+    # Determine which Bs are closest to each A
+    bbox = A.bounds + [-tolerance, -tolerance, tolerance, tolerance]       
+    hits = bbox.apply(lambda row: list(B.loc[list(B.sindex.intersection(row))]['indexB']), axis=1)        
+    closest = pd.DataFrame({
+        # index of points table
+        "point": np.repeat(hits.index, hits.apply(len)),
+        # name of link
+        "indexB": np.concatenate(hits.values)
+        })
+    
+    # Merge the closest dataframe with the lines dataframe on the line names
+    closest = pd.merge(closest, B, on="indexB")
+
+    # Join back to the original points to get their geometry
+    # rename the point geometry as "points"
+    closest = closest.join(A.geometry.rename("points"), on="point")
+    
+    # Convert back to a GeoDataFrame, so we can do spatial ops
+    closest = gpd.GeoDataFrame(closest, geometry="geometry", crs=crs)  
+    
+    # Calculate distance between the point and nearby links
+    closest["snap_distance"] = closest.geometry.distance(gpd.GeoSeries(closest.points, crs=crs))
+    
+    # Sort on ascending snap distance, so that closest goes to top
+    closest = closest.sort_values(by=["snap_distance"]) 
+       
+    # group by the index of the points and take the first, which is the closest line
+    closest = closest.groupby("point").first()      
+    
+    # construct a GeoDataFrame of the closest elements of B
+    closest = gpd.GeoDataFrame(closest, geometry="geometry", crs=crs)
+    
+    # Reset B index
+    B.set_index('indexB', inplace=True)
+    B.index.name = None
+    
+    # snap to points
+    if B['geometry'].geom_type.isin(['Point']).all():
+        snapped_points = closest.rename(columns={"indexB":"node"})
+        snapped_points = snapped_points[["node", "snap_distance", "geometry"]]
+        snapped_points.index.name = None      
+        
+    # snap to lines
+    if B['geometry'].geom_type.isin(['LineString', 'MultiLineString']).all():
+        closest = closest.rename(columns={"indexB":"link"})        
+        # position of nearest point from start of the line
+        pos = closest.geometry.project(gpd.GeoSeries(closest.points))        
+        # get new point location geometry
+        snapped_points = closest.geometry.interpolate(pos)
+        snapped_points = gpd.GeoDataFrame(data=closest ,geometry=snapped_points, crs=crs)
+        # determine whether the snapped point is closer to the start or end node
+        snapped_points["line_position"] = closest.geometry.project(snapped_points, normalized=True)
+        snapped_points.loc[snapped_points["line_position"]<0.5, "node"] = closest["start_node_name"]
+        snapped_points.loc[snapped_points["line_position"]>=0.5, "node"] = closest["end_node_name"]
+        snapped_points = snapped_points[["link", "node", "snap_distance", "line_position", "geometry"]]
+        snapped_points.index.name = None
+        
+    return snapped_points
+
+def _backgound(A, B):
+    
+    hull_geom = A.unary_union.convex_hull
+    hull_data = gpd.GeoDataFrame(pd.DataFrame([{'geometry': hull_geom}]), crs=A.crs)
+    
+    background_geom = hull_data.overlay(B, how='difference').unary_union
+   
+    background = gpd.GeoDataFrame(pd.DataFrame([{'geometry': background_geom}]), crs=A.crs)
+    background.index = ['BACKGROUND']
+    
+    return background
+
+
+def intersect(A, B, B_value=None, include_background=False, background_value=0):
+    """
+    Intersect Points, Lines or Polygons in A with Points, Lines, or Polygons in B.
+    Return statistics on the intersection.
+    
+    The function returns information about the intersection for each geometry 
+    in A. Each geometry in B is assigned a value based on a column of data in 
+    that GeoDataFrame.  Note the CRS of A must equal the CRS of B.
+    
+    Parameters
+    ----------
+    A : geopandas GeoDataFrame
+        GeoDataFrame containing Point, LineString, or Polygon geometries
+    B : geopandas GeoDataFrame
+        GeoDataFrame containing  Point, LineString, or Polygon geometries
+    B_value : str or None (optional)
+        Column name in B used to assign a value to each geometry.
+        Default is None.
+    include_background : bool (optional) 
+         Include background, defined as space covered by A that is not covered by B 
+         (overlay difference between A and B). The background geometry is added
+         to B and is given the name 'BACKGROUND'. Default is False.
+    background_value : int or float (optional)
+        The value given to background space. This value is used in the intersection 
+        statistics if a B_value column name is provided. Default is 0.
+        
+    Returns
+    -------
+    pandas DataFrame
+        Intersection statistics (index = A.index, columns = defined below)
+        Columns include:
+            - n: number of intersecting B geometries
+            - intersections: list of intersecting B indices
+            
+        If B_value is given:
+            - values: list of intersecting B values
+            - sum: sum of the intersecting B values
+            - min: minimum of the intersecting B values
+            - max: maximum of the intersecting B values
+            - mean: mean of the intersecting B values
+            
+        If A contains Lines and B contains Polygons:
+            - weighted_mean: weighted mean of intersecting B values
+
+    """
+    if not has_shapely or not has_geopandas:
+        raise ModuleNotFoundError('shapley and geopandas are required')
+        
+    isinstance(A, gpd.GeoDataFrame)
+    assert (A['geometry'].geom_type).isin(['Point', 'LineString', 
+                                           'MultiLineString', 'Polygon', 
+                                           'MultiPolygon']).all()
+    isinstance(B, gpd.GeoDataFrame)
+    assert (B['geometry'].geom_type).isin(['Point', 'LineString', 
+                                           'MultiLineString', 'Polygon', 
+                                           'MultiPolygon']).all()
+    if isinstance(B_value, str):
+        assert B_value in B.columns
+    isinstance(include_background, bool)
+    isinstance(background_value, (int, float))
+    assert A.crs == B.crs
+    
+    if include_background:
+        background = _backgound(A, B)
+        if B_value is not None:
+            background[B_value] = background_value
+        B = B.append(background)
+        
+    intersects = gpd.sjoin(A, B, predicate='intersects')
+    intersects.index.name = 'name'
+    
+    n = intersects.groupby('name')['geometry'].count()
+    B_indices = intersects.groupby('name')['index_right'].apply(list)
+    stats = pd.DataFrame(index=A.index, data={'intersections': B_indices,
+                                              'n': n,})
+    stats['n'] = stats['n'].fillna(0)
+    stats['n'] = stats['n'].apply(int)
+    stats.loc[stats['intersections'].isnull(), 'intersections'] = stats.loc[stats['intersections'].isnull(), 'intersections'] .apply(lambda x: [])
+    
+    if B_value is not None:
+        stats['values'] = intersects.groupby('name')[B_value].apply(list)
+        stats['sum'] = intersects.groupby('name')[B_value].sum()
+        stats['min'] = intersects.groupby('name')[B_value].min()
+        stats['max'] = intersects.groupby('name')[B_value].max()
+        stats['mean'] = intersects.groupby('name')[B_value].mean()
+        
+        stats = stats.reindex(['intersections', 'values', 'n', 'sum', 'min', 'max', 'mean'], axis=1)
+        stats.loc[stats['values'].isnull(), 'values'] = stats.loc[stats['values'].isnull(), 'values'] .apply(lambda x: [])
+        
+    weighted_mean = False
+    if (A['geometry'].geom_type).isin(['LineString', 'MultiLineString']).all():
+        if (B['geometry'].geom_type).isin(['Polygon', 'MultiPolygon']).all():
+            weighted_mean = True
+            
+    if weighted_mean and B_value is not None:
+        stats['weighted_mean'] = 0
+        A_length = A.length
+        covered_length = pd.Series(0, index = A.index)
+        
+        for i in B.index:
+            B_geom = gpd.GeoDataFrame(B.loc[[i],:], crs=B.crs)
+            val = float(B_geom[B_value])
+            A_subset = A.loc[stats['intersections'].apply(lambda x: i in x),:]
+            #print(i, lines_subset)
+            A_clip = gpd.clip(A_subset, B_geom) 
+            A_clip_length = A_clip.length
+            A_clip_index = A_clip.index
+            
+            if A_clip_length.shape[0] > 0:
+                fraction_length = A_clip_length/A_length[A_clip_index]
+                covered_length[A_clip_index] = covered_length[A_clip_index] + fraction_length
+                weighed_val = fraction_length*val
+                stats.loc[A_clip_index, 'weighted_mean'] = stats.loc[A_clip_index, 'weighted_mean'] + weighed_val
+        
+        # Normalize weighted mean by covered length (can be over 1 if polygons overlap)
+        # Can be less than 1 if there are gaps (when background is not used)
+        stats['weighted_mean'] = stats['weighted_mean']/covered_length
+        
+        # Covered_length is NaN if length A is 0, set weighted mean to mean
+        stats.loc[covered_length.isna(), 'weighted_mean'] = stats.loc[covered_length.isna(), 'mean']
+        
+        # No intersection, set weighted mean to NaN
+        stats.loc[stats['n']==0, 'weighted_mean'] = np.NaN
+        
+    stats.index.name = None
+    
+    return stats
diff --git a/wntr/gis/network.py b/wntr/gis/network.py
new file mode 100644
index 000000000..e6bd19c80
--- /dev/null
+++ b/wntr/gis/network.py
@@ -0,0 +1,309 @@
+"""
+The wntr.gis.network module contains methods to convert between water network 
+models and GIS formatted data
+"""
+
+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
+
+
+class WaterNetworkGIS:
+    """
+    Water network GIS class 
+    
+    Contains methods to create GeoDataFrames from WaterNetworkModel and 
+    create WaterNetworkModel from GeoDataFrames.
+    
+    Parameters
+    ----------
+    gis_data : dict, optional
+        Dictionary of GeoDataFrames containing data to populate an instance 
+        of WaterNetworkGIS.  Valid dictionary keys are 'junction', 'tanks',
+        'reservoirs', 'pipes', 'pumps', and 'valves'
+    
+    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.tanks = gpd.GeoDataFrame()
+        self.reservoirs = gpd.GeoDataFrame()
+        self.pipes = gpd.GeoDataFrame()
+        self.pumps = gpd.GeoDataFrame()
+        self.valves = gpd.GeoDataFrame()
+        
+        if isinstance(gis_data, dict):
+            if 'junctions' in gis_data.keys():
+                assert isinstance(gis_data['junctions'], gpd.GeoDataFrame)
+                self.junctions = gis_data['junctions']
+                    
+            if 'tanks' in gis_data.keys():
+                assert isinstance(gis_data['tanks'], gpd.GeoDataFrame)
+                self.tanks = gis_data['tanks']
+                    
+            if 'reservoirs' in gis_data.keys():
+                assert isinstance(gis_data['reservoirs'], gpd.GeoDataFrame)
+                self.reservoirs = gis_data['reservoirs']
+                    
+            if 'pipes' in gis_data.keys():
+                assert isinstance(gis_data['pipes'], gpd.GeoDataFrame)
+                self.pipes = gis_data['pipes']
+                    
+            if 'pumps' in gis_data.keys():
+                assert isinstance(gis_data['pumps'], gpd.GeoDataFrame)
+                self.pumps = gis_data['pumps']
+                    
+            if 'valves' in gis_data.keys():
+                assert isinstance(gis_data['valves'], gpd.GeoDataFrame)
+                self.valves = gis_data['valves']
+
+        
+    def create_gis(self, wn, crs: str = None, pumps_as_points: bool = False, 
+                   valves_as_points: bool = False,) -> None:
+        """
+        Create GIS data from a water network model.
+        
+        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
+        pumps_as_points : bool, optional
+            Represent pumps as points (True) or lines (False), by default False
+        valves_as_points : bool, optional
+            Represent valves as points (True) or lines (False), by default False
+        """
+        
+        def _extract_geodataframe(df, crs=None, links_as_points=False):
+            # Drop any column with all NaN
+            df = df.loc[:, ~df.isna().all()]
+            
+            if df.shape[0] > 0:
+                # Define geom
+                if 'node_type' in df.columns:
+                    geom = [Point((x,y)) for x,y in df['coordinates']]
+                elif 'link_type' in df.columns:
+                    geom = []
+                    for link_name in df['name']:
+                        link = wn.get_link(link_name)
+                        if links_as_points: #Point
+                            geom.append(Point(link.start_node.coordinates))
+                        else: # LineString
+                            ls = list()
+                            ls.append(link.start_node.coordinates)
+                            for v in link.vertices:
+                                ls.append(v)
+                            ls.append(link.end_node.coordinates)
+                            geom.append(LineString(ls))
+                    
+                # Drop column if not a str, float, int, or bool
+                # This could be extended to keep additional data type (list, 
+                # tuple, network elements like Patterns, Curves)
+                drop_cols = []
+                for col in df.columns:
+                    if not isinstance(df.iloc[0][col], (str, float, int, bool)):
+                        drop_cols.append(col) 
+                df = df.drop(columns=drop_cols)
+                
+                # Set index
+                if len(df) > 0:
+                    df.set_index('name', inplace=True)
+                    df.index.name = None
+                
+                df = gpd.GeoDataFrame(df, crs=crs, geometry=geom)
+            else:
+                df = gpd.GeoDataFrame()
+                
+            return df
+        
+        # Convert the WaterNetworkModel to a dictionary
+        wn_dict = wn.to_dict()
+        # Create dataframes for node and link attributes
+        df_nodes = pd.DataFrame(wn_dict['nodes'])
+        df_links = pd.DataFrame(wn_dict['links'])
+        
+        # Junctions
+        df = df_nodes[df_nodes['node_type'] == 'Junction']
+        self.junctions = _extract_geodataframe(df, crs)
+        
+        # Tanks
+        df = df_nodes[df_nodes['node_type'] == 'Tank']
+        self.tanks = _extract_geodataframe(df, crs)
+            
+        # Reservoirs
+        df = df_nodes[df_nodes['node_type'] == 'Reservoir']
+        self.reservoirs = _extract_geodataframe(df, crs)
+            
+        # Pipes
+        df = df_links[df_links['link_type'] == 'Pipe']
+        self.pipes = _extract_geodataframe(df, crs, False)
+            
+        # Pumps
+        df = df_links[df_links['link_type'] == 'Pump']
+        self.pumps = _extract_geodataframe(df, crs, pumps_as_points)
+            
+        # Valves
+        df = df_links[df_links['link_type'] == 'Valve']
+        self.valves = _extract_geodataframe(df, crs, valves_as_points) 
+        
+    def create_wn(self, append=None):
+        """
+        Create or append a WaterNetworkModel from GeoDataFrames
+        
+       Parameters
+       ----------
+        append : WaterNetworkModel or None, optional
+            Existing WaterNetworkModel to append.  If None, a new WaterNetworkModel 
+            is created.
+        """
+        # Convert the WaterNetworkGIS to a dictionary
+        wn_dict = {}
+        wn_dict['nodes'] = []
+        wn_dict['links'] = []
+
+        for element in [self.junctions, self.tanks, self.reservoirs]:
+            if element.shape[0] > 0:
+                assert (element['geometry'].geom_type).isin(['Point']).all()
+                df = element.reset_index()
+                df.rename(columns={'index':'name', 'geometry':'coordinates'}, inplace=True)
+                df['coordinates'] = [[x,y] for x,y in zip(df['coordinates'].x, 
+                                                          df['coordinates'].y)]
+                wn_dict['nodes'].extend(df.to_dict('records'))
+
+        for element in [self.pipes, self.pumps, self.valves]:
+            if element.shape[0] > 0:
+                assert 'start_node_name' in element.columns
+                assert 'end_node_name' in element.columns
+                df = element.reset_index()
+                df.rename(columns={'index':'name'}, inplace=True)
+                # TODO: create vertices from LineString geometry
+                df.drop(columns=['geometry'], inplace=True)
+                wn_dict['links'].extend(df.to_dict('records'))
+        
+        # Create WaterNetworkModel from dictionary
+        from wntr.network import from_dict
+        wn = from_dict(wn_dict, append)
+        
+        return wn
+                
+    def add_node_attributes(self, values, name):
+        """
+        Add attribute to junctions, tanks, or reservoirs GeoDataFrames
+
+        Parameters
+        ----------
+        values : dict or Series or row of a DataFrame
+            Attribute values
+        name : str
+            Attribute name
+        """
+        for node_name, value in values.items():
+            if node_name in self.junctions.index:
+                if name not in self.junctions.columns:
+                    self.junctions[name] = np.nan
+                self.junctions.loc[node_name, name] = value
+            elif node_name in self.tanks.index:
+                if name not in self.tanks.columns:
+                    self.tanks[name] = np.nan
+                self.tanks.loc[node_name, name] = value
+            elif node_name in self.reservoirs.index:
+                if name not in self.reservoirs.columns:
+                    self.reservoirs[name] = np.nan
+                self.reservoirs.loc[node_name, name] = value
+            
+    def add_link_attributes(self, values, name):
+        """
+        Add attribute to pipes, pumps, or valves GeoDataFrames
+
+        Parameters
+        ----------
+        values : dict or Series or row of a DataFrame
+            Attribute values
+        name : str
+            Attribute name 
+        """
+        for link_name, value in values.items():
+            if link_name in self.pipes.index:
+                if name not in self.pipes.columns:
+                    self.pipes[name] = np.nan
+                self.pipes.loc[link_name, name] = value
+            elif link_name in self.valves.index:
+                if name not in self.valves.columns:
+                    self.valves[name] = np.nan
+                self.valves.loc[link_name, name] = value
+            elif link_name in self.pumps.index:
+                if name not in self.pumps.columns:
+                    self.pumps[name] = np.nan
+                self.pumps.loc[link_name, name] = value
+
+    def write(self, prefix: str, driver="GeoJSON") -> None:
+        """
+        Write the WaterNetworkGIS object to GIS file(s) with names 
+        constructed from parameters.
+
+        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()
+        
+        if len(self.junctions) > 0:
+            filename = prefix + "_junctions" + extension
+            self.junctions.to_file(filename, driver=driver)
+                
+        if len(self.tanks) > 0:
+            filename = prefix + "_tanks" + extension
+            self.tanks.to_file(filename, driver=driver)
+            
+        if len(self.reservoirs) > 0:
+            filename = prefix + "_reservoirs" + extension
+            self.reservoirs.to_file(filename, driver=driver)
+            
+        if len(self.pipes) > 0:
+            filename = prefix + "_pipes" + extension
+            self.pipes.to_file(filename, driver=driver)
+            
+        if len(self.pumps) > 0:
+            filename = prefix + "_pumps" + extension
+            self.pumps.to_file(filename, driver=driver)
+            
+        if len(self.valves) > 0:
+            filename = prefix + "_valves" + extension
+            self.valves.to_file(filename, driver=driver)
+    
\ No newline at end of file
diff --git a/wntr/metrics/topographic.py b/wntr/metrics/topographic.py
index e741c82cb..2f5006d9a 100644
--- a/wntr/metrics/topographic.py
+++ b/wntr/metrics/topographic.py
@@ -221,8 +221,8 @@ def valve_segments(G, valve_layer):
     matrix_link_names = ['L_'+k for u,v,k in links]
 
     # Pipe-node connectivity matrix
-    A = nx.incidence_matrix(uG).todense().T
-    AC = pd.DataFrame(A, columns=matrix_node_names, index=matrix_link_names, dtype=int)
+    AC = nx.incidence_matrix(uG).astype('int').todense().T
+    AC = pd.DataFrame(AC, columns=matrix_node_names, index=matrix_link_names, dtype=int)
     
     # Valve-node connectivity matrix
     VC = pd.DataFrame(0, columns=matrix_node_names, index=matrix_link_names)
@@ -231,7 +231,8 @@ def valve_segments(G, valve_layer):
     
     # Valve deficient matrix (anti-valve matrix)
     VD = AC - VC
-
+    del AC, VC
+    
     # Direct connectivity matrix
     NI = pd.DataFrame(np.identity(len(matrix_node_names)),
                       index = matrix_node_names, columns = matrix_node_names)
@@ -239,7 +240,9 @@ def valve_segments(G, valve_layer):
                       index = matrix_link_names, columns = matrix_link_names)
     DC_left = pd.concat([NI, VD], sort=False)
     DC_right = pd.concat([VD.T, LI], sort=False)
+    del LI, VD
     DC = pd.concat([DC_left, DC_right], axis=1, sort=False)
+    del DC_left, DC_right
     DC = DC.astype(int)
 
     # initialization for looping routine
diff --git a/wntr/network/__init__.py b/wntr/network/__init__.py
index fa10cd745..d6f2d48c3 100644
--- a/wntr/network/__init__.py
+++ b/wntr/network/__init__.py
@@ -1,6 +1,6 @@
 """
 The wntr.network package contains methods to define a water network model,
-network controls, and graph representation of the network.
+network controls, and water network model I/O.
 """
 from .base import Node, Link, NodeType, LinkType, LinkStatus
 from .elements import Junction, Reservoir, Tank, Pipe, Pump, Valve, Pattern, \
@@ -10,5 +10,10 @@
 from .options import Options
 from .controls import Comparison, ControlPriority, TimeOfDayCondition, \
     SimTimeCondition, ValueCondition, TankLevelCondition, RelativeCondition, \
-    OrCondition, AndCondition, ControlAction, Control, ControlChecker, ControlChangeTracker, Rule
-from .io import to_dict, from_dict, read_inpfile, write_inpfile, read_json, write_json
\ No newline at end of file
+    OrCondition, AndCondition, ControlAction, Control, ControlChecker, \
+    ControlChangeTracker, Rule
+from .io import to_dict, from_dict, to_gis, from_gis, \
+    read_inpfile, write_inpfile, \
+    read_json, write_json, \
+    read_geojson, write_geojson, \
+    read_shapefile, write_shapefile
diff --git a/wntr/network/elements.py b/wntr/network/elements.py
index 0a05a964f..16b53825d 100644
--- a/wntr/network/elements.py
+++ b/wntr/network/elements.py
@@ -562,9 +562,9 @@ def overflow(self):
     def overflow(self, value):
         if isinstance(value, six.string_types):
             value = value.upper()
-            if value in ["YES", "TRUE"]:
+            if value in ["YES", "TRUE", "1"]:
                 value = True
-            elif value in ["NO", "FALSE"]:
+            elif value in ["NO", "FALSE", "0"]:
                 value = False
             else:
                 raise ValueError('The overflow entry must "YES" or "NO"')
@@ -913,7 +913,7 @@ def cv(self):
     @cv.setter
     def cv(self, value): 
         self._cv = value
-
+        
     @property
     def bulk_coeff(self):
         """float or None : if not None, then a pipe specific bulk reaction coefficient"""
@@ -1218,7 +1218,7 @@ def _compare(self, other):
         if not super(HeadPump, self)._compare(other):
             return False
         if self.pump_type == other.pump_type and \
-           self.get_pump_curve() == other.get_pump_curve():
+           self.pump_curve_name == other.pump_curve_name:
             return True
         return False
     
@@ -1429,7 +1429,7 @@ def _compare(self, other):
         if not super(PowerPump, self)._compare(other):
             return False
         if self.pump_type == other.pump_type and \
-           self.power == other.power:
+            abs(self.power - other.power)<1e-9:
             return True
         return False
     
diff --git a/wntr/network/io.py b/wntr/network/io.py
index 4b8605cd3..babd4bb80 100644
--- a/wntr/network/io.py
+++ b/wntr/network/io.py
@@ -1,9 +1,9 @@
 # coding: utf-8
 
 """
-The wntr.network.elements module includes elements of a water network model, 
-including junction, tank, reservoir, pipe, pump, valve, pattern, timeseries, 
-demands, curves, and sources.
+The wntr.network.io module includes functions that convert the water network
+model to other data formats, create a water network model from file, and write 
+the water network model to a file.
 
 .. rubric:: Contents
 
@@ -11,10 +11,16 @@
 
     to_dict
     from_dict
+    to_gis
+    from_gis
     write_json
     read_json
     write_inpfile
     read_inpfile
+    write_geojson
+    read_geojson
+    write_shapefile
+    read_shapefile
 
 """
 import logging
@@ -23,22 +29,31 @@
 import wntr.epanet
 from wntr.epanet.util import FlowUnits
 import wntr.network.model
-
+from wntr.gis.network import WaterNetworkGIS
+try:
+    import geopandas as gpd
+    has_geopandas = True
+except ModuleNotFoundError:
+    gpd = None
+    has_geopandas = False
+    
 logger = logging.getLogger(__name__)
 
 
 def to_dict(wn) -> dict:
-    """Dictionary representation of the WaterNetworkModel.
+    """
+    Convert a WaterNetworkModel into a dictionary
     
     Parameters
     ----------
     wn : WaterNetworkModel
-        the water network model to convert to a dictionary
+        Water network model
 
     Returns
     -------
     dict
-        dictionary representation of the water network model
+        Dictionary representation of the WaterNetworkModel
+        
     """
     from wntr import __version__
 
@@ -63,23 +78,22 @@ def to_dict(wn) -> dict:
     return d
 
 
-def from_dict(
-    d: dict, append=None,
-):
+def from_dict(d: dict, append=None):
     """
-    Create or append a water network model object from a dictionary.
+    Create or append a WaterNetworkModel from a dictionary
 
     Parameters
     ----------
     d : dict
-        dictionary representation of the water network model
-    append : WaterNetworkModel, optional
-        an existing water network model to append results to
+        Dictionary representation of the water network model
+    append : WaterNetworkModel or None, optional
+        Existing WaterNetworkModel to append.  If None, a new WaterNetworkModel 
+        is created.
 
     Returns
     -------
     WaterNetworkModel
-        a new water network model or the appended existing model
+    
     """
     from wntr.epanet.io import _read_control_line, _EpanetRule
 
@@ -108,9 +122,7 @@ def from_dict(
         wn.options.__init__(**d["options"])
     if "curves" in d:
         for curve in d["curves"]:
-            wn.add_curve(
-                name=curve["name"], curve_type=curve["curve_type"], xy_tuples_list=curve["points"]
-            )
+            wn.add_curve(name=curve["name"], curve_type=curve["curve_type"], xy_tuples_list=curve["points"])
     if "patterns" in d:
         for pattern in d["patterns"]:
             wn.add_pattern(name=pattern["name"], pattern=pattern["multipliers"])
@@ -142,6 +154,9 @@ def from_dict(
                 j.pressure_exponent = node.setdefault("pressure_exponent")
                 j.required_pressure = node.setdefault("required_pressure")
                 j.tag = node.setdefault("tag")
+                # custom additional attributes
+                for attr in list(set(node.keys()) - set(dir(j))):
+                    setattr( j, attr, node[attr] )
                 if dl is not None and len(dl) > 1:
                     for i in range(1, len(dl)):
                         base_val = dl[i].setdefault("base_val", 0.0)
@@ -170,6 +185,9 @@ def from_dict(
                     t.mixing_model = node.setdefault("mixing_model")
                 t.tag = node.setdefault("tag")
                 t.bulk_coeff = node.setdefault("bulk_coeff")
+                # custom additional attributes
+                for attr in list(set(node.keys()) - set(dir(t))):
+                    setattr( t, attr, node[attr] )
             elif node["node_type"] == "Reservoir":
                 wn.add_reservoir(
                     name,
@@ -180,6 +198,9 @@ def from_dict(
                 r = wn.get_node(name)
                 r.initial_quality = node.setdefault("initial_quality", 0.0)
                 r.tag = node.setdefault("tag")
+                # custom additional attributes
+                for attr in list(set(node.keys()) - set(dir(r))):
+                    setattr( r, attr, node[attr] )
             else:
                 raise ValueError("Illegal node type '{}'".format(node["node_type"]))
     if "links" in d:
@@ -202,6 +223,9 @@ def from_dict(
                 p.tag = link.setdefault("tag")
                 p.vertices = link.setdefault("vertices", list())
                 p.wall_coeff = link.setdefault("wall_coeff")
+                # custom additional attributes
+                for attr in list(set(link.keys()) - set(dir(p))):
+                    setattr( p, attr, link[attr] )
             elif link["link_type"] == "Pump":
                 pump_type = link.setdefault("pump_type", "POWER")
                 wn.add_pump(
@@ -223,6 +247,9 @@ def from_dict(
                 p.initial_setting = link.setdefault("initial_setting")
                 p.tag = link.setdefault("tag")
                 p.vertices = link.setdefault("vertices", list())
+                # custom additional attributes
+                for attr in list(set(link.keys()) - set(dir(p))):
+                    setattr( p, attr, link[attr] )
             elif link["link_type"] == "Valve":
                 valve_type = link["valve_type"]
                 wn.add_valve(
@@ -238,6 +265,9 @@ def from_dict(
                 v = wn.get_link(name)
                 if valve_type.lower() == "gpv":
                     v.headloss_curve_name = link.setdefault("headloss_curve_name")
+                # custom additional attributes
+                for attr in list(set(link.keys()) - set(dir(v))):
+                    setattr( v, attr, link[attr] )
             else:
                 raise ValueError("Illegal link type '{}'".format(link["link_type"]))
     if "sources" in d:
@@ -286,6 +316,57 @@ def from_dict(
                 raise ValueError("Illegal control type '{}'".format(ctrl_type))
     return wn
 
+def to_gis(wn, crs=None, pumps_as_points=False, valves_as_points=False):
+    """
+    Convert a WaterNetworkModel into GeoDataFrames
+    
+    Parameters
+    ----------
+    wn : WaterNetworkModel
+        Water network model
+    crs : str, optional
+        Coordinate reference system, by default None
+    pumps_as_points : bool, optional
+        Represent pumps as points (True) or lines (False), by default False
+    valves_as_points : bool, optional
+        Represent valves as points (True) or lines (False), by default False
+        
+    Returns
+    -------
+    WaterNetworkGIS object that contains GeoDataFrames
+        
+    """
+    gis_data = WaterNetworkGIS()
+    gis_data.create_gis(wn, crs, pumps_as_points, valves_as_points)
+    
+    return gis_data
+
+def from_gis(gis_data, append=None):
+    """
+    Create or append a WaterNetworkModel from GeoDataFrames
+    
+    Parameters
+    ----------
+    gis_data : WaterNetworkGIS or dictionary of GeoDataFrames
+        GeoDataFrames containing water network attributes. If gis_data is a 
+        dictionary, then the keys are junctions, tanks, reservoirs, pipes, 
+        pumps, and valves. If the pumps or valves are Points, they will be 
+        converted to Lines with the same start and end node location.
+    append : WaterNetworkModel or None, optional
+        Existing WaterNetworkModel to append.  If None, a new WaterNetworkModel 
+        is created.
+        
+    Returns
+    -------
+    WaterNetworkModel
+        
+    """
+    if isinstance(gis_data, dict):
+        gis_data = WaterNetworkGIS(gis_data)
+
+    wn = gis_data.create_wn(append=append)
+    
+    return wn
 
 def write_json(
     wn, path_or_buf, **kw_json,
@@ -298,7 +379,8 @@ def write_json(
     path_or_buf : str or IO stream
         Name of the file or file pointer
     kw_json : keyword arguments
-        arguments to pass directly to `json.dump`
+        Arguments to pass directly to `json.dump`
+        
     """
     if isinstance(path_or_buf, str):
         with open(path_or_buf, "w") as fout:
@@ -307,22 +389,24 @@ def write_json(
         json.dump(to_dict(wn), path_or_buf, **kw_json)
 
 
-def read_json(
-    path_or_buf, append=None, **kw_json,
-):
+def read_json(path_or_buf, append=None, **kw_json):
     """
-    Create a WaterNetworkModel from a JSON file.
+    Create or append a WaterNetworkModel from a JSON file
 
     Parameters
     ----------
     f : str
         Name of the file or file pointer
+    append : WaterNetworkModel or None, optional
+        Existing WaterNetworkModel to append.  If None, a new WaterNetworkModel 
+        is created.
     kw_json : keyword arguments
-        keyword arguments to pass to `json.load`
+        Keyword arguments to pass to `json.load`
 
     Returns
     -------
     WaterNetworkModel
+    
     """
     if isinstance(path_or_buf, str):
         with open(path_or_buf, "r") as fin:
@@ -336,7 +420,7 @@ def write_inpfile(
     wn, filename: str, units=None, version: float = 2.2, force_coordinates: bool = False,
 ):
     """
-    Writes the current water network model to an EPANET INP file
+    Write the WaterNetworkModel to an EPANET INP file
 
     .. note::
 
@@ -346,21 +430,16 @@ def write_inpfile(
         the user to force EPANET 2.0 compatibility at the expense of pressured-dependent analysis 
         options being turned off.
 
-
     Parameters
     ----------
-    wn : WaterNetworkModel
-        model object to write INP file for
-
+    wn : wntr WaterNetworkModel
+        Water network model
     filename : string
         Name of the inp file.
-
     units : str, int or FlowUnits
         Name of the units being written to the inp file.
-
     version : float, {2.0, **2.2**}
         Optionally specify forcing EPANET 2.0 compatibility.
-
     force_coordinates : bool
         This only applies if `self.options.graphics.map_filename` is not `None`,
         and will force the COORDINATES section to be written even if a MAP file is
@@ -373,24 +452,218 @@ def write_inpfile(
         wn._inpfile = wntr.epanet.InpFile()
     if units is None:
         units = wn._options.hydraulic.inpfile_units
-    wn._inpfile.write(
-        filename, wn, units=units, version=version, force_coordinates=force_coordinates
-    )
+    wn._inpfile.write(filename, wn, units=units, version=version, force_coordinates=force_coordinates)
 
 
-def read_inpfile(
-    filename, append=None,
-):
+def read_inpfile(filename, append=None):
     """
-    Defines water network model components from an EPANET INP file
+    Create or append a WaterNetworkModel from an EPANET INP file
 
     Parameters
     ----------
     filename : string
         Name of the INP file.
-
+    append : WaterNetworkModel or None, optional
+        Existing WaterNetworkModel to append.  If None, a new WaterNetworkModel 
+        is created.
+        
+    Returns
+    -------
+    WaterNetworkModel
+    
     """
     inpfile = wntr.epanet.InpFile()
     wn = inpfile.read(filename, wn=append)
     wn._inpfile = inpfile
+    
+    return wn
+
+def _write_geopandas(wn, prefix, crs=None, pumps_as_points=True,
+                     valves_as_points=True, driver="GeoJSON"):
+    wn_gis = wn.to_gis(crs, pumps_as_points=pumps_as_points, 
+                       valves_as_points=valves_as_points)
+    wn_gis.write(prefix=prefix, driver=driver)
+
+def _read_geopandas(files, index_col='index', append=None, field_name_map=None):
+    
+    gis_data = WaterNetworkGIS()
+    
+    if 'junctions' in files.keys():
+        gis_data.junctions = gpd.read_file(files['junctions']).set_index(index_col)
+        if field_name_map is not None:
+            gis_data.junctions.rename(columns=field_name_map['junctions'], inplace=True)
+            
+    if 'tanks' in files.keys():
+        gis_data.tanks = gpd.read_file(files['tanks']).set_index(index_col)
+        if field_name_map is not None:
+            gis_data.tanks.rename(columns=field_name_map['tanks'], inplace=True)
+            
+    if 'reservoirs' in files.keys():
+        gis_data.reservoirs = gpd.read_file(files['reservoirs']).set_index(index_col)
+        if field_name_map is not None:
+            gis_data.reservoirs.rename(columns=field_name_map['reservoirs'], inplace=True)
+            
+    if 'pipes' in files.keys():
+        gis_data.pipes = gpd.read_file(files['pipes']).set_index(index_col)
+        if field_name_map is not None:
+            gis_data.pipes.rename(columns=field_name_map['pipes'], inplace=True)
+    
+    if 'pumps' in files.keys():
+        gis_data.pumps = gpd.read_file(files['pumps']).set_index(index_col)
+        if field_name_map is not None:
+            gis_data.pumps.rename(columns=field_name_map['pumps'], inplace=True)
+    
+    if 'valves' in files.keys():
+        gis_data.valves = gpd.read_file(files['valves']).set_index(index_col)
+        if field_name_map is not None:
+            gis_data.valves.rename(columns=field_name_map['valves'], inplace=True)
+    
+    wn = gis_data.create_wn(append=append)
+        
     return wn
+
+def write_geojson(
+    wn,
+    prefix: str,
+    crs=None,
+    pumps_as_points=True,
+    valves_as_points=True,
+):
+    """
+    Write the WaterNetworkModel to a set of GeoJSON files, one file for each network element.
+
+    Note that this is a simple export of the water network data.
+    To add results of a simulation or analysis, do:
+
+    .. code::
+
+        wn_gis = wn.to_gis()
+        wn_gis.add_node_attributes(some_data_to_add, 'name_of_attribute')
+        wn_gis.write(...)
+
+    Parameters
+    ----------
+    wn : wntr WaterNetworkModel
+        Water network model
+    prefix : str
+        File prefix
+    crs : str, optional
+        Coordinate reference system, by default None
+    pumps_as_points : bool, optional
+        Represent pumps as points (True) or lines (False), by default False
+    valves_as_points : bool, optional
+        Represent valves as points (True) or lines (False), by default False
+        
+    """
+    _write_geopandas(wn, prefix, crs, pumps_as_points, valves_as_points, "GeoJSON")
+
+def read_geojson(files, index_col='index', append=None):
+    """
+    Create or append a WaterNetworkModel from GeoJSON files
+
+    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
+    append : WaterNetworkModel or None, optional
+        Existing WaterNetworkModel to append.  If None, a new WaterNetworkModel 
+        is created.
+
+    Returns
+    -------
+    WaterNetworkModel
+    
+    """
+    return _read_geopandas(files, index_col, append)
+      
+def write_shapefile(
+    wn,
+    prefix: str,
+    crs=None,
+    pumps_as_points=True,
+    valves_as_points=True
+):
+    """
+    Write the WaterNetworkModel to a set of ESRI Shapefiles, one directory for each network element.
+
+    Note that this is a simple export of the water network data.
+    To add results of a simulation or analysis, do:
+
+    .. code::
+
+        wn_gis = wn.to_gis()
+        wn_gis.add_node_attributes(some_data_to_add, 'name_of_attribute')
+        wn_gis.write(...)
+
+    Parameters
+    ----------
+    wn : wntr WaterNetworkModel
+        Water network model
+    prefix : str
+        File and directory prefix
+    crs : str, optional
+        Coordinate reference system, by default None
+    pumps_as_points : bool, optional
+        Represent pumps as points (True) or lines (False), by default False
+    valves_as_points : bool, optional
+        Represent valves as points (True) or lines (False), by default False
+        
+    """
+    _write_geopandas(wn, prefix, crs, pumps_as_points, valves_as_points, None)
+
+def read_shapefile(dirs, index_col='index', append=None):
+    """
+    
+    Create or append a WaterNetworkModel from ESRI Shapefiles
+
+    Parameters
+    ----------
+    dirs : dictionary
+        Dictionary of Shapefile directory names, where the keys are in the set 
+        ('junction', 'tanks', 'reservoirs', 'pipes', 'pumps', 'valves') and 
+        values are the corresponding Shapefile directories
+    index_col : str, optional
+        Column that contains the element name
+    append : WaterNetworkModel or None, optional
+        Existing WaterNetworkModel to append.  If None, a new WaterNetworkModel 
+        is created.
+
+    Returns
+    -------
+    WaterNetworkModel
+    
+    """
+    # ESRI Shapefiles truncate field names to 10 characters. The field_name_map
+    # maps truncated names to long names.  The following code assumes the 
+    # first 10 characters are unique.
+    element_attributes = {
+        'junctions': dir(wntr.network.elements.Junction),
+        'tanks': dir(wntr.network.elements.Tank),
+        'reservoirs': dir(wntr.network.elements.Reservoir),
+        'pipes': dir(wntr.network.elements.Pipe),
+        'pumps': dir(wntr.network.elements.Pump) + 
+                 dir(wntr.network.elements.PowerPump) + 
+                 dir(wntr.network.elements.HeadPump),
+        'valves': dir(wntr.network.elements.Valve) + 
+                  dir(wntr.network.elements.PRValve) + 
+                  dir(wntr.network.elements.PSValve) + 
+                  dir(wntr.network.elements.PBValve) + 
+                  dir(wntr.network.elements.FCValve) + 
+                  dir(wntr.network.elements.TCValve) + 
+                  dir(wntr.network.elements.GPValve)}
+    
+    field_name_map = {}
+    for element,attribute in element_attributes.items():
+        field_name_map[element] = {}
+        for field_name in attribute:
+            if (len(field_name) > 10) and (not field_name.startswith('_')):
+                field_name_map[element][field_name[0:10]] = field_name
+    
+    # TODO: pipe property is cv instead of check_valve, this should be updated
+    field_name_map['pipes']['check_valv'] = 'check_valve'
+    
+    return _read_geopandas(dirs, index_col, append, field_name_map)
diff --git a/wntr/network/model.py b/wntr/network/model.py
index 8d232c1be..cbde0bc4a 100644
--- a/wntr/network/model.py
+++ b/wntr/network/model.py
@@ -27,11 +27,30 @@
 
 from .base import AbstractModel, Link, LinkStatus, Registry
 from .controls import Control, Rule
-from .elements import (Curve, Demands, FCValve, GPValve, HeadPump, Junction,
-                       Pattern, PBValve, Pipe, PowerPump, PRValve, PSValve,
-                       Pump, Reservoir, Source, Tank, TCValve, TimeSeries,
-                       Valve)
+from .elements import (
+    Curve,
+    Demands,
+    FCValve,
+    GPValve,
+    HeadPump,
+    Junction,
+    Pattern,
+    PBValve,
+    Pipe,
+    PowerPump,
+    PRValve,
+    PSValve,
+    Pump,
+    Reservoir,
+    Source,
+    Tank,
+    TCValve,
+    TimeSeries,
+    Valve,
+)
+
 from .options import Options
+from .io import read_inpfile
 
 logger = logging.getLogger(__name__)
 
@@ -72,14 +91,14 @@ def __init__(self, inp_file_name=None):
 
         self._inpfile = None
         if inp_file_name:
-            self.read_inpfile(inp_file_name)
-            
+            read_inpfile(inp_file_name, append=self)
+
         # To be deleted and/or renamed and/or moved
         # Time parameters
         self.sim_time = 0.0
         self._prev_sim_time = None  # the last time at which results were accepted
-    
-    def _compare(self, other):
+
+    def _compare(self, other, level=1):
         """
         Parameters
         ----------
@@ -89,12 +108,14 @@ def _compare(self, other):
         -------
         bool
         """
-        if self.num_junctions  != other.num_junctions  or \
-           self.num_reservoirs != other.num_reservoirs or \
-           self.num_tanks      != other.num_tanks      or \
-           self.num_pipes      != other.num_pipes      or \
-           self.num_pumps      != other.num_pumps      or \
-           self.num_valves     != other.num_valves:
+        if (
+            self.num_junctions != other.num_junctions
+            or self.num_reservoirs != other.num_reservoirs
+            or self.num_tanks != other.num_tanks
+            or self.num_pipes != other.num_pipes
+            or self.num_pumps != other.num_pumps
+            or self.num_valves != other.num_valves
+        ):
             return False
         for name, node in self.nodes():
             if not node._compare(other.get_node(name)):
@@ -102,22 +123,24 @@ def _compare(self, other):
         for name, link in self.links():
             if not link._compare(other.get_link(name)):
                 return False
-        for name, pat in self.patterns():
-            if pat != other.get_pattern(name):
-                return False
-        for name, curve in self.curves():
-            if curve != other.get_curve(name):
-                return False
-        for name, source in self.sources():
-            if source != other.get_source(name):
-                return False
-        if self.options != other.options:
-            return False
-        for name, control in self.controls():
-            if not control._compare(other.get_control(name)):
+            
+        if level > 0:
+            for name, pat in self.patterns():
+                if pat != other.get_pattern(name):
+                    return False
+            for name, curve in self.curves():
+                if curve != other.get_curve(name):
+                    return False
+            for name, source in self.sources():
+                if source != other.get_source(name):
+                    return False
+            if self.options != other.options:
                 return False
+            for name, control in self.controls():
+                if not control._compare(other.get_control(name)):
+                    return False
         return True
-    
+
     @property
     def _shifted_time(self):
         """
@@ -143,17 +166,17 @@ def _clock_time(self):
         """
         Return the current time of day in seconds from 12 AM
         """
-        return self._shifted_time % (24*3600)
+        return self._shifted_time % (24 * 3600)
 
     @property
     def _clock_day(self):
         """Return the clock-time day of the simulation"""
         return int(self._shifted_time / 86400)
 
-    ### # 
+    ### #
     ### Iteratable attributes
     @property
-    def options(self): 
+    def options(self):
         """The model's options object
         
         Returns
@@ -162,9 +185,9 @@ def options(self):
         
         """
         return self._options
-    
+
     @property
-    def nodes(self): 
+    def nodes(self):
         """The node registry (as property) or a generator for iteration (as function call)
         
         Returns
@@ -173,9 +196,9 @@ def nodes(self):
         
         """
         return self._node_reg
-    
+
     @property
-    def links(self): 
+    def links(self):
         """The link registry (as property) or a generator for iteration (as function call)
         
         Returns
@@ -184,9 +207,9 @@ def links(self):
         
         """
         return self._link_reg
-    
+
     @property
-    def patterns(self): 
+    def patterns(self):
         """The pattern registry (as property) or a generator for iteration (as function call)
 
         Returns
@@ -195,9 +218,9 @@ def patterns(self):
 
         """
         return self._pattern_reg
-    
+
     @property
-    def curves(self): 
+    def curves(self):
         """The curve registry (as property) or a generator for iteration (as function call)
         
         Returns
@@ -206,7 +229,7 @@ def curves(self):
         
         """
         return self._curve_reg
-    
+
     def sources(self):
         """Returns a generator to iterate over all sources
 
@@ -216,7 +239,7 @@ def sources(self):
         """
         for source_name, source in self._sources.items():
             yield source_name, source
-        
+
     def controls(self):
         """Returns a generator to iterate over all controls
 
@@ -226,37 +249,37 @@ def controls(self):
         """
         for control_name, control in self._controls.items():
             yield control_name, control
-                
-    ### # 
+
+    ### #
     ### Element iterators
     @property
-    def junctions(self): 
+    def junctions(self):
         """Iterator over all junctions"""
         return self._node_reg.junctions
-    
+
     @property
-    def tanks(self): 
+    def tanks(self):
         """Iterator over all tanks"""
         return self._node_reg.tanks
-    
+
     @property
-    def reservoirs(self): 
-        """Iterator over all reservoirs"""        
+    def reservoirs(self):
+        """Iterator over all reservoirs"""
         return self._node_reg.reservoirs
-    
+
     @property
-    def pipes(self): 
-        """Iterator over all pipes"""        
+    def pipes(self):
+        """Iterator over all pipes"""
         return self._link_reg.pipes
-    
+
     @property
-    def pumps(self): 
-        """Iterator over all pumps"""        
+    def pumps(self):
+        """Iterator over all pumps"""
         return self._link_reg.pumps
-    
+
     @property
-    def valves(self): 
-        """Iterator over all valves"""        
+    def valves(self):
+        """Iterator over all valves"""
         return self._link_reg.valves
 
     @property
@@ -271,34 +294,34 @@ def power_pumps(self):
 
     @property
     def prvs(self):
-        """Iterator over all pressure reducing valves (PRVs)"""        
+        """Iterator over all pressure reducing valves (PRVs)"""
         return self._link_reg.prvs
 
     @property
     def psvs(self):
-        """Iterator over all pressure sustaining valves (PSVs)"""        
+        """Iterator over all pressure sustaining valves (PSVs)"""
         return self._link_reg.psvs
 
     @property
     def pbvs(self):
-        """Iterator over all pressure breaker valves (PBVs)"""        
+        """Iterator over all pressure breaker valves (PBVs)"""
         return self._link_reg.pbvs
 
     @property
     def tcvs(self):
-        """Iterator over all throttle control valves (TCVs)"""        
+        """Iterator over all throttle control valves (TCVs)"""
         return self._link_reg.tcvs
 
     @property
     def fcvs(self):
-        """Iterator over all flow control valves (FCVs)"""        
+        """Iterator over all flow control valves (FCVs)"""
         return self._link_reg.fcvs
 
     @property
     def gpvs(self):
-        """Iterator over all general purpose valves (GPVs)"""        
+        """Iterator over all general purpose valves (GPVs)"""
         return self._link_reg.gpvs
-    
+
     """
     ### # 
     ### Create blank, unregistered objects (for direct assignment)
@@ -311,11 +334,12 @@ def new_timeseries(self):
     def new_pattern(self):
         return Pattern(None, time_options=self._options.time)
     """
-    
-    ### # 
+
+    ### #
     ### Add elements to the model
-    def add_junction(self, name, base_demand=0.0, demand_pattern=None, 
-                     elevation=0.0, coordinates=None, demand_category=None):
+    def add_junction(
+        self, name, base_demand=0.0, demand_pattern=None, elevation=0.0, coordinates=None, demand_category=None
+    ):
         """
         Adds a junction to the water network model
 
@@ -334,12 +358,21 @@ def add_junction(self, name, base_demand=0.0, demand_pattern=None,
         demand_category  : string
             Name of the demand category
         """
-        self._node_reg.add_junction(name, base_demand, demand_pattern, 
-                                    elevation, coordinates, demand_category)
-
-    def add_tank(self, name, elevation=0.0, init_level=3.048,
-                 min_level=0.0, max_level=6.096, diameter=15.24,
-                 min_vol=0.0, vol_curve=None, overflow=False, coordinates=None):
+        self._node_reg.add_junction(name, base_demand, demand_pattern, elevation, coordinates, demand_category)
+
+    def add_tank(
+        self,
+        name,
+        elevation=0.0,
+        init_level=3.048,
+        min_level=0.0,
+        max_level=6.096,
+        diameter=15.24,
+        min_vol=0.0,
+        vol_curve=None,
+        overflow=False,
+        coordinates=None,
+    ):
         """
         Adds a tank to the water network model
 
@@ -367,9 +400,9 @@ def add_tank(self, name, elevation=0.0, init_level=3.048,
             X-Y coordinates of the node location.
             
         """
-        self._node_reg.add_tank(name, elevation, init_level, min_level, 
-                                max_level, diameter, min_vol, vol_curve, 
-                                overflow, coordinates)
+        self._node_reg.add_tank(
+            name, elevation, init_level, min_level, max_level, diameter, min_vol, vol_curve, overflow, coordinates
+        )
 
     def add_reservoir(self, name, base_head=0.0, head_pattern=None, coordinates=None):
         """
@@ -389,9 +422,18 @@ def add_reservoir(self, name, base_head=0.0, head_pattern=None, coordinates=None
         """
         self._node_reg.add_reservoir(name, base_head, head_pattern, coordinates)
 
-    def add_pipe(self, name, start_node_name, end_node_name, length=304.8,
-                 diameter=0.3048, roughness=100, minor_loss=0.0, initial_status='OPEN', 
-                 check_valve=False):
+    def add_pipe(
+        self,
+        name,
+        start_node_name,
+        end_node_name,
+        length=304.8,
+        diameter=0.3048,
+        roughness=100,
+        minor_loss=0.0,
+        initial_status="OPEN",
+        check_valve=False,
+    ):
         """
         Adds a pipe to the water network model
 
@@ -418,12 +460,21 @@ def add_pipe(self, name, start_node_name, end_node_name, length=304.8,
             False if the pipe does not have a check valve.
         
         """
-        self._link_reg.add_pipe(name, start_node_name, end_node_name, length, 
-                                diameter, roughness, minor_loss, initial_status, 
-                                check_valve)
+        self._link_reg.add_pipe(
+            name, start_node_name, end_node_name, length, diameter, roughness, minor_loss, initial_status, check_valve
+        )
 
-    def add_pump(self, name, start_node_name, end_node_name, pump_type='POWER',
-                 pump_parameter=50.0, speed=1.0, pattern=None, initial_status='OPEN'):
+    def add_pump(
+        self,
+        name,
+        start_node_name,
+        end_node_name,
+        pump_type="POWER",
+        pump_parameter=50.0,
+        speed=1.0,
+        pattern=None,
+        initial_status="OPEN",
+    ):
         """
         Adds a pump to the water network model
 
@@ -448,12 +499,21 @@ def add_pump(self, name, start_node_name, end_node_name, pump_type='POWER',
             Pump initial status. Options are 'OPEN' or 'CLOSED'.
         
         """
-        self._link_reg.add_pump(name, start_node_name, end_node_name, pump_type, 
-                                pump_parameter, speed, pattern, initial_status)
-    
-    def add_valve(self, name, start_node_name, end_node_name,
-                 diameter=0.3048, valve_type='PRV', minor_loss=0.0, 
-                 initial_setting=0.0, initial_status='ACTIVE'):
+        self._link_reg.add_pump(
+            name, start_node_name, end_node_name, pump_type, pump_parameter, speed, pattern, initial_status
+        )
+
+    def add_valve(
+        self,
+        name,
+        start_node_name,
+        end_node_name,
+        diameter=0.3048,
+        valve_type="PRV",
+        minor_loss=0.0,
+        initial_setting=0.0,
+        initial_status="ACTIVE",
+    ):
         """
         Adds a valve to the water network model
 
@@ -480,8 +540,9 @@ def add_valve(self, name, start_node_name, end_node_name,
         initial_status: string or LinkStatus
             Valve initial status. Options are 'OPEN',  'CLOSED', or 'ACTIVE'.
         """
-        self._link_reg.add_valve(name, start_node_name, end_node_name, diameter, 
-                                 valve_type, minor_loss, initial_setting, initial_status)
+        self._link_reg.add_valve(
+            name, start_node_name, end_node_name, diameter, valve_type, minor_loss, initial_setting, initial_status
+        )
 
     def add_pattern(self, name, pattern=None):
         """
@@ -522,7 +583,7 @@ def add_pattern(self, name, pattern=None):
         
         """
         self._pattern_reg.add_pattern(name, pattern)
-            
+
     def add_curve(self, name, curve_type, xy_tuples_list):
         """
         Adds a curve to the water network model
@@ -537,7 +598,7 @@ def add_curve(self, name, curve_type, xy_tuples_list):
             List of X-Y coordinate tuples on the curve.
         """
         self._curve_reg.add_curve(name, curve_type, xy_tuples_list)
-        
+
     def add_source(self, name, node_name, source_type, quality, pattern=None):
         """
         Adds a source to the water network model
@@ -564,8 +625,8 @@ def add_source(self, name, node_name, source_type, quality, pattern=None):
             pattern = self.get_pattern(pattern)
         source = Source(self, name, node_name, source_type, quality, pattern)
         self._sources[source.name] = source
-        self._pattern_reg.add_usage(source.strength_timeseries.pattern_name, (source.name, 'Source'))
-        self._node_reg.add_usage(source.node_name, (source.name, 'Source'))
+        self._pattern_reg.add_usage(source.strength_timeseries.pattern_name, (source.name, "Source"))
+        self._node_reg.add_usage(source.node_name, (source.name, "Source"))
 
     def add_control(self, name, control_object):
         """
@@ -579,28 +640,37 @@ def add_control(self, name, control_object):
             Control or Rule object.
         """
         if name in self._controls:
-            raise ValueError('The name provided for the control is already used. Please either remove the control with that name first or use a different name for this control.')
+            raise ValueError(
+                "The name provided for the control is already used. Please either remove the control with that name first or use a different name for this control."
+            )
         self._controls[name] = control_object
-    
-    
-    ### # 
+
+    ### #
     ### Remove elements from the model
     def remove_node(self, name, with_control=False, force=False):
         """Removes a node from the water network model"""
         node = self.get_node(name)
         if not force:
             if with_control:
-                x=[]
+                x = []
                 for control_name, control in self._controls.items():
                     if node in control.requires():
-                        logger.warning(control._control_type_str()+' '+control_name+' is being removed along with node '+name)
+                        logger.warning(
+                            control._control_type_str()
+                            + " "
+                            + control_name
+                            + " is being removed along with node "
+                            + name
+                        )
                         x.append(control_name)
                 for i in x:
                     self.remove_control(i)
             else:
                 for control_name, control in self._controls.items():
                     if node in control.requires():
-                        raise RuntimeError('Cannot remove node {0} without first removing control/rule {1}'.format(name, control_name))
+                        raise RuntimeError(
+                            "Cannot remove node {0} without first removing control/rule {1}".format(name, control_name)
+                        )
         self._node_reg.__delitem__(name)
 
     def remove_link(self, name, with_control=False, force=False):
@@ -608,27 +678,35 @@ def remove_link(self, name, with_control=False, force=False):
         link = self.get_link(name)
         if not force:
             if with_control:
-                x=[]
+                x = []
                 for control_name, control in self._controls.items():
                     if link in control.requires():
-                        logger.warning(control._control_type_str()+' '+control_name+' is being removed along with link '+name)
+                        logger.warning(
+                            control._control_type_str()
+                            + " "
+                            + control_name
+                            + " is being removed along with link "
+                            + name
+                        )
                         x.append(control_name)
                 for i in x:
                     self.remove_control(i)
             else:
                 for control_name, control in self._controls.items():
                     if link in control.requires():
-                        raise RuntimeError('Cannot remove link {0} without first removing control/rule {1}'.format(name, control_name))
+                        raise RuntimeError(
+                            "Cannot remove link {0} without first removing control/rule {1}".format(name, control_name)
+                        )
         self._link_reg.__delitem__(name)
 
-    def remove_pattern(self, name): 
+    def remove_pattern(self, name):
         """Removes a pattern from the water network model"""
         self._pattern_reg.__delitem__(name)
-        
-    def remove_curve(self, name): 
+
+    def remove_curve(self, name):
         """Removes a curve from the water network model"""
         self._curve_reg.__delitem__(name)
-        
+
     def remove_source(self, name):
         """Removes a source from the water network model
 
@@ -637,15 +715,17 @@ def remove_source(self, name):
         name : string
            The name of the source object to be removed
         """
-        logger.warning('You are deleting a source. This could have unintended \
+        logger.warning(
+            "You are deleting a source. This could have unintended \
             side effects. If you are replacing values, use get_source(name) \
-            and modify it instead.')
+            and modify it instead."
+        )
         source = self._sources[name]
-        self._pattern_reg.remove_usage(source.strength_timeseries.pattern_name, (source.name, 'Source'))
-        self._node_reg.remove_usage(source.node_name, (source.name, 'Source'))            
+        self._pattern_reg.remove_usage(source.strength_timeseries.pattern_name, (source.name, "Source"))
+        self._node_reg.remove_usage(source.node_name, (source.name, "Source"))
         del self._sources[name]
-        
-    def remove_control(self, name): 
+
+    def remove_control(self, name):
         """Removes a control from the water network model"""
         del self._controls[name]
 
@@ -663,10 +743,10 @@ def _discard_control(self, name):
             del self._controls[name]
         except KeyError:
             pass
-    
-    ### # 
+
+    ### #
     ### Get elements from the model
-    def get_node(self, name): 
+    def get_node(self, name):
         """Get a specific node
         
         Parameters
@@ -680,8 +760,8 @@ def get_node(self, name):
         
         """
         return self._node_reg[name]
-    
-    def get_link(self, name): 
+
+    def get_link(self, name):
         """Get a specific link
         
         Parameters
@@ -695,8 +775,8 @@ def get_link(self, name):
         
         """
         return self._link_reg[name]
-    
-    def get_pattern(self, name): 
+
+    def get_pattern(self, name):
         """Get a specific pattern
         
         Parameters
@@ -710,8 +790,8 @@ def get_pattern(self, name):
         
         """
         return self._pattern_reg[name]
-    
-    def get_curve(self, name): 
+
+    def get_curve(self, name):
         """Get a specific curve
         
         Parameters
@@ -725,7 +805,7 @@ def get_curve(self, name):
         
         """
         return self._curve_reg[name]
-    
+
     def get_source(self, name):
         """Get a specific source
         
@@ -740,8 +820,8 @@ def get_source(self, name):
         
         """
         return self._sources[name]
-    
-    def get_control(self, name): 
+
+    def get_control(self, name):
         """Get a specific control or rule
         
         Parameters
@@ -755,11 +835,11 @@ def get_control(self, name):
         
         """
         return self._controls[name]
-    
+
     ### #
     ### Name lists
     @property
-    def node_name_list(self): 
+    def node_name_list(self):
         """Get a list of node names
         
         Returns
@@ -770,7 +850,7 @@ def node_name_list(self):
         return list(self._node_reg.keys())
 
     @property
-    def junction_name_list(self): 
+    def junction_name_list(self):
         """Get a list of junction names
         
         Returns
@@ -781,7 +861,7 @@ def junction_name_list(self):
         return list(self._node_reg.junction_names)
 
     @property
-    def tank_name_list(self): 
+    def tank_name_list(self):
         """Get a list of tanks names
         
         Returns
@@ -792,7 +872,7 @@ def tank_name_list(self):
         return list(self._node_reg.tank_names)
 
     @property
-    def reservoir_name_list(self): 
+    def reservoir_name_list(self):
         """Get a list of reservoir names
         
         Returns
@@ -803,7 +883,7 @@ def reservoir_name_list(self):
         return list(self._node_reg.reservoir_names)
 
     @property
-    def link_name_list(self): 
+    def link_name_list(self):
         """Get a list of link names
         
         Returns
@@ -814,7 +894,7 @@ def link_name_list(self):
         return list(self._link_reg.keys())
 
     @property
-    def pipe_name_list(self): 
+    def pipe_name_list(self):
         """Get a list of pipe names
         
         Returns
@@ -935,7 +1015,7 @@ def gpv_name_list(self):
         return list(self._link_reg.gpv_names)
 
     @property
-    def pattern_name_list(self): 
+    def pattern_name_list(self):
         """Get a list of pattern names
         
         Returns
@@ -946,7 +1026,7 @@ def pattern_name_list(self):
         return list(self._pattern_reg.keys())
 
     @property
-    def curve_name_list(self): 
+    def curve_name_list(self):
         """Get a list of curve names
         
         Returns
@@ -957,7 +1037,7 @@ def curve_name_list(self):
         return list(self._curve_reg.keys())
 
     @property
-    def source_name_list(self): 
+    def source_name_list(self):
         """Get a list of source names
         
         Returns
@@ -968,7 +1048,7 @@ def source_name_list(self):
         return list(self._sources.keys())
 
     @property
-    def control_name_list(self): 
+    def control_name_list(self):
         """Get a list of control/rule names
         
         Returns
@@ -977,69 +1057,69 @@ def control_name_list(self):
         
         """
         return list(self._controls.keys())
-    
-    ### # 
+
+    ### #
     ### Counts
     @property
-    def num_nodes(self): 
+    def num_nodes(self):
         """The number of nodes"""
         return len(self._node_reg)
-    
+
     @property
-    def num_junctions(self): 
+    def num_junctions(self):
         """The number of junctions"""
         return len(self._node_reg.junction_names)
-    
+
     @property
-    def num_tanks(self): 
+    def num_tanks(self):
         """The number of tanks"""
         return len(self._node_reg.tank_names)
-    
+
     @property
-    def num_reservoirs(self): 
+    def num_reservoirs(self):
         """The number of reservoirs"""
         return len(self._node_reg.reservoir_names)
-    
+
     @property
-    def num_links(self): 
+    def num_links(self):
         """The number of links"""
         return len(self._link_reg)
-    
+
     @property
-    def num_pipes(self): 
+    def num_pipes(self):
         """The number of pipes"""
         return len(self._link_reg.pipe_names)
-    
+
     @property
-    def num_pumps(self): 
+    def num_pumps(self):
         """The number of pumps"""
         return len(self._link_reg.pump_names)
-    
+
     @property
-    def num_valves(self): 
+    def num_valves(self):
         """The number of valves"""
         return len(self._link_reg.valve_names)
-    
+
     @property
-    def num_patterns(self): 
+    def num_patterns(self):
         """The number of patterns"""
         return len(self._pattern_reg)
-    
+
     @property
-    def num_curves(self): 
+    def num_curves(self):
         """The number of curves"""
         return len(self._curve_reg)
-    
+
     @property
-    def num_sources(self): 
+    def num_sources(self):
         """The number of sources"""
         return len(self._sources)
-    
+
     @property
-    def num_controls(self): 
+    def num_controls(self):
         """The number of controls"""
         return len(self._controls)
-    
+
     ### #
     ### Helper functions
     def describe(self, level=0):
@@ -1063,93 +1143,97 @@ def describe(self, level=0):
         -------
         A dictionary with component counts
         """
-        
-        d = {'Nodes': self.num_nodes,
-             'Links': self.num_links,
-             'Patterns': self.num_patterns,
-             'Curves': self.num_curves,
-             'Sources': self.num_sources,
-             'Controls': self.num_controls}
-        
+
+        d = {
+            "Nodes": self.num_nodes,
+            "Links": self.num_links,
+            "Patterns": self.num_patterns,
+            "Curves": self.num_curves,
+            "Sources": self.num_sources,
+            "Controls": self.num_controls,
+        }
+
         if level >= 1:
-            d['Nodes'] = {
-                    'Junctions': self.num_junctions,
-                    'Tanks': self.num_tanks,
-                    'Reservoirs': self.num_reservoirs}
-            d['Links'] = {
-                    'Pipes': self.num_pipes,
-                    'Pumps': self.num_pumps,
-                    'Valves': self.num_valves}
-            d['Curves'] = {
-                    'Pump': len(self._curve_reg._pump_curves), 
-                    'Efficiency': len(self._curve_reg._efficiency_curves),  
-                    'Headloss': len(self._curve_reg._headloss_curves), 
-                    'Volume': len(self._curve_reg._volume_curves)}
-            
+            d["Nodes"] = {"Junctions": self.num_junctions, "Tanks": self.num_tanks, "Reservoirs": self.num_reservoirs}
+            d["Links"] = {"Pipes": self.num_pipes, "Pumps": self.num_pumps, "Valves": self.num_valves}
+            d["Curves"] = {
+                "Pump": len(self._curve_reg._pump_curves),
+                "Efficiency": len(self._curve_reg._efficiency_curves),
+                "Headloss": len(self._curve_reg._headloss_curves),
+                "Volume": len(self._curve_reg._volume_curves),
+            }
+
         if level >= 2:
-            d['Links']['Pumps'] = {
-                    'Head': len(list(self.head_pumps())),
-                    'Power': len(list(self.power_pumps()))}
-            d['Links']['Valves'] = {
-                    'PRV': len(list(self.prvs())),
-                    'PSV': len(list(self.psvs())),
-                    'PBV': len(list(self.pbvs())),
-                    'TCV': len(list(self.tcvs())),
-                    'FCV': len(list(self.fcvs())),
-                    'GPV': len(list(self.gpvs()))}
-                
+            d["Links"]["Pumps"] = {"Head": len(list(self.head_pumps())), "Power": len(list(self.power_pumps()))}
+            d["Links"]["Valves"] = {
+                "PRV": len(list(self.prvs())),
+                "PSV": len(list(self.psvs())),
+                "PBV": len(list(self.pbvs())),
+                "TCV": len(list(self.tcvs())),
+                "FCV": len(list(self.fcvs())),
+                "GPV": len(list(self.gpvs())),
+            }
+
         return d
-    
+
     def to_dict(self):
-        """Dictionary representation of the WaterNetworkModel.
+        """
+        Dictionary representation of the WaterNetworkModel.
         
         Returns
         -------
         dict
-            Dictionary representation of the WaterNetworkModel
         """
         return wntr.network.io.to_dict(self)
 
     def from_dict(self, d: dict):
         """
-        Append the model with elements from a water network model dictionary.
+        Append the model with elements from a dictionary.
 
         Parameters
         ----------
         d : dict
-            dictionary representation of the water network model to append to existing model
+            Dictionary representation of the WaterNetworkModel
         """
         wntr.network.io.from_dict(d, append=self)
 
-    def write_json(self, f, **kw_json):
+    def to_gis(self, crs=None, pumps_as_points=False, valves_as_points=False):
         """
-        Write the WaterNetworkModel to a JSON file
-
+        Convert a WaterNetworkModel into GeoDataFrames
+        
         Parameters
         ----------
-        f : str
-            Name of the file or file pointer
-        kw_json : keyword arguments
-            arguments to pass directly to `json.dump`
+        crs : str, optional
+            Coordinate reference system, by default None
+        pumps_as_points : bool, optional
+            Represent pumps as points (True) or lines (False), by default False
+        valves_as_points : bool, optional
+            Represent valves as points (True) or lines (False), by default False
+            
+        Returns
+        -------
+        WaterNetworkGIS object that contains junctions, tanks, reservoirs, pipes, 
+        pumps, and valves GeoDataFrames
         """
-        wntr.network.io.write_json(self, f, **kw_json)
+        return wntr.network.io.to_gis(self, crs, pumps_as_points, valves_as_points)
     
-    def read_json(self, f, **kw_json):
+    def from_gis(self, gis_data):
         """
-        Create a WaterNetworkModel from a JSON file.
-
+        Append the model with elements from GeoDataFrames
+        
         Parameters
         ----------
-        f : str
-            Name of the file or file pointer
-        kw_json : keyword arguments
-            keyword arguments to pass to `json.load`
-
+        gis_data : WaterNetworkGIS or dictionary of GeoDataFrames
+            GeoDataFrames containing water network attributes. If gis_data is a 
+            dictionary, then the keys are junctions, tanks, reservoirs, pipes, 
+            pumps, and valves. If the pumps or valves are Points, they will be 
+            converted to Lines with the same start and end node location.
+            
         Returns
         -------
         WaterNetworkModel
         """
-        return wntr.network.io.read_json(f, append=self, **kw_json)
+        return wntr.network.io.from_gis(gis_data, append=self)
 
     def get_graph(self, node_weight=None, link_weight=None, modify_direction=False):
         """
@@ -1172,45 +1256,41 @@ def get_graph(self, node_weight=None, link_weight=None, modify_direction=False):
         networkx MultiDiGraph
         """
         G = nx.MultiDiGraph()
-        
+
         for name, node in self.nodes():
             G.add_node(name)
-            nx.set_node_attributes(G, name='pos', values={name: node.coordinates})
-            nx.set_node_attributes(G, name='type', values={name: node.node_type})
-            
+            nx.set_node_attributes(G, name="pos", values={name: node.coordinates})
+            nx.set_node_attributes(G, name="type", values={name: node.node_type})
+
             if node_weight is not None:
-                try: # weight nodes
+                try:  # weight nodes
                     value = node_weight[name]
-                    nx.set_node_attributes(G, name='weight', values={name: value})
+                    nx.set_node_attributes(G, name="weight", values={name: value})
                 except:
                     pass
-            
+
         for name, link in self.links():
             start_node = link.start_node_name
             end_node = link.end_node_name
             G.add_edge(start_node, end_node, key=name)
-            nx.set_edge_attributes(G, name='type', 
-                        values={(start_node, end_node, name): link.link_type})
-                
+            nx.set_edge_attributes(G, name="type", values={(start_node, end_node, name): link.link_type})
+
             if link_weight is not None:
-                try: # weight links
+                try:  # weight links
                     value = link_weight[name]
-                    if modify_direction and value < 0: # change the direction of the link and value
+                    if modify_direction and value < 0:  # change the direction of the link and value
                         G.remove_edge(start_node, end_node, name)
                         G.add_edge(end_node, start_node, name)
-                        nx.set_edge_attributes(G, name='type', 
-                                values={(end_node, start_node, name): link.link_type})
-                        nx.set_edge_attributes(G, name='weight', 
-                                values={(end_node, start_node, name): -value})
+                        nx.set_edge_attributes(G, name="type", values={(end_node, start_node, name): link.link_type})
+                        nx.set_edge_attributes(G, name="weight", values={(end_node, start_node, name): -value})
                     else:
-                        nx.set_edge_attributes(G, name='weight', 
-                            values={(start_node, end_node, name): value})
+                        nx.set_edge_attributes(G, name="weight", values={(start_node, end_node, name): value})
                 except:
                     pass
-            
+
         return G
-    
-    def assign_demand(self, demand, pattern_prefix='ResetDemand'):
+
+    def assign_demand(self, demand, pattern_prefix="ResetDemand"):
         """
         Assign demands using values in a DataFrame. 
         
@@ -1234,24 +1314,24 @@ def assign_demand(self, demand, pattern_prefix='ResetDemand'):
             should use a different pattern prefix name.
         """
         for junc_name in demand.columns:
-            
+
             # Extract the node demand pattern and resample to match the pattern timestep
             demand_pattern = demand.loc[:, junc_name]
-            demand_pattern.index = pd.TimedeltaIndex(demand_pattern.index, 's')
-            resample_offset = str(int(self.options.time.pattern_timestep))+'S'
+            demand_pattern.index = pd.TimedeltaIndex(demand_pattern.index, "s")
+            resample_offset = str(int(self.options.time.pattern_timestep)) + "S"
             demand_pattern = demand_pattern.resample(resample_offset).mean() / self.options.hydraulic.demand_multiplier
 
             # Add the pattern
-            # If the pattern name already exists, this fails 
+            # If the pattern name already exists, this fails
             pattern_name = pattern_prefix + junc_name
             self.add_pattern(pattern_name, demand_pattern.tolist())
-            
+
             # Reset base demand
             junction = self.get_node(junc_name)
             junction.demand_timeseries_list.clear()
             junction.demand_timeseries_list.append((1.0, pattern_name))
 
-    def get_links_for_node(self, node_name, flag='ALL'):
+    def get_links_for_node(self, node_name, flag="ALL"):
         """
         Returns a list of links connected to a node
 
@@ -1270,20 +1350,33 @@ def get_links_for_node(self, node_name, flag='ALL'):
         -------
         A list of link names connected to the node
         """
-        link_types = {'Pipe', 'Pump', 'Valve'}
+        link_types = {"Pipe", "Pump", "Valve"}
         link_data = self._node_reg.get_usage(node_name)
         if link_data is None:
             return []
         else:
-            if flag.upper() == 'ALL':
-                return [link_name for link_name, link_type in link_data if link_type in link_types and node_name in {self.get_link(link_name).start_node_name, self.get_link(link_name).end_node_name}]
-            elif flag.upper() == 'INLET':
-                return [link_name for link_name, link_type in link_data if link_type in link_types and node_name == self.get_link(link_name).end_node_name]
-            elif flag.upper() == 'OUTLET':
-                return [link_name for link_name, link_type in link_data if link_type in link_types and node_name == self.get_link(link_name).start_node_name]
+            if flag.upper() == "ALL":
+                return [
+                    link_name
+                    for link_name, link_type in link_data
+                    if link_type in link_types
+                    and node_name in {self.get_link(link_name).start_node_name, self.get_link(link_name).end_node_name}
+                ]
+            elif flag.upper() == "INLET":
+                return [
+                    link_name
+                    for link_name, link_type in link_data
+                    if link_type in link_types and node_name == self.get_link(link_name).end_node_name
+                ]
+            elif flag.upper() == "OUTLET":
+                return [
+                    link_name
+                    for link_name, link_type in link_data
+                    if link_type in link_types and node_name == self.get_link(link_name).start_node_name
+                ]
             else:
-                logger.error('Unrecognized flag: {0}'.format(flag))
-                raise ValueError('Unrecognized flag: {0}'.format(flag))
+                logger.error("Unrecognized flag: {0}".format(flag))
+                raise ValueError("Unrecognized flag: {0}".format(flag))
 
     def query_node_attribute(self, attribute, operation=None, value=None, node_type=None):
         """
@@ -1410,7 +1503,7 @@ def convert_controls_to_rules(self, priority=3):
                 act = control.actions()[0]
                 cond = control.condition
                 rule = Rule(cond, act, priority=priority)
-                self.add_control(name.replace(' ', '_')+'_Rule', rule)
+                self.add_control(name.replace(" ", "_") + "_Rule", rule)
                 self.remove_control(name)
 
     def reset_initial_values(self):
@@ -1430,7 +1523,7 @@ def reset_initial_values(self):
             node._is_isolated = False
 
         for name, node in self.nodes(Tank):
-            node._head = node.init_level+node.elevation
+            node._head = node.init_level + node.elevation
             node._prev_head = node.head
             node._demand = None
             node._leak_demand = None
@@ -1471,74 +1564,31 @@ def reset_initial_values(self):
 
         for name, control in self.controls():
             control._reset()
-    
-    def read_inpfile(self, filename):
-        """
-        Defines water network model components from an EPANET INP file
-
-        Parameters
-        ----------
-        filename : string
-            Name of the INP file.
-
-        """
-        return wntr.network.io.read_inpfile(filename, append=self)
-
-    def write_inpfile(self, filename, units=None, version=2.2, force_coordinates=False):
-        """
-        Writes the current water network model to an EPANET INP file
-
-        .. note::
-
-            By default, WNTR now uses EPANET version 2.2 for the EPANET simulator engine. Thus,
-            The WaterNetworkModel will also write an EPANET 2.2 formatted INP file by default as well.
-            Because the PDD analysis options will break EPANET 2.0, the ``version`` option will allow
-            the user to force EPANET 2.0 compatibility at the expense of pressured-dependent analysis 
-            options being turned off.
-
-
-        Parameters
-        ----------
-        filename : string
-            Name of the inp file.
-
-        units : str, int or FlowUnits
-            Name of the units being written to the inp file.
-
-        version : float, {2.0, **2.2**}
-            Optionally specify forcing EPANET 2.0 compatibility.
-
-        force_coordinates : bool
-            This only applies if `self.options.graphics.map_filename` is not `None`,
-            and will force the COORDINATES section to be written even if a MAP file is
-            provided. False by default, but coordinates **are** written by default since
-            the MAP file is `None` by default.
-
-        """
-        wntr.network.io.write_inpfile(self, filename, units=units, version=version, force_coordinates=force_coordinates)
 
 
 class PatternRegistry(Registry):
     """A registry for patterns."""
+
     def _finalize_(self, model):
         super()._finalize_(model)
         self._pattern_reg = None
 
     class DefaultPattern(object):
         """An object that always points to the current default pattern for a model"""
+
         def __init__(self, options):
             self._options = options
-        
+
         def __str__(self):
-            return str(self._options.hydraulic.pattern) if self._options.hydraulic.pattern is not None else ''
-        
+            return str(self._options.hydraulic.pattern) if self._options.hydraulic.pattern is not None else ""
+
         def __repr__(self):
-            return 'DefaultPattern()'
-        
+            return "DefaultPattern()"
+
         @property
         def name(self):
             """The name of the default pattern, or ``None`` if no pattern is assigned"""
-            return str(self._options.hydraulic.pattern) if self._options.hydraulic.pattern is not None else ''
+            return str(self._options.hydraulic.pattern) if self._options.hydraulic.pattern is not None else ""
 
     def __getitem__(self, key):
         try:
@@ -1581,17 +1631,19 @@ def add_pattern(self, name, pattern=None):
         ValueError
             If adding a pattern with `name` that already exists.
         """
-        assert isinstance(name, str) and len(name) < 32 and name.find(' ') == -1, "name must be a string with less than 32 characters and contain no spaces"
+        assert (
+            isinstance(name, str) and len(name) < 32 and name.find(" ") == -1
+        ), "name must be a string with less than 32 characters and contain no spaces"
         assert isinstance(pattern, (list, np.ndarray, Pattern)), "pattern must be a list or Pattern"
-                          
+
         if not isinstance(pattern, Pattern):
-            pattern = Pattern(name, multipliers=pattern, time_options=self._options.time)            
-        else: #elif pattern.time_options is None:
+            pattern = Pattern(name, multipliers=pattern, time_options=self._options.time)
+        else:  # elif pattern.time_options is None:
             pattern.time_options = self._options.time
         if pattern.name in self._data.keys():
-            raise ValueError('Pattern name already exists')
+            raise ValueError("Pattern name already exists")
         self[name] = pattern
-    
+
     @property
     def default_pattern(self):
         """A new default pattern object"""
@@ -1600,6 +1652,7 @@ def default_pattern(self):
 
 class CurveRegistry(Registry):
     """A registry for curves."""
+
     def __init__(self, model):
         super(CurveRegistry, self).__init__(model)
         self._pump_curves = OrderedSet()
@@ -1613,11 +1666,11 @@ def _finalize_(self, model):
 
     def __setitem__(self, key, value):
         if not isinstance(key, six.string_types):
-            raise ValueError('Registry keys must be strings')
+            raise ValueError("Registry keys must be strings")
         self._data[key] = value
         if value is not None:
             self.set_curve_type(key, value.curve_type)
-    
+
     def set_curve_type(self, key, curve_type):
         """
         Sets curve type.
@@ -1627,17 +1680,17 @@ def set_curve_type(self, key, curve_type):
         if curve_type is None:
             return
         curve_type = curve_type.upper()
-        if curve_type == 'HEAD':
+        if curve_type == "HEAD":
             self._pump_curves.add(key)
-        elif curve_type == 'HEADLOSS':
+        elif curve_type == "HEADLOSS":
             self._headloss_curves.add(key)
-        elif curve_type == 'VOLUME':
+        elif curve_type == "VOLUME":
             self._volume_curves.add(key)
-        elif curve_type == 'EFFICIENCY':
+        elif curve_type == "EFFICIENCY":
             self._efficiency_curves.add(key)
         else:
-            raise ValueError('curve_type must be HEAD, HEADLOSS, VOLUME, or EFFICIENCY')
-        
+            raise ValueError("curve_type must be HEAD, HEADLOSS, VOLUME, or EFFICIENCY")
+
     def add_curve(self, name, curve_type, xy_tuples_list):
         """
         Adds a curve to the water network model.
@@ -1651,13 +1704,15 @@ def add_curve(self, name, curve_type, xy_tuples_list):
         xy_tuples_list : list of (x, y) tuples
             List of X-Y coordinate tuples on the curve.
         """
-        assert isinstance(name, str) and len(name) < 32 and name.find(' ') == -1, "name must be a string with less than 32 characters and contain no spaces"
+        assert (
+            isinstance(name, str) and len(name) < 32 and name.find(" ") == -1
+        ), "name must be a string with less than 32 characters and contain no spaces"
         assert isinstance(curve_type, (type(None), str)), "curve_type must be a string"
         assert isinstance(xy_tuples_list, (list, np.ndarray)), "xy_tuples_list must be a list of (x,y) tuples"
-        
+
         curve = Curve(name, curve_type, xy_tuples_list)
         self[name] = curve
-        
+
     def untyped_curves(self):
         """Generator to get all curves without type
         
@@ -1670,17 +1725,19 @@ def untyped_curves(self):
             
         """
         defined = set(self._data.keys())
-        untyped = defined.difference(self._pump_curves, self._efficiency_curves, 
-                                     self._headloss_curves, self._volume_curves)
+        untyped = defined.difference(
+            self._pump_curves, self._efficiency_curves, self._headloss_curves, self._volume_curves
+        )
         for key in untyped:
             yield key, self._data[key]
 
-    @property    
+    @property
     def untyped_curve_names(self):
         """List of names of all curves without types"""
         defined = set(self._data.keys())
-        untyped = defined.difference(self._pump_curves, self._efficiency_curves, 
-                                     self._headloss_curves, self._volume_curves)
+        untyped = defined.difference(
+            self._pump_curves, self._efficiency_curves, self._headloss_curves, self._volume_curves
+        )
         return list(untyped)
 
     def pump_curves(self):
@@ -1696,7 +1753,7 @@ def pump_curves(self):
         """
         for key in self._pump_curves:
             yield key, self._data[key]
-    
+
     @property
     def pump_curve_names(self):
         """List of names of all pump curves"""
@@ -1753,7 +1810,7 @@ def volume_curves(self):
         """
         for key in self._volume_curves:
             yield key, self._data[key]
-    
+
     @property
     def volume_curve_names(self):
         """List of names of all volume curves"""
@@ -1762,6 +1819,7 @@ def volume_curve_names(self):
 
 class SourceRegistry(Registry):
     """A registry for sources."""
+
     def _finalize_(self, model):
         super()._finalize_(model)
         self._sources = None
@@ -1769,15 +1827,14 @@ def _finalize_(self, model):
     def __delitem__(self, key):
         try:
             if self._usage and key in self._usage and len(self._usage[key]) > 0:
-                raise RuntimeError('cannot remove %s %s, still used by %s'%( 
-                                   self.__class__.__name__,
-                                   key,
-                                   self._usage[key]))
+                raise RuntimeError(
+                    "cannot remove %s %s, still used by %s" % (self.__class__.__name__, key, self._usage[key])
+                )
             elif key in self._usage:
                 self._usage.pop(key)
             source = self._data.pop(key)
-            self._pattern_reg.remove_usage(source.strength_timeseries.pattern_name, (source.name, 'Source'))
-            self._node_reg.remove_usage(source.node_name, (source.name, 'Source'))            
+            self._pattern_reg.remove_usage(source.strength_timeseries.pattern_name, (source.name, "Source"))
+            self._node_reg.remove_usage(source.node_name, (source.name, "Source"))
             return source
         except KeyError:
             # Do not raise an exception if there is no key of that name
@@ -1786,19 +1843,20 @@ def __delitem__(self, key):
 
 class NodeRegistry(Registry):
     """A registry for nodes."""
+
     def __init__(self, model):
         super(NodeRegistry, self).__init__(model)
         self._junctions = OrderedSet()
         self._reservoirs = OrderedSet()
         self._tanks = OrderedSet()
-    
+
     def _finalize_(self, model):
         super()._finalize_(model)
         self._node_reg = None
-    
+
     def __setitem__(self, key, value):
         if not isinstance(key, six.string_types):
-            raise ValueError('Registry keys must be strings')
+            raise ValueError("Registry keys must be strings")
         self._data[key] = value
         if isinstance(value, Junction):
             self._junctions.add(key)
@@ -1806,14 +1864,13 @@ def __setitem__(self, key, value):
             self._tanks.add(key)
         elif isinstance(value, Reservoir):
             self._reservoirs.add(key)
-    
+
     def __delitem__(self, key):
         try:
             if self._usage and key in self._usage and len(self._usage[key]) > 0:
-                raise RuntimeError('cannot remove %s %s, still used by %s'%(
-                                   self.__class__.__name__,
-                                   key,
-                                   str(self._usage[key])))
+                raise RuntimeError(
+                    "cannot remove %s %s, still used by %s" % (self.__class__.__name__, key, str(self._usage[key]))
+                )
             elif key in self._usage:
                 self._usage.pop(key)
             node = self._data.pop(key)
@@ -1823,15 +1880,15 @@ def __delitem__(self, key):
             if isinstance(node, Junction):
                 for pat_name in node.demand_timeseries_list.pattern_list():
                     if pat_name:
-                        self._curve_reg.remove_usage(pat_name, (node.name, 'Junction'))
+                        self._curve_reg.remove_usage(pat_name, (node.name, "Junction"))
             if isinstance(node, Reservoir) and node.head_pattern_name:
-                self._curve_reg.remove_usage(node.head_pattern_name, (node.name, 'Reservoir'))
+                self._curve_reg.remove_usage(node.head_pattern_name, (node.name, "Reservoir"))
             if isinstance(node, Tank) and node.vol_curve_name:
-                self._curve_reg.remove_usage(node.vol_curve_name, (node.name, 'Tank'))
+                self._curve_reg.remove_usage(node.vol_curve_name, (node.name, "Tank"))
             return node
         except KeyError:
-            return 
-    
+            return
+
     def __call__(self, node_type=None):
         """
         Returns a generator to iterate over all nodes of a specific node type.
@@ -1851,24 +1908,32 @@ def __call__(self, node_type=None):
         -------
         A generator in the format (name, object).
         """
-        if node_type==None:
+        if node_type == None:
             for node_name, node in self._data.items():
                 yield node_name, node
-        elif node_type==Junction:
+        elif node_type == Junction:
             for node_name in self._junctions:
                 yield node_name, self._data[node_name]
-        elif node_type==Tank:
+        elif node_type == Tank:
             for node_name in self._tanks:
                 yield node_name, self._data[node_name]
-        elif node_type==Reservoir:
+        elif node_type == Reservoir:
             for node_name in self._reservoirs:
                 yield node_name, self._data[node_name]
         else:
-            raise RuntimeError('node_type, '+str(node_type)+', not recognized.')
-
-    def add_junction(self, name, base_demand=0.0, demand_pattern=None, 
-                     elevation=0.0, coordinates=None, demand_category=None,
-                     emitter_coeff=None, initial_quality=None):
+            raise RuntimeError("node_type, " + str(node_type) + ", not recognized.")
+
+    def add_junction(
+        self,
+        name,
+        base_demand=0.0,
+        demand_pattern=None,
+        elevation=0.0,
+        coordinates=None,
+        demand_category=None,
+        emitter_coeff=None,
+        initial_quality=None,
+    ):
         """
         Adds a junction to the water network model.
 
@@ -1891,18 +1956,22 @@ def add_junction(self, name, base_demand=0.0, demand_pattern=None,
         initial_quality : float, optional
             Initial quality at this junction
         """
-        assert isinstance(name, str) and len(name) < 32 and name.find(' ') == -1, "name must be a string with less than 32 characters and contain no spaces"
+        assert (
+            isinstance(name, str) and len(name) < 32 and name.find(" ") == -1
+        ), "name must be a string with less than 32 characters and contain no spaces"
         assert isinstance(base_demand, (int, float)), "base_demand must be a float"
-        assert isinstance(demand_pattern, (type(None), str, PatternRegistry.DefaultPattern, Pattern)), "demand_pattern must be a string or Pattern"
+        assert isinstance(
+            demand_pattern, (type(None), str, PatternRegistry.DefaultPattern, Pattern)
+        ), "demand_pattern must be a string or Pattern"
         assert isinstance(elevation, (int, float)), "elevation must be a float"
-        assert isinstance(coordinates, (type(None), (tuple,list,))), "coordinates must be a tuple"
+        assert isinstance(coordinates, (type(None), (tuple, list,))), "coordinates must be a tuple"
         assert isinstance(demand_category, (type(None), str)), "demand_category must be a string"
         assert isinstance(emitter_coeff, (type(None), int, float)), "emitter_coeff must be a float"
         assert isinstance(initial_quality, (type(None), int, float)), "initial_quality must be a float"
-        
+
         base_demand = float(base_demand)
         elevation = float(elevation)
-        
+
         junction = Junction(name, self)
         junction.elevation = elevation
         junction.add_demand(base_demand, demand_pattern, demand_category)
@@ -1914,10 +1983,19 @@ def add_junction(self, name, base_demand=0.0, demand_pattern=None,
         if initial_quality is not None:
             junction.initial_quality = initial_quality
 
-    def add_tank(self, name, elevation=0.0, init_level=3.048,
-                 min_level=0.0, max_level=6.096, diameter=15.24,
-                 min_vol=0.0, vol_curve=None, overflow=False, 
-                 coordinates=None):
+    def add_tank(
+        self,
+        name,
+        elevation=0.0,
+        init_level=3.048,
+        min_level=0.0,
+        max_level=6.096,
+        diameter=15.24,
+        min_vol=0.0,
+        vol_curve=None,
+        overflow=False,
+        coordinates=None,
+    ):
         """
         Adds a tank to the water network model.
 
@@ -1947,7 +2025,9 @@ def add_tank(self, name, elevation=0.0, init_level=3.048,
             X-Y coordinates of the node location.
             
         """
-        assert isinstance(name, str) and len(name) < 32 and name.find(' ') == -1, "name must be a string with less than 32 characters and contain no spaces"
+        assert (
+            isinstance(name, str) and len(name) < 32 and name.find(" ") == -1
+        ), "name must be a string with less than 32 characters and contain no spaces"
         assert isinstance(elevation, (int, float)), "elevation must be a float"
         assert isinstance(init_level, (int, float)), "init_level must be a float"
         assert isinstance(min_level, (int, float)), "min_level must be a float"
@@ -1968,22 +2048,42 @@ def add_tank(self, name, elevation=0.0, init_level=3.048,
             raise ValueError("Initial tank level must be greater than or equal to the tank minimum level.")
         if init_level > max_level:
             raise ValueError("Initial tank level must be less than or equal to the tank maximum level.")
-        if vol_curve is not None and vol_curve != '*':
+        if vol_curve is not None and vol_curve != "*":
             if not isinstance(vol_curve, six.string_types):
-                raise ValueError('Volume curve name must be a string')
+                raise ValueError("Volume curve name must be a string")
             elif not vol_curve in self._curve_reg.volume_curve_names:
-                raise ValueError('The volume curve ' + vol_curve + ' is not one of the curves in the ' +
-                                 'list of volume curves. Valid volume curves are:' + 
-                                 str(self._curve_reg.volume_curve_names))
+                raise ValueError(
+                    "The volume curve "
+                    + vol_curve
+                    + " is not one of the curves in the "
+                    + "list of volume curves. Valid volume curves are:"
+                    + str(self._curve_reg.volume_curve_names)
+                )
             vcurve = np.array(self._curve_reg[vol_curve].points)
-            if min_level < vcurve[0,0]:
-                raise ValueError(('The volume curve ' + vol_curve + ' has a minimum value ({0:5.2f}) \n' +
-                                 'greater than the minimum level for tank "' + name + '" ({1:5.2f})\n' +
-                                 'please correct the user input.').format(vcurve[0,0],min_level))
-            elif max_level > vcurve[-1,0]:
-                raise ValueError(('The volume curve ' + vol_curve + ' has a maximum value ({0:5.2f}) \n' +
-                                 'less than the maximum level for tank "' + name + '" ({1:5.2f})\n' +
-                                 'please correct the user input.').format(vcurve[-1,0],max_level))
+            if min_level < vcurve[0, 0]:
+                raise ValueError(
+                    (
+                        "The volume curve "
+                        + vol_curve
+                        + " has a minimum value ({0:5.2f}) \n"
+                        + 'greater than the minimum level for tank "'
+                        + name
+                        + '" ({1:5.2f})\n'
+                        + "please correct the user input."
+                    ).format(vcurve[0, 0], min_level)
+                )
+            elif max_level > vcurve[-1, 0]:
+                raise ValueError(
+                    (
+                        "The volume curve "
+                        + vol_curve
+                        + " has a maximum value ({0:5.2f}) \n"
+                        + 'less than the maximum level for tank "'
+                        + name
+                        + '" ({1:5.2f})\n'
+                        + "please correct the user input."
+                    ).format(vcurve[-1, 0], max_level)
+                )
 
         tank = Tank(name, self)
         tank.elevation = elevation
@@ -2014,11 +2114,13 @@ def add_reservoir(self, name, base_head=0.0, head_pattern=None, coordinates=None
             X-Y coordinates of the node location.
         
         """
-        assert isinstance(name, str) and len(name) < 32 and name.find(' ') == -1, "name must be a string with less than 32 characters and contain no spaces"
+        assert (
+            isinstance(name, str) and len(name) < 32 and name.find(" ") == -1
+        ), "name must be a string with less than 32 characters and contain no spaces"
         assert isinstance(base_head, (int, float)), "base_head must be float"
         assert isinstance(head_pattern, (type(None), str)), "head_pattern must be a string"
         assert isinstance(coordinates, (type(None), (tuple, list))), "coordinates must be a tuple"
-        
+
         base_head = float(base_head)
 
         reservoir = Reservoir(name, self)
@@ -2032,17 +2134,17 @@ def add_reservoir(self, name, base_head=0.0, head_pattern=None, coordinates=None
     def junction_names(self):
         """List of names of all junctions"""
         return self._junctions
-    
+
     @property
     def tank_names(self):
         """List of names of all junctions"""
         return self._tanks
-    
+
     @property
     def reservoir_names(self):
         """List of names of all junctions"""
         return self._reservoirs
-    
+
     def junctions(self):
         """Generator to get all junctions
         
@@ -2056,7 +2158,7 @@ def junctions(self):
         """
         for node_name in self._junctions:
             yield node_name, self._data[node_name]
-    
+
     def tanks(self):
         """Generator to get all tanks
         
@@ -2070,7 +2172,7 @@ def tanks(self):
         """
         for node_name in self._tanks:
             yield node_name, self._data[node_name]
-    
+
     def reservoirs(self):
         """Generator to get all reservoirs
         
@@ -2088,7 +2190,20 @@ def reservoirs(self):
 
 class LinkRegistry(Registry):
     """A registry for links."""
-    __subsets = ['_pipes', '_pumps', '_head_pumps', '_power_pumps', '_prvs', '_psvs', '_pbvs', '_tcvs', '_fcvs', '_gpvs', '_valves']
+
+    __subsets = [
+        "_pipes",
+        "_pumps",
+        "_head_pumps",
+        "_power_pumps",
+        "_prvs",
+        "_psvs",
+        "_pbvs",
+        "_tcvs",
+        "_fcvs",
+        "_gpvs",
+        "_valves",
+    ]
 
     def __init__(self, model):
         super(LinkRegistry, self).__init__(model)
@@ -2103,14 +2218,14 @@ def __init__(self, model):
         self._fcvs = OrderedSet()
         self._gpvs = OrderedSet()
         self._valves = OrderedSet()
-    
+
     def _finalize_(self, model):
         super()._finalize_(model)
         self._link_reg = None
 
     def __setitem__(self, key, value):
         if not isinstance(key, six.string_types):
-            raise ValueError('Registry keys must be strings')
+            raise ValueError("Registry keys must be strings")
         self._data[key] = value
         if isinstance(value, Pipe):
             self._pipes.add(key)
@@ -2134,32 +2249,31 @@ def __setitem__(self, key, value):
                 self._fcvs.add(key)
             elif isinstance(value, GPValve):
                 self._gpvs.add(key)
-    
+
     def __delitem__(self, key):
         try:
             if self._usage and key in self._usage and len(self._usage[key]) > 0:
-                raise RuntimeError('cannot remove %s %s, still used by %s', 
-                                   self.__class__.__name__,
-                                   key,
-                                   self._usage[key])
+                raise RuntimeError(
+                    "cannot remove %s %s, still used by %s", self.__class__.__name__, key, self._usage[key]
+                )
             elif key in self._usage:
                 self._usage.pop(key)
             link = self._data.pop(key)
             self._node_reg.remove_usage(link.start_node_name, (link.name, link.link_type))
             self._node_reg.remove_usage(link.end_node_name, (link.name, link.link_type))
             if isinstance(link, GPValve):
-                self._curve_reg.remove_usage(link.headloss_curve_name, (link.name, 'Valve'))
+                self._curve_reg.remove_usage(link.headloss_curve_name, (link.name, "Valve"))
             if isinstance(link, Pump):
-                self._curve_reg.remove_usage(link.speed_pattern_name, (link.name, 'Pump'))
+                self._curve_reg.remove_usage(link.speed_pattern_name, (link.name, "Pump"))
             if isinstance(link, HeadPump):
-                self._curve_reg.remove_usage(link.pump_curve_name, (link.name, 'Pump'))
+                self._curve_reg.remove_usage(link.pump_curve_name, (link.name, "Pump"))
             for ss in self.__subsets:
                 # Go through the _pipes, _prvs, ..., and remove this link
                 getattr(self, ss).discard(key)
             return link
         except KeyError:
             return
-    
+
     def __call__(self, link_type=None):
         """
         Returns a generator to iterate over all nodes of a specific node type.
@@ -2179,23 +2293,33 @@ def __call__(self, link_type=None):
         -------
         A generator in the format (name, object).
         """
-        if link_type==None:
+        if link_type == None:
             for name, node in self._data.items():
                 yield name, node
-        elif link_type==Pipe:
+        elif link_type == Pipe:
             for name in self._pipes:
                 yield name, self._data[name]
-        elif link_type==Pump:
+        elif link_type == Pump:
             for name in self._pumps:
                 yield name, self._data[name]
-        elif link_type==Valve:
+        elif link_type == Valve:
             for name in self._valves:
                 yield name, self._data[name]
         else:
-            raise RuntimeError('link_type, '+str(link_type)+', not recognized.')
-
-    def add_pipe(self, name, start_node_name, end_node_name, length=304.8,
-                 diameter=0.3048, roughness=100, minor_loss=0.0, initial_status='OPEN', check_valve=False):
+            raise RuntimeError("link_type, " + str(link_type) + ", not recognized.")
+
+    def add_pipe(
+        self,
+        name,
+        start_node_name,
+        end_node_name,
+        length=304.8,
+        diameter=0.3048,
+        roughness=100,
+        minor_loss=0.0,
+        initial_status="OPEN",
+        check_valve=False,
+    ):
         """
         Adds a pipe to the water network model.
 
@@ -2222,15 +2346,22 @@ def add_pipe(self, name, start_node_name, end_node_name, length=304.8,
             False if the pipe does not have a check valve.
         
         """
-        assert isinstance(name, str) and len(name) < 32 and name.find(' ') == -1, "name must be a string with less than 32 characters and contain no spaces"
-        assert isinstance(start_node_name, str) and len(start_node_name) < 32 and start_node_name.find(' ') == -1, "start_node_name must be a string with less than 32 characters and contain no spaces"
-        assert isinstance(end_node_name, str) and len(end_node_name) < 32 and end_node_name.find(' ') == -1, "end_node_name must be a string with less than 32 characters and contain no spaces"
+        assert (
+            isinstance(name, str) and len(name) < 32 and name.find(" ") == -1
+        ), "name must be a string with less than 32 characters and contain no spaces"
+        assert (
+            isinstance(start_node_name, str) and len(start_node_name) < 32 and start_node_name.find(" ") == -1
+        ), "start_node_name must be a string with less than 32 characters and contain no spaces"
+        assert (
+            isinstance(end_node_name, str) and len(end_node_name) < 32 and end_node_name.find(" ") == -1
+        ), "end_node_name must be a string with less than 32 characters and contain no spaces"
         assert isinstance(length, (int, float)), "length must be a float"
         assert isinstance(diameter, (int, float)), "diameter must be a float"
         assert isinstance(roughness, (int, float)), "roughness must be a float"
         assert isinstance(minor_loss, (int, float)), "minor_loss must be a float"
-        assert isinstance(initial_status, (str, LinkStatus)), "initial_status must be a string or LinkStatus"
-        assert isinstance(check_valve, bool), "check_valve must be a Boolean"
+        assert isinstance(initial_status, (int, str, LinkStatus)), "initial_status must be an int, string or LinkStatus"
+        assert isinstance(check_valve, (str, int, bool)), "check_valve must be a Boolean"
+        check_valve = bool(int(check_valve))
         
         length = float(length)
         diameter = float(diameter)
@@ -2238,7 +2369,7 @@ def add_pipe(self, name, start_node_name, end_node_name, length=304.8,
         minor_loss = float(minor_loss)
         if isinstance(initial_status, str):
             initial_status = LinkStatus[initial_status]
-            
+
         pipe = Pipe(name, start_node_name, end_node_name, self)
         pipe.length = length
         pipe.diameter = diameter
@@ -2249,8 +2380,17 @@ def add_pipe(self, name, start_node_name, end_node_name, length=304.8,
         pipe.check_valve = check_valve
         self[name] = pipe
 
-    def add_pump(self, name, start_node_name, end_node_name, pump_type='POWER',
-                 pump_parameter=50.0, speed=1.0, pattern=None, initial_status='OPEN'):
+    def add_pump(
+        self,
+        name,
+        start_node_name,
+        end_node_name,
+        pump_type="POWER",
+        pump_parameter=50.0,
+        speed=1.0,
+        pattern=None,
+        initial_status="OPEN",
+    ):
         """
         Adds a pump to the water network model.
 
@@ -2275,21 +2415,27 @@ def add_pump(self, name, start_node_name, end_node_name, pump_type='POWER',
             Pump initial status. Options are 'OPEN' or 'CLOSED'.
         
         """
-        assert isinstance(name, str) and len(name) < 32 and name.find(' ') == -1, "name must be a string with less than 32 characters and contain no spaces"
-        assert isinstance(start_node_name, str) and len(start_node_name) < 32 and start_node_name.find(' ') == -1, "start_node_name must be a string with less than 32 characters and contain no spaces"
-        assert isinstance(end_node_name, str) and len(end_node_name) < 32 and end_node_name.find(' ') == -1, "end_node_name must be a string with less than 32 characters and contain no spaces"
+        assert (
+            isinstance(name, str) and len(name) < 32 and name.find(" ") == -1
+        ), "name must be a string with less than 32 characters and contain no spaces"
+        assert (
+            isinstance(start_node_name, str) and len(start_node_name) < 32 and start_node_name.find(" ") == -1
+        ), "start_node_name must be a string with less than 32 characters and contain no spaces"
+        assert (
+            isinstance(end_node_name, str) and len(end_node_name) < 32 and end_node_name.find(" ") == -1
+        ), "end_node_name must be a string with less than 32 characters and contain no spaces"
         assert isinstance(pump_type, str), "pump_type must be a string"
         assert isinstance(pump_parameter, (int, float, str)), "pump_parameter must be a float or string"
         assert isinstance(speed, (int, float)), "speed must be a float"
         assert isinstance(pattern, (type(None), str)), "pattern must be a string"
-        assert isinstance(initial_status, (str, LinkStatus)), "initial_status must be a string or LinkStatus"
-        
+        assert isinstance(initial_status, (int, str, LinkStatus)), "initial_status must be an int, string or LinkStatus"
+
         if isinstance(initial_status, str):
             initial_status = LinkStatus[initial_status]
-        if pump_type.upper() == 'POWER':
+        if pump_type.upper() == "POWER":
             pump = PowerPump(name, start_node_name, end_node_name, self)
             pump.power = pump_parameter
-        elif pump_type.upper() == 'HEAD':
+        elif pump_type.upper() == "HEAD":
             pump = HeadPump(name, start_node_name, end_node_name, self)
             pump.pump_curve_name = pump_parameter
         else:
@@ -2298,10 +2444,18 @@ def add_pump(self, name, start_node_name, end_node_name, pump_type='POWER',
         pump.initial_status = initial_status
         pump.speed_pattern_name = pattern
         self[name] = pump
-    
-    def add_valve(self, name, start_node_name, end_node_name,
-                 diameter=0.3048, valve_type='PRV', minor_loss=0.0, 
-                 initial_setting=0.0, initial_status='ACTIVE'):
+
+    def add_valve(
+        self,
+        name,
+        start_node_name,
+        end_node_name,
+        diameter=0.3048,
+        valve_type="PRV",
+        minor_loss=0.0,
+        initial_setting=0.0,
+        initial_status="ACTIVE",
+    ):
         """
         Adds a valve to the water network model.
 
@@ -2329,60 +2483,77 @@ def add_valve(self, name, start_node_name, end_node_name,
             Valve initial status. Options are 'OPEN',  'CLOSED', or 'ACTIVE'
             
         """
-        assert isinstance(name, str) and len(name) < 32 and name.find(' ') == -1, "name must be a string with less than 32 characters and contain no spaces"
-        assert isinstance(start_node_name, str) and len(start_node_name) < 32 and start_node_name.find(' ') == -1, "start_node_name must be a string with less than 32 characters and contain no spaces"
-        assert isinstance(end_node_name, str) and len(end_node_name) < 32 and end_node_name.find(' ') == -1, "end_node_name must be a string with less than 32 characters and contain no spaces"
+        assert (
+            isinstance(name, str) and len(name) < 32 and name.find(" ") == -1
+        ), "name must be a string with less than 32 characters and contain no spaces"
+        assert (
+            isinstance(start_node_name, str) and len(start_node_name) < 32 and start_node_name.find(" ") == -1
+        ), "start_node_name must be a string with less than 32 characters and contain no spaces"
+        assert (
+            isinstance(end_node_name, str) and len(end_node_name) < 32 and end_node_name.find(" ") == -1
+        ), "end_node_name must be a string with less than 32 characters and contain no spaces"
         assert isinstance(diameter, (int, float)), "diameter must be a float"
         assert isinstance(valve_type, str), "valve_type must be a string"
         assert isinstance(minor_loss, (int, float)), "minor_loss must be a float"
         assert isinstance(initial_setting, (int, float, str)), "initial_setting must be a float or string"
         assert isinstance(initial_status, (str, LinkStatus)), "initial_status must be a string or LinkStatus"
-        
+
         if isinstance(initial_status, str):
             initial_status = LinkStatus[initial_status]
         start_node = self._node_reg[start_node_name]
         end_node = self._node_reg[end_node_name]
-        
+
         valve_type = valve_type.upper()
-        
+
         # A PRV, PSV or FCV cannot be directly connected to a reservoir or tank (use a length of pipe to separate the two)
-        if valve_type in ['PRV', 'PSV', 'FCV']:
-            if type(start_node)==Tank or type(end_node)==Tank or type(start_node)==Reservoir or type(end_node)==Reservoir:
-                msg = '%ss cannot be directly connected to a tank.  Add a pipe to separate the valve from the tank.' % valve_type
-                logger.error(msg)   
+        if valve_type in ["PRV", "PSV", "FCV"]:
+            if (
+                type(start_node) == Tank
+                or type(end_node) == Tank
+                or type(start_node) == Reservoir
+                or type(end_node) == Reservoir
+            ):
+                msg = (
+                    "%ss cannot be directly connected to a tank.  Add a pipe to separate the valve from the tank."
+                    % valve_type
+                )
+                logger.error(msg)
                 raise RuntimeError(msg)
-            if type(start_node)==Reservoir or type(end_node)==Reservoir:
-                msg = '%ss cannot be directly connected to a reservoir.  Add a pipe to separate the valve from the reservoir.' % valve_type
-                logger.error(msg)   
+            if type(start_node) == Reservoir or type(end_node) == Reservoir:
+                msg = (
+                    "%ss cannot be directly connected to a reservoir.  Add a pipe to separate the valve from the reservoir."
+                    % valve_type
+                )
+                logger.error(msg)
                 raise RuntimeError(msg)
-        
+
         # TODO check the following: PRVs cannot share the same downstream node or be linked in series
-            
+
         # TODO check the following: Two PSVs cannot share the same upstream node or be linked in series
-        
+
         # TODO check the following: A PSV cannot be connected to the downstream node of a PRV
-        
-        if valve_type == 'PRV':
+
+        if valve_type == "PRV":
             valve = PRValve(name, start_node_name, end_node_name, self)
             valve.initial_setting = initial_setting
             valve._setting = initial_setting
-        elif valve_type == 'PSV':
+        elif valve_type == "PSV":
             valve = PSValve(name, start_node_name, end_node_name, self)
             valve.initial_setting = initial_setting
             valve._setting = initial_setting
-        elif valve_type == 'PBV':
+        elif valve_type == "PBV":
             valve = PBValve(name, start_node_name, end_node_name, self)
             valve.initial_setting = initial_setting
             valve._setting = initial_setting
-        elif valve_type == 'FCV':
+        elif valve_type == "FCV":
             valve = FCValve(name, start_node_name, end_node_name, self)
             valve.initial_setting = initial_setting
             valve._setting = initial_setting
-        elif valve_type == 'TCV':
+        elif valve_type == "TCV":
             valve = TCValve(name, start_node_name, end_node_name, self)
             valve.initial_setting = initial_setting
             valve._setting = initial_setting
-        elif valve_type == 'GPV':
+        elif valve_type == "GPV":
             valve = GPValve(name, start_node_name, end_node_name, self)
             valve.headloss_curve_name = initial_setting
         valve.initial_status = initial_status
@@ -2409,12 +2580,12 @@ def check_valves(self):
     def pipe_names(self):
         """A list of all pipe names"""
         return self._pipes
-    
+
     @property
     def valve_names(self):
         """A list of all valve names"""
         return self._valves
-    
+
     @property
     def pump_names(self):
         """A list of all pump names"""
@@ -2473,7 +2644,7 @@ def pipes(self):
         """
         for name in self._pipes:
             yield name, self._data[name]
-    
+
     def pumps(self):
         """Generator to get all pumps
         
@@ -2487,7 +2658,7 @@ def pumps(self):
         """
         for name in self._pumps:
             yield name, self._data[name]
-    
+
     def valves(self):
         """Generator to get all valves
         
diff --git a/wntr/sim/epanet.py b/wntr/sim/epanet.py
index d5a395595..934324578 100644
--- a/wntr/sim/epanet.py
+++ b/wntr/sim/epanet.py
@@ -1,4 +1,5 @@
 from wntr.sim.core import WaterNetworkSimulator
+from wntr.network.io import write_inpfile
 import wntr.epanet.io
 import warnings
 import logging
@@ -99,7 +100,7 @@ def run_sim(self, file_prefix='temp', save_hyd=False, use_hyd=False, hydfile=Non
         if isinstance(version, str):
             version = float(version)
         inpfile = file_prefix + '.inp'
-        self._wn.write_inpfile(inpfile, units=self._wn.options.hydraulic.inpfile_units, version=version)
+        write_inpfile(self._wn, inpfile, units=self._wn.options.hydraulic.inpfile_units, version=version)
         enData = wntr.epanet.toolkit.ENepanet(version=version)
         rptfile = file_prefix + '.rpt'
         outfile = file_prefix + '.bin'


Test patch:
diff --git a/.github/workflows/build_tests.yml b/.github/workflows/build_tests.yml
index 532adec37..a41844c7c 100644
--- a/.github/workflows/build_tests.yml
+++ b/.github/workflows/build_tests.yml
@@ -15,8 +15,6 @@ jobs:
     runs-on: ${{ matrix.os }}
     strategy:
       matrix:
-        # python-version: ['3.7', '3.10']
-        # os: [ubuntu-latest]
         python-version: ['3.7', '3.8', '3.9', '3.10']
         os: [windows-latest, macOS-latest, ubuntu-latest] 
     steps:
@@ -47,8 +45,6 @@ jobs:
     runs-on: ${{ matrix.os }}
     strategy:
       matrix:
-        # python-version: ['3.7', '3.10']
-        # os: [ubuntu-latest]
         python-version: ['3.7', '3.8', '3.9', '3.10']
         os: [windows-latest, macOS-latest, ubuntu-latest]    
     steps:
@@ -77,7 +73,6 @@ jobs:
     strategy:
       matrix:
         python-version: ['3.7', '3.10']
-        # os: [ubuntu-latest]
         os: [windows-latest, macOS-latest, ubuntu-latest]
     steps:
     - uses: actions/checkout@v2
@@ -91,14 +86,23 @@ jobs:
         python -m pip install --upgrade pip
         pip install -r requirements.txt
         python setup.py develop
-    # pip install coveralls
+    - name: Run Tests on Windows
+      if: matrix.os == 'windows-latest'
+      run: | # the following skips windows rst doctests that rely on geopandas (gis.rst) because skipif is not working properly
+        coverage erase
+        coverage run --context=${{ matrix.os }}.py${{ matrix.python-version }} --source=wntr --omit="*/tests/*" -m pytest  --doctest-modules --doctest-glob="*.rst" wntr
+        coverage run --context=${{ matrix.os }}.py${{ matrix.python-version }} --source=wntr --omit="*/tests/*" --append -m pytest --doctest-glob="*.rst" --ignore-glob="*gis.rst" documentation
+      env:
+        COVERAGE_FILE: .coverage.${{ matrix.python-version }}.${{ matrix.os }}
     - name: Run Tests
-      run: |
+      if: matrix.os != 'windows-latest'
+      run: | 
         coverage erase
         coverage run --context=${{ matrix.os }}.py${{ matrix.python-version }} --source=wntr --omit="*/tests/*" -m pytest  --doctest-modules --doctest-glob="*.rst" wntr
         coverage run --context=${{ matrix.os }}.py${{ matrix.python-version }} --source=wntr --omit="*/tests/*" --append -m pytest --doctest-glob="*.rst" documentation
       env:
         COVERAGE_FILE: .coverage.${{ matrix.python-version }}.${{ matrix.os }}
+        
     - name: Save coverage
       uses: actions/upload-artifact@v2
       with:
diff --git a/wntr/tests/test_epanet_io.py b/wntr/tests/test_epanet_io.py
index debb821d6..afa5082f3 100644
--- a/wntr/tests/test_epanet_io.py
+++ b/wntr/tests/test_epanet_io.py
@@ -18,7 +18,7 @@ def setUpClass(self):
 
         inp_file = join(test_datadir, "io.inp")
         self.wn = self.wntr.network.WaterNetworkModel(inp_file)
-        self.wn.write_inpfile("temp.inp", "GPM")
+        self.wntr.network.write_inpfile(self.wn, "temp.inp", "GPM")
         self.wn2 = self.wntr.network.WaterNetworkModel("temp.inp")
 
     @classmethod
@@ -100,7 +100,7 @@ def setUpClass(self):
         self.wntr = wntr
         inp_file = join(test_datadir, "Net6_plus.inp")  # UNITS = GPM
         self.wn = wntr.network.WaterNetworkModel(inp_file)
-        self.wn.write_inpfile("temp.inp", units="LPM")
+        self.wntr.network.write_inpfile(self.wn, "temp.inp", units="LPM")
         self.wn2 = self.wntr.network.WaterNetworkModel(inp_file)
 
     @classmethod
@@ -449,7 +449,7 @@ def test_convert_controls_to_rules(self):
         inp_file = join(test_datadir, "Net6_plus.inp")  # UNITS = GPM
         self.wn = self.wntr.network.WaterNetworkModel(inp_file)
         self.wn.convert_controls_to_rules()
-        self.wn.write_inpfile('temp2.inp')
+        self.wntr.network.write_inpfile(self.wn, 'temp2.inp')
         self.wn2 = self.wntr.network.WaterNetworkModel('temp2.inp')
         for control_name in self.wn.control_name_list:
             ctrl1 = self.wn.get_control(control_name)
@@ -460,11 +460,11 @@ def test_convert_single_line_rules(self):
         inp_file = join(test_datadir, "Net6_plus.inp")  # UNITS = GPM
         self.wn = self.wntr.network.WaterNetworkModel(inp_file)
         self.wn.convert_controls_to_rules()
-        self.wn.write_inpfile('temp2.inp')
+        self.wntr.network.write_inpfile(self.wn, 'temp2.inp')
         self.wn2 = self.wntr.network.WaterNetworkModel('temp2.inp')
         for control_name in [n for n in self.wn2.control_name_list]:
             self.wn2.remove_control(control_name)
-        self.wn2.write_inpfile('temp3.inp')
+        self.wntr.network.write_inpfile(self.wn2, 'temp3.inp')
         self.wn3 = self.wntr.network.WaterNetworkModel('temp3.inp')
         rules = self._rules.splitlines()
         new_rules = self.wntr.epanet.io._EpanetRule.parse_rules_lines(rules, self.wntr.epanet.util.FlowUnits.GPM)
@@ -488,7 +488,7 @@ def setUpClass(self):
         self.wn = wntr.network.WaterNetworkModel(inp_file)
         self.wn.options.hydraulic.demand_model = "PDA"
         self.wn.options.hydraulic.required_pressure = 1.0
-        self.wn.write_inpfile("temp2.inp", units="LPM")
+        self.wntr.network.write_inpfile(self.wn, "temp2.inp", units="LPM")
         self.wn2 = self.wntr.network.WaterNetworkModel("temp2.inp")
 
     @classmethod
@@ -513,13 +513,13 @@ def setUpClass(self):
         sim = self.wntr.sim.EpanetSimulator(self.wn)
         self.results = sim.run_sim()
 
-        self.wn.write_inpfile("temp.inp")
+        self.wntr.network.write_inpfile(self.wn, "temp.inp")
         self.wn2 = self.wntr.network.WaterNetworkModel("temp.inp")
 
         sim = self.wntr.sim.EpanetSimulator(self.wn2)
         self.results2 = sim.run_sim()
 
-        self.wn.write_inpfile("temp.inp")
+        self.wntr.network.write_inpfile(self.wn, "temp.inp")
         self.wn22 = self.wntr.network.WaterNetworkModel("temp.inp")
         self.wn22.options.hydraulic.demand_model = "PDA"
 
@@ -612,7 +612,7 @@ def setUpClass(self):
         sim = self.wntr.sim.EpanetSimulator(self.wn)
         self.results = sim.run_sim()
 
-        self.wn.write_inpfile("temp.inp", units="CMH")
+        self.wntr.network.write_inpfile(self.wn, "temp.inp", units="CMH")
         self.wn2 = self.wntr.network.WaterNetworkModel("temp.inp")
 
         sim = self.wntr.sim.EpanetSimulator(self.wn2)
diff --git a/wntr/tests/test_gis.py b/wntr/tests/test_gis.py
new file mode 100644
index 000000000..51a4a50dc
--- /dev/null
+++ b/wntr/tests/test_gis.py
@@ -0,0 +1,249 @@
+import sys
+import unittest
+import warnings
+import os
+from os.path import abspath, dirname, isfile, join
+
+import numpy as np
+import pandas as pd
+import networkx as nx
+import wntr
+from pandas.testing import assert_frame_equal, assert_series_equal
+
+try:
+    from shapely.geometry import LineString, Point, Polygon, shape
+    has_shapely = True
+except ModuleNotFoundError:
+    has_shapely = False
+
+try:
+    import geopandas as gpd
+    has_geopandas = True
+except ModuleNotFoundError:
+    gpd = None
+    has_geopandas = False
+    
+testdir = dirname(abspath(str(__file__)))
+datadir = join(testdir, "networks_for_testing")
+ex_datadir = join(testdir, "..", "..", "examples", "networks")
+
+
+@unittest.skipIf(not has_geopandas,
+                 "Cannot test GIS capabilities: geopandas is missing")
+class TestGIS(unittest.TestCase):
+    @classmethod
+    def setUpClass(self):
+        import wntr
+
+        self.wntr = wntr
+
+        inp_file = join(ex_datadir, "Net1.inp")
+        self.wn = self.wntr.network.WaterNetworkModel(inp_file)
+        sim = wntr.sim.EpanetSimulator(self.wn)
+        self.results = sim.run_sim()
+        self.gis_data = self.wn.to_gis()
+        
+        polygon_pts = [[(25,80), (65,80), (65,60), (25,60)],
+                       [(25,60), (80,60), (80,30), (25,30)],
+                       [(40,50), (60,65), (60,15), (40,15)]]
+        polygon_data = []
+        for i, pts in enumerate(polygon_pts):
+            geometry = Polygon(pts)
+            polygon_data.append({'name': str(i+1), 
+                                 'value': (i+1)*10,
+                                 'geometry': geometry.convex_hull})
+            
+        df = pd.DataFrame(polygon_data)
+        df.set_index("name", inplace=True)
+        self.polygons = gpd.GeoDataFrame(df, crs=None)
+
+        points = [(52,71), (75,40), (27,37)]
+        point_data = []
+        for i, pts in enumerate(points):
+            geometry = Point(pts)
+            point_data.append({'geometry': geometry})
+            
+        df = pd.DataFrame(point_data)
+        self.points = gpd.GeoDataFrame(df, crs=None)
+        
+    @classmethod
+    def tearDownClass(self):
+        pass
+    
+    def test_wn_to_gis(self):
+        # Check type
+        isinstance(self.gis_data.junctions, gpd.GeoDataFrame)
+        isinstance(self.gis_data.tanks, gpd.GeoDataFrame)
+        isinstance(self.gis_data.reservoirs, gpd.GeoDataFrame)
+        isinstance(self.gis_data.pipes, gpd.GeoDataFrame)
+        isinstance(self.gis_data.pumps, gpd.GeoDataFrame)
+        isinstance(self.gis_data.valves, gpd.GeoDataFrame)
+        
+        # Check size
+        assert self.gis_data.junctions.shape[0] == self.wn.num_junctions
+        assert self.gis_data.tanks.shape[0] == self.wn.num_tanks
+        assert self.gis_data.reservoirs.shape[0] == self.wn.num_reservoirs
+        assert self.gis_data.pipes.shape[0] == self.wn.num_pipes
+        assert self.gis_data.pumps.shape[0] == self.wn.num_pumps
+        #assert self.gis_data.valves.shape[0] == self.wn.num_valves
+        
+        # Check minimal set of attributes
+        assert set(['node_type', 'elevation', 'geometry']).issubset(self.gis_data.junctions.columns)
+        assert set(['node_type', 'elevation', 'geometry']).issubset(self.gis_data.tanks.columns)
+        assert set(['node_type', 'geometry']).issubset(self.gis_data.reservoirs.columns)
+        assert set(['link_type', 'start_node_name', 'end_node_name', 'geometry']).issubset(self.gis_data.pipes.columns)
+        assert set(['link_type', 'start_node_name', 'end_node_name', 'geometry']).issubset(self.gis_data.pumps.columns)
+        #assert set(['link_type', 'start_node_name', 'end_node_name', 'geometry']).issubset(self.gis_data.valves.columns) # Net1 has no valves
+        
+    def test_gis_to_wn(self):
+        
+        wn2 = wntr.network.io.from_gis(self.gis_data)
+        G1 = self.wn.get_graph()
+        G2 = wn2.get_graph()
+        
+        assert nx.is_isomorphic(G1, G2)
+                         
+    def test_intersect_points_with_polygons(self):
+        
+        stats = wntr.gis.intersect(self.gis_data.junctions, self.polygons, 'value')
+        
+        # Junction 22 intersects poly2 val=20, intersects poly3 val=30
+        # weighted mean = (1*20+0.5*30)/2 = 17.5
+        expected = pd.Series({'intersections': ['2','3'], 'values': [20,30]})
+        expected['n'] = len(expected['values'])
+        expected['sum'] = float(sum(expected['values']))
+        expected['min'] = float(min(expected['values']))
+        expected['max'] = float(max(expected['values']))
+        expected['mean'] = expected['sum']/expected['n']
+        expected = expected.reindex(stats.columns)
+        
+        assert_series_equal(stats.loc['22',:], expected, check_dtype=False, check_names=False)
+        
+        # Junction 31: no intersections
+        expected = pd.Series({'intersections': [], 'values': [], 'n': 0, 
+                              'sum': np.nan, 'min': np.nan, 'max': np.nan,
+                              'mean': np.nan, })
+        
+        assert_series_equal(stats.loc['31',:], expected, check_dtype=False, check_names=False)
+        
+    def test_intersect_lines_with_polygons(self):
+        
+        bv = 0
+        stats = wntr.gis.intersect(self.gis_data.pipes, self.polygons, 'value', True, bv)
+
+        ax = self.polygons.plot(column='value', alpha=0.5)
+        ax = wntr.graphics.plot_network(self.wn, ax=ax)
+        
+        # Pipe 22 intersects poly2 100%, val=20, intersects poly3 50%, val=30
+        expected_weighted_mean = (20*1+30*0.5)/1.5
+        expected = pd.Series({'intersections': ['3','2'], 'values': [30,20], 'weighted_mean': expected_weighted_mean})
+        expected['n'] = len(expected['values'])
+        expected['sum'] = float(sum(expected['values']))
+        expected['min'] = float(min(expected['values']))
+        expected['max'] = float(max(expected['values']))
+        expected['mean'] = expected['sum']/expected['n']
+        expected = expected.reindex(stats.columns)
+        
+        assert_series_equal(stats.loc['22',:], expected, check_dtype=False, check_names=False)
+        
+        # Pipe 31: no intersections
+        expected = pd.Series({'intersections': ['BACKGROUND'], 'values': [bv],
+                              'n': 1, 'sum': bv, 'min': bv, 'max': bv,
+                              'mean': bv, 'weighted_mean': bv})
+        
+        assert_series_equal(stats.loc['31',:], expected, check_dtype=False, check_names=False)
+        
+        # Pipe 122
+        self.assertEqual(stats.loc['122','intersections'], ['BACKGROUND', '3', '2'])
+        # total length = 30
+        expected_weighted_mean = (bv*(5/30) + 30*(25/30) + 20*(10/30))/(40/30)
+        self.assertAlmostEqual(stats.loc['122','weighted_mean'], expected_weighted_mean, 2)
+        
+    
+    def test_intersect_polygons_with_lines(self):
+        
+        stats = wntr.gis.intersect(self.polygons, self.gis_data.pipes)
+        
+        expected = pd.DataFrame([{'intersections': ['10', '111', '11', '110', '112', '12'], 'n': 6},
+                                 {'intersections': ['111', '112', '121', '21', '122', '22', '113'], 'n': 7},
+                                 {'intersections': ['112', '21', '122', '22'], 'n': 4}])
+        expected.index=['1', '2', '3']
+        
+        assert_frame_equal(stats, expected, check_dtype=False)
+        
+    def test_intersect_polygons_with_lines_zero_length(self):
+        
+        wn = wntr.morph.break_pipe(self.wn, '22', '22_A', '22_A','22_B', split_at_point=0.5)
+        wn.add_pipe('22_0', '22_A', '22_B', length=0)
+        
+        gis_data = wn.to_gis()
+        assert gis_data.pipes.length['22_0'] == 0
+        
+        # No value
+        stats = wntr.gis.intersect(gis_data.pipes, self.polygons)
+        assert stats.shape == (14,2)
+        
+        # With value
+        stats = wntr.gis.intersect(gis_data.pipes, self.polygons, 'value')
+        assert stats.shape == (14,8)
+        
+        assert stats.loc['22_0', 'weighted_mean'] == stats.loc['22_0', 'mean'] # zero length pipe
+        assert stats.loc['22_A', 'weighted_mean'] == 20 # overlaps with 20 and 30, but zero length with 30
+        assert stats.loc['22', 'weighted_mean'] == 25 # overlaps with 20 and 30 across the entire pipe
+        
+        assert (stats.loc[stats['n']>0, 'weighted_mean'] >= stats.loc[stats['n']>0, 'min']).all()
+        assert (stats.loc[stats['n']>0, 'weighted_mean'] <= stats.loc[stats['n']>0, 'max']).all()
+        
+        
+    def test_add_attributes_and_write(self):
+        
+        self.gis_data.add_node_attributes(self.results.node['pressure'].loc[3600,:], 'Pressure_1hr')
+        self.gis_data.add_link_attributes(self.results.link['flowrate'].loc[3600,:], 'Flowrate_1hr')
+       
+        assert 'Pressure_1hr' in self.gis_data.junctions.columns
+        assert 'Pressure_1hr' in self.gis_data.tanks.columns
+        assert 'Pressure_1hr' in self.gis_data.reservoirs.columns
+        assert 'Flowrate_1hr' in self.gis_data.pipes.columns
+        assert 'Flowrate_1hr' in self.gis_data.pumps.columns
+        assert 'Flowrate_1hr' not in self.gis_data.valves.columns # Net1 has no valves
+    
+    def test_write(self):
+        prefix = 'temp_Net1'
+        components = ['junctions', 'tanks', 'reservoirs', 'pipes', 'pumps', 'valves']
+        for component in components:
+            filename = abspath(join(testdir, prefix+'_'+component+'.geojson'))
+            if isfile(filename):
+                os.remove(filename)
+            
+        self.gis_data.write(abspath(join(testdir, prefix)))
+
+        for component in components:
+            if component == 'valves':
+                continue # Net1 has no valves
+            filename = abspath(join(testdir, prefix+'_'+component+'.geojson'))
+            self.assertTrue(isfile(filename))
+
+    def test_snap_points_to_points(self):
+        
+        snapped_points = wntr.gis.snap(self.points, self.gis_data.junctions, tolerance=5.0)        
+        # distance = np.sqrt(2)*2, 5, np.sqrt(2)*3
+        expected = pd.DataFrame([{'node': '12', 'snap_distance': 2.23607, 'geometry': Point([50.0,70.0])},
+                                 {'node': '23', 'snap_distance': 5.0,      'geometry': Point([70.0,40.0])},
+                                 {'node': '21', 'snap_distance': 4.242641, 'geometry': Point([30.0,40.0])}])
+        
+        assert_frame_equal(pd.DataFrame(snapped_points), expected, check_dtype=False)
+        
+
+    def test_snap_points_to_lines(self):
+        
+        snapped_points = wntr.gis.snap(self.points, self.gis_data.pipes, tolerance=5.0)
+        
+        # distance = 1,5,3
+        expected = pd.DataFrame([{'link': '12', 'node': '12', 'snap_distance': 1, 'line_position': 0.1, 'geometry': Point([52.0,70.0])},
+                                 {'link':  '22', 'node': '23', 'snap_distance': 5.0, 'line_position': 1.0, 'geometry': Point([70.0,40.0])},
+                                 {'link': '121', 'node': '21', 'snap_distance': 3.0, 'line_position': 0.1, 'geometry': Point([30.0,37.0])}])
+        
+        assert_frame_equal(pd.DataFrame(snapped_points), expected, check_dtype=False)
+
+if __name__ == "__main__":
+    unittest.main()
\ No newline at end of file
diff --git a/wntr/tests/test_morph.py b/wntr/tests/test_morph.py
index d431f07f2..e9876ca6b 100644
--- a/wntr/tests/test_morph.py
+++ b/wntr/tests/test_morph.py
@@ -207,7 +207,7 @@ def test_skeletonize(self):
 
             # Write skel_wn to an inp file, read it back in, then extract the demands
             skel_inp_file = "temp.inp"
-            skel_wn.write_inpfile(skel_inp_file, "GPM")
+            wntr.network.write_inpfile(skel_wn, skel_inp_file, "GPM")
             skel_wn_io = wntr.network.WaterNetworkModel(skel_inp_file)
             demand_io = wntr.metrics.expected_demand(skel_wn_io)
             total_demand_io = demand_io.loc[0, :].sum()
diff --git a/wntr/tests/test_network.py b/wntr/tests/test_network.py
index ba716487b..2273cefed 100644
--- a/wntr/tests/test_network.py
+++ b/wntr/tests/test_network.py
@@ -939,7 +939,7 @@ def test_convert_controls_to_rules(self):
                 )
 
 
-class TestNetworkDict(unittest.TestCase):
+class TestNetworkIO_Dict(unittest.TestCase):
     @classmethod
     def setUpClass(self):
         import wntr
@@ -961,20 +961,122 @@ def test_dict_roundtrip(self):
             assert(wn._compare(B))
 
     def test_json_roundtrip(self):
-        import json
         for inp_file in self.inp_files:
             wn = self.wntr.network.WaterNetworkModel(inp_file)
             wn.convert_controls_to_rules()
-            wn.write_json('temp.json')
+            self.wntr.network.write_json(wn, 'temp.json')
             B = self.wntr.network.read_json('temp.json')
             assert(wn._compare(B))
 
     def test_json_pattern_dump(self):
         wn = wntr.network.WaterNetworkModel()
         wn.add_pattern('pat0', [0,1,0,1,0,1,0])
-        wn.write_json(f'temp.json')
+        self.wntr.network.write_json(wn, f'temp.json')
         
+class TestNetworkIO_GIS(unittest.TestCase):
+    @classmethod
+    def setUpClass(self):
+        import wntr
+
+        self.wntr = wntr
+
+        inp_file = join(ex_datadir, "Net6.inp")
+        self.inp_files = [join(ex_datadir, f) for f in ["Net1.inp", "Net2.inp", "Net3.inp", "Net6.inp"]]
 
+    @classmethod
+    def tearDownClass(self):
+        pass
+    
+    def test_gis_additional_attributes(self):
+        wn = self.wntr.network.WaterNetworkModel(self.inp_files[0])
+        A = wn.to_gis(pumps_as_points=True, valves_as_points=True)
+        A.junctions['new'] = range(wn.num_junctions)
+        A.junctions['pressure_exponent'] = 0.45
+        A.tanks['new'] = 4
+        A.reservoirs['new'] = 5
+        A.pipes['new'] = 6
+        A.pumps['new'] = 7
+        A.valves['new'] = 8
+        B = self.wntr.network.from_gis(A)
+
+        node_attr = B.query_node_attribute('new')
+        assert((node_attr[wn.junction_name_list] == range(wn.num_junctions)).all())
+        assert((node_attr[wn.tank_name_list] == 4).all())
+        assert((node_attr[wn.reservoir_name_list] == 5).all())
+        
+        link_attr = B.query_link_attribute('new')
+        assert((link_attr[wn.pipe_name_list] == 6).all())
+        assert((link_attr[wn.pump_name_list] == 7).all())
+        assert((link_attr[wn.valve_name_list] == 8).all())
+        
+        pressure_exponent = B.query_node_attribute('pressure_exponent')
+        assert((pressure_exponent[wn.junction_name_list] == 0.45).all())
+        
+    def test_gis_roundtrip(self):
+        for inp_file in self.inp_files:
+            wn = self.wntr.network.WaterNetworkModel(inp_file)
+            A = wn.to_gis()
+            B = self.wntr.network.from_gis(A)
+            assert(wn._compare(B, level=0))
+    
+    def test_gis_point_roundtrip(self):
+        for inp_file in self.inp_files:
+            wn = self.wntr.network.WaterNetworkModel(inp_file)
+            A = wn.to_gis(pumps_as_points=True, valves_as_points=True)
+            B = self.wntr.network.from_gis(A)
+            assert(wn._compare(B, level=0))
+    
+    def test_gis_append_roundtrip(self):
+        for inp_file in self.inp_files:
+            wn = self.wntr.network.WaterNetworkModel(inp_file)
+            A = wn.to_gis(pumps_as_points=True, valves_as_points=True)
+            B = self.wntr.network.from_gis({'junctions': A.junctions})
+            B = self.wntr.network.from_gis({'tanks': A.tanks,
+                                            'reservoirs': A.reservoirs}, B)
+            B = self.wntr.network.from_gis({'pipes': A.pipes}, B)
+            B = self.wntr.network.from_gis({'pumps': A.pumps,
+                                            'valves': A.valves}, B)
+            assert(wn._compare(B, level=0))
+    
+    def test_geojson_roundtrip(self):
+        for inp_file in self.inp_files:
+            wn = self.wntr.network.WaterNetworkModel(inp_file)
+            wntr.network.write_geojson(wn, 'temp', pumps_as_points=True, valves_as_points=False)
+            files = {}
+            if wn.num_junctions > 0:
+                files['junctions'] = 'temp_junctions.geojson'
+            if wn.num_tanks > 0:
+                files['tanks'] = 'temp_tanks.geojson'
+            if wn.num_reservoirs > 0:
+                files['reservoirs'] = 'temp_reservoirs.geojson'
+            if wn.num_pipes > 0:
+                files['pipes'] = 'temp_pipes.geojson'
+            if wn.num_pumps > 0:
+                files['pumps'] = 'temp_pumps.geojson'
+            if wn.num_valves > 0:
+                files['valves'] = 'temp_valves.geojson'
+            B = self.wntr.network.read_geojson(files)
+            assert(wn._compare(B, level=0))
+
+    def test_shapefile_roundtrip(self):
+        for inp_file in self.inp_files:
+            wn = self.wntr.network.WaterNetworkModel(inp_file)
+            wntr.network.write_shapefile(wn, 'temp')
+            files = {}
+            if wn.num_junctions > 0:
+                files['junctions'] = 'temp_junctions'
+            if wn.num_tanks > 0:
+                files['tanks'] = 'temp_tanks'
+            if wn.num_reservoirs > 0:
+                files['reservoirs'] = 'temp_reservoirs'
+            if wn.num_tanks > 0:
+                files['pipes'] = 'temp_pipes'
+            if wn.num_pumps > 0:
+                files['pumps'] = 'temp_pumps'
+            if wn.num_valves > 0:
+                files['valves'] = 'temp_valves'
+            B = self.wntr.network.read_shapefile(files)
+            assert(wn._compare(B, level=0))
 
 if __name__ == "__main__":
     unittest.main()
diff --git a/wntr/tests/test_network_minor_loss.py b/wntr/tests/test_network_minor_loss.py
index 2e0e83ee1..d684938fc 100644
--- a/wntr/tests/test_network_minor_loss.py
+++ b/wntr/tests/test_network_minor_loss.py
@@ -23,7 +23,7 @@ def test_pipe_minor_loss(self):
         sim = wntr.sim.WNTRSimulator(wn)
 
         results1 = sim.run_sim()
-        wn.write_inpfile("temp.inp", "CMH")
+        wntr.network.write_inpfile(wn, "temp.inp", "CMH")
 
         wn2 = wntr.network.WaterNetworkModel("temp.inp")
         sim = wntr.sim.EpanetSimulator(wn2)
diff --git a/wntr/tests/test_sim_reset_conditions.py b/wntr/tests/test_sim_reset_conditions.py
index f7ec656b0..ad5d37a73 100644
--- a/wntr/tests/test_sim_reset_conditions.py
+++ b/wntr/tests/test_sim_reset_conditions.py
@@ -72,7 +72,7 @@ def test_wntrsim_continue(self):
             results2.node["pressure"].plot(ax=ax)
 
             wn.options.time.duration = 100 * 3600
-            wn.write_inpfile("temp_Net3_epanet.inp")
+            wntr.network.write_inpfile(wn, "temp_Net3_epanet.inp")
 
             results3 = sim.run_sim()
             results3.node["pressure"].plot(ax=ax, title="WNTRSim_continue")