Adding vertices to morph.node

.gitignore

@@ -26,4 +26,3 @@ wntr/tests/performance_results
 
 examples/signal.csv
 examples/temp.hyd
-

documentation/gis.rst

@@ -773,8 +773,8 @@ to assign the demographic data, specifically the mean income, to junctions and p
     111        [0, 3]  [63326.0, 54040.0]  2  117366.0  54040.0  63326.0  58683.0      60953.558
     112        [3, 5]  [54040.0, 57620.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
+    121        [2, 3]  [91452.0, 54040.0]  2  145492.0  54040.0  91452.0  72746.0      73586.482
+    122        [2, 3]  [91452.0, 54040.0]  2  145492.0  54040.0  91452.0  72746.0      66314.037
 
 The data, water network model, and census tracts can be plotted as follows. The 
 resulting :numref:`fig-intersect-demographics` illustrates Net1 with the intersection of junctions and pipes with the census tracts (polygons). 

wntr/epanet/io.py

@@ -1270,7 +1270,7 @@ def _read_demands(self):
         has_been_read = set()
         for lnum, line in self.sections['[DEMANDS]']:
             ldata = line.split(';')
-            if len(ldata) > 1:
+            if len(ldata) > 1 and (ldata[1] != ""):
                 category = ldata[1]
             else:
                 category = None
@@ -1970,12 +1970,10 @@ def _write_vertices(self, f, wn):
         entry = '{:10s} {:20.9f} {:20.9f}\n'
         label = '{:10s} {:10s} {:10s}\n'
         f.write(label.format(';Link', 'X-Coord', 'Y-Coord').encode('latin-1'))
-        lnames = list(wn.pipe_name_list)
-        # lnames.sort()
-        for pipe_name in lnames:
-            pipe = wn.links[pipe_name]
-            for vert in pipe._vertices:
-                f.write(entry.format(pipe_name, vert[0], vert[1]).encode('latin-1'))
+        for name, link in wn.links():
+                for vert in link._vertices:
+                    f.write(entry.format(name, vert[0], vert[1]).encode('latin-1'))
+
         f.write('\n'.encode('latin-1'))
 
     def _read_labels(self):

wntr/gis/geospatial.py

@@ -206,7 +206,7 @@ def intersect(A, B, B_value=None, include_background=False, background_value=0):
         assert B_value in B.columns
     isinstance(include_background, bool)
     isinstance(background_value, (int, float))
-    assert A.crs == B.crs
+    assert A.crs == B.crs, "A and B must have the same crs."
 
     if include_background:
         background = _backgound(A, B)
@@ -224,7 +224,6 @@ def intersect(A, B, B_value=None, include_background=False, background_value=0):
 
     n = intersects.groupby('_tmp_index_name')['geometry'].count()
     B_indices = intersects.groupby('_tmp_index_name')['index_right'].apply(list)
-    B_indices.sort_values()
     stats = pd.DataFrame(index=A.index, data={'intersections': B_indices,
                                               'n': n,})
     stats['n'] = stats['n'].fillna(0)

wntr/morph/node.py

@@ -40,6 +40,11 @@ def scale_node_coordinates(wn, scale, return_copy=True):
     for name, node in wn2.nodes():
         pos = node.coordinates
         node.coordinates = (pos[0]*scale, pos[1]*scale)
+    for name, link in wn2.links():
+        if link.vertices:
+            for k in range(len(link.vertices)):
+                vertex = link.vertices[k]
+                link.vertices[k] = (vertex[0]*scale, vertex[1]*scale)
 
     return wn2
 
@@ -73,7 +78,12 @@ def translate_node_coordinates(wn, offset_x, offset_y, return_copy=True):
     for name, node in wn2.nodes():
         pos = node.coordinates
         node.coordinates = (pos[0]+offset_x, pos[1]+offset_y)
-
+    for name, link in wn2.links():
+        if link.vertices:
+            for k in range(len(link.vertices)):
+                vertex = link.vertices[k]
+                link.vertices[k] = (vertex[0]+offset_x, vertex[1]+offset_y)
+
     return wn2
 
 
@@ -107,6 +117,10 @@ def rotate_node_coordinates(wn, theta, return_copy=True):
     for name, node in wn2.nodes():
         pos = node.coordinates
         node.coordinates = tuple(np.dot(R,pos))
+    for name, link in wn2.links():
+        if link.vertices:
+            for k in range(len(link.vertices)):
+                link.vertices[k] = tuple(np.dot(R,link.vertices[k]))
 
     return wn2
 
@@ -145,6 +159,13 @@ def convert_node_coordinates_UTM_to_longlat(wn, zone_number, zone_letter, return
         lat, long = utm.to_latlon(pos[0], pos[1], zone_number, zone_letter)
         node.coordinates = (long, lat)
 
+    for name, link in wn2.links():
+        if link.vertices:
+            for k in range(len(link.vertices)):
+                vertex = link.vertices[k]
+                lat, long = utm.to_latlon(vertex[0], vertex[1], zone_number, zone_letter)
+                link.vertices[k] = (long, lat)
+
     return wn2
 
 
@@ -181,6 +202,16 @@ def convert_node_coordinates_longlat_to_UTM(wn, return_copy=True):
         easting = utm_coords[0]
         northing = utm_coords[1]
         node.coordinates = (easting, northing)
+    for name, link in wn2.links():
+        if link.vertices:
+            for k in range(len(link.vertices)):
+                vertex = link.vertices[k]
+                longitude = vertex[0]
+                latitude = vertex[1]
+                utm_coords = utm.from_latlon(latitude, longitude)
+                easting = utm_coords[0]
+                northing = utm_coords[1]
+                link.vertices[k] = (easting, northing)
 
     return wn2
 
@@ -238,14 +269,14 @@ def convert_node_coordinates_to_longlat(wn, longlat_map, return_copy=True):
 
     return wn2
 
+
 def _convert_with_map(wn, node_map, flag, return_copy):
 
     if utm is None:
         raise ImportError('utm package is required')
 
     if not len(node_map.keys()) == 2:
-        print('map must have exactly 2 entries')
-        return
+        raise Exception('map must have exactly 2 entries')
 
     if return_copy: # Get a copy of the WaterNetworkModel
         wn2 = copy.deepcopy(wn)
@@ -302,5 +333,16 @@ def _convert_with_map(wn, node_map, flag, return_copy):
             node.coordinates = (long, lat)
         elif flag == 'UTM':
             node.coordinates = (easting, northing)
+
+    for name, link in wn2.links():
+        if link.vertices:
+            for k in range(len(link.vertices)):
+                pos = link.vertices[k]
+                (easting, northing) = (np.array(pos) - cpA)*ratio + cpB
+                if flag == 'LONGLAT':
+                    lat, long = utm.to_latlon(easting, northing, zone_number, zone_letter)
+                    link.vertices[k] = (long, lat)
+                elif flag == 'UTM':
+                    node.coordinates = (easting, northing)
 
     return wn2

wntr/network/base.py

@@ -398,8 +398,11 @@ def _compare(self, other):
             return False
         if self.initial_status != other.initial_status:
             return False
-        if self.initial_setting != other.initial_setting:
-            return False
+        if (self.initial_setting is not None) ^ (other.initial_setting is not None):
+                return False 
+        elif (self.initial_setting is not None) and (other.initial_setting is not None):
+            if abs(self.initial_setting - other.initial_setting) > 1e-9:
+                return False
         if self.start_node_name != other.start_node_name:
             return False
         if self.end_node_name != other.end_node_name:

wntr/tests/test_epanet_io.py

@@ -100,8 +100,11 @@ 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.get_link("LINK-3774").vertices.append((305.31, 206.755)) # add vertex for testing
         self.wntr.network.write_inpfile(self.wn, "temp.inp", units="LPM")
-        self.wn2 = self.wntr.network.WaterNetworkModel(inp_file)
+        self.wn2 = self.wntr.network.WaterNetworkModel("temp.inp")
+        # adjusting for comparison tests
+        self.wn2.options.hydraulic.inpfile_units = 'GPM'
 
     @classmethod
     def tearDownClass(self):

wntr/tests/test_morph.py

@@ -20,65 +20,96 @@ def test_scale_node_coordinates(self):
         wn = wntr.network.WaterNetworkModel(inp_file)
         node = wn.get_node("123")
         coord = node.coordinates
+        vertex = (8.5, 27)
+        wn.get_link('10').vertices.append(vertex)
 
         wn2 = wntr.morph.scale_node_coordinates(wn, 100)
         node2 = wn2.get_node("123")
         coord2 = node2.coordinates
+        vertex2 = wn2.get_link('10').vertices[0]
 
         self.assertEqual(coord[0] * 100, coord2[0])
         self.assertEqual(coord[1] * 100, coord2[1])
+        self.assertEqual(vertex[0] * 100, vertex2[0])
+        self.assertEqual(vertex[1] * 100, vertex2[1])
 
     def test_translate_node_coordinates(self):
 
         inp_file = join(netdir, "Net3.inp")
         wn = wntr.network.WaterNetworkModel(inp_file)
         node = wn.get_node("123")
         coord = node.coordinates
+        vertex = (8.5, 27)
+        wn.get_link('10').vertices.append(vertex)
 
         wn2 = wntr.morph.translate_node_coordinates(wn, 5, 10)
         node2 = wn2.get_node("123")
         coord2 = node2.coordinates
+        vertex2 = wn2.get_link('10').vertices[0]
 
         self.assertEqual(coord[0] + 5, coord2[0])
         self.assertEqual(coord[1] + 10, coord2[1])
+        self.assertEqual(vertex[0] + 5, vertex2[0])
+        self.assertEqual(vertex[1] + 10, vertex2[1])
 
     def test_rotate_node_coordinates(self):
 
         wn = wntr.network.WaterNetworkModel()
         wn.add_junction("J1", base_demand=5, elevation=100.0, coordinates=(2, 0))
+        wn.add_junction("J2", base_demand=5, elevation=100.0, coordinates=(8, 0))
+        wn.add_pipe("P1", "J1", "J2")
+        vertex = (4, 0)
+        wn.get_link('P1').vertices.append(vertex)
 
         wn2 = wntr.morph.rotate_node_coordinates(wn, 45)
         node2 = wn2.get_node("J1")
         coord2 = node2.coordinates
+        vertex2 = wn2.get_link('P1').vertices[0]
 
         self.assertAlmostEqual(np.sqrt(2), coord2[0], 6)
         self.assertAlmostEqual(np.sqrt(2), coord2[1], 6)
+        self.assertAlmostEqual(vertex[0] * np.sqrt(2) / 2, vertex2[0])
+        self.assertAlmostEqual(vertex[0] * np.sqrt(2) / 2 , vertex2[1])
 
     def test_UTM_to_longlat_to_UTM(self):
 
         wn = wntr.network.WaterNetworkModel()
         wn.add_junction(
             "J1", base_demand=5, elevation=100.0, coordinates=(351521.07, 3886097.33)
         )  # easting, northing
+        wn.add_junction(
+            "J2", base_demand=5, elevation=100.0, coordinates=(351700.00, 3886097.33)
+        )
+        wn.add_pipe("P1", "J1", "J2")
+        vertex = (351600.00, 3886097.33)
+        wn.get_link('P1').vertices.append(vertex)
 
         wn2 = wntr.morph.convert_node_coordinates_UTM_to_longlat(wn, 13, "S")
         node2 = wn2.get_node("J1")
         coord2 = node2.coordinates
+        vertex2 = wn2.get_link('P1').vertices[0]
 
         self.assertAlmostEqual(-106.629181, coord2[0], 6)  # longitude
         self.assertAlmostEqual(35.106766, coord2[1], 6)  # latitude
+        self.assertAlmostEqual(-106.628315, vertex2[0], 6)  # longitude
+        self.assertAlmostEqual(35.106778, vertex2[1], 6)  # latitude
 
         wn3 = wntr.morph.convert_node_coordinates_longlat_to_UTM(wn2)
         node3 = wn3.get_node("J1")
         coord3 = node3.coordinates
+        vertex3 = wn3.get_link('P1').vertices[0]
 
         self.assertAlmostEqual(351521.07, coord3[0], 1)  # easting
         self.assertAlmostEqual(3886097.33, coord3[1], 1)  # northing
+        self.assertAlmostEqual(351600.00, vertex3[0], 1)  # easting
+        self.assertAlmostEqual(3886097.33, vertex3[1], 1)  # northing
 
     def test_convert_node_coordinates_to_longlat(self):
 
         inp_file = join(netdir, "Net3.inp")
         wn = wntr.network.WaterNetworkModel(inp_file)
+        vertex = (8.5, 27)
+        wn.get_link('10').vertices.append(vertex)
 
         longlat_map = {"Lake": (-106.6587, 35.0623), "219": (-106.5248, 35.191)}
         wn2 = wntr.morph.convert_node_coordinates_to_longlat(wn, longlat_map)
@@ -87,6 +118,9 @@ def test_convert_node_coordinates_to_longlat(self):
             coord = node.coordinates
             self.assertAlmostEqual(longlat_map[node_name][0], coord[0], 4)
             self.assertAlmostEqual(longlat_map[node_name][1], coord[1], 4)
+        vertex2 = wn2.get_link('10').vertices[0]
+        self.assertAlmostEqual(-106.6555, vertex2[0], 4)
+        self.assertAlmostEqual(35.0638, vertex2[1], 4)
 
         # opposite rotation
         longlat_map = {"Lake": (-106.6851, 35.1344), "219": (-106.5073, 35.0713)}
@@ -96,6 +130,9 @@ def test_convert_node_coordinates_to_longlat(self):
             coord = node.coordinates
             self.assertAlmostEqual(longlat_map[node_name][0], coord[0], 4)
             self.assertAlmostEqual(longlat_map[node_name][1], coord[1], 4)
+        vertex2 = wn2.get_link('10').vertices[0]
+        self.assertAlmostEqual(-106.6826, vertex2[0], 4)
+        self.assertAlmostEqual(35.1325, vertex2[1], 4)
 
     def test_split_pipe(self):