Closed line draws as edge of no length; breaks graph.py at add_edges_from()

[hasecbinusr]

Attempted to read a shapefile of OSM road data with the following code:

G=nx.read_shp('Test/Roads_Only_NoUnclass.shp')
pos = {k: v for k,v in enumerate(G.nodes())}
X=nx.Graph() #Empty graph
X.add_nodes_from(pos.keys()) #Add nodes preserving coordinates
l=[set(x) for x in G.edges()] #To speed things up in case of large objects
edg=[tuple(k for k,v in pos.items() if v in sl) for sl in l] #Map the G.edges start and endpoints onto pos
nx.draw_networkx_nodes(X,pos,node_size=1,node_color='r')
X.add_edges_from(edg)
nx.draw_networkx_edges(X,pos)
plt.xlim(1399000, 1471300) #This changes and is problem specific
plt.ylim(883800, 947800) #This changes and is problem specific
plt.xlabel('Longitude')
plt.ylabel('Latitude')
plt.title('From shapefiles to NetworkX')

This resulted in the following error:

---------------------------------------------------------------------------
NetworkXError                             Traceback (most recent call last)
<ipython-input-33-fd234537f6bd> in <module>()
      6 edg=[tuple(k for k,v in pos.items() if v in sl) for sl in l] #Map the G.edges start and endpoints onto pos
      7 nx.draw_networkx_nodes(X,pos,node_size=1,node_color='r')
----> 8 X.add_edges_from(edg)
      9 nx.draw_networkx_edges(X,pos)
     10 plt.xlim(1399000, 1471300) #This changes and is problem specific

...\networkx\classes\graph.py in add_edges_from(self, ebunch, attr_dict, **attr)
    863             else:
    864                 raise NetworkXError(
--> 865                     "Edge tuple %s must be a 2-tuple or 3-tuple." % (e,))
    866             if u not in self.node:
    867                 self.adj[u] = self.adjlist_dict_factory()

NetworkXError: Edge tuple (946,) must be a 2-tuple or 3-tuple.

Some error tracing revealed that closed lines have the same start point and end point, resulting in a 1-tuple in G.edges():

...
{(1406536.6150999963, 889657.0859999992), (1406483.8554999977, 889674.7773000002)}
{(1456536.862300001, 889682.3059)}  # Red node above
{(1457006.9228999987, 889683.0931999981), (1456907.8184999973, 889655.3469999991)}
...

1-tuple is created in read_shp() at line 94:

93                    e1, e2, attr = edge
94                    net.add_edge(e1, e2)
95                    net[e1][e2].update(attr)

Perhaps the following change to ignore these rings?

                    e1, e2, attr = edge
+                   if e1 != e2:  # make sure we don't create a null-length edge
+                       net.add_edge(e1, e2)
+                       net[e1][e2].update(attr)
-                   net.add_edge(e1, e2)
-                   net[e1][e2].update(attr)

[hasecbinusr]

Or to allow the user to choose:

+def read_shp(path, simplify=True, geom_attrs=True, ignore_rings=True):
-def read_shp(path, simplify=True, geom_attrs=True):
    print(simplify)
    """Generates a networkx.DiGraph from shapefiles. Point geometries are
    translated into nodes, lines into edges. Coordinate tuples are used as
    keys. Attributes are preserved, line geometries are simplified into start
    and end coordinates. Accepts a single shapefile or directory of many
    shapefiles.
    "The Esri Shapefile or simply a shapefile is a popular geospatial vector
    data format for geographic information systems software [1]_."
    Parameters
    ----------
    path : file or string
       File, directory, or filename to read.
    simplify:  bool
        If True, simplify line geometries to start and end coordinates.
        If False, and line feature geometry has multiple segments, the
        non-geometric attributes for that feature will be repeated for each
        edge comprising that feature.
    geom_attrs: bool
        If True, include the Wkb, Wkt and Json geometry attributes with
        each edge.
        NOTE:  if these attributes are available, write_shp will use them
        to write the geometry.  If nodes store the underlying coordinates for
        the edge geometry as well (as they do when they are read via
        this method) and they change, your geomety will be out of sync.
+   ignore_rings:  bool
+       If True, ignores non-intersecting rings that would create an edge of no length
    Returns
    -------
    G : NetworkX graph
    Examples
    --------
    >>> G=nx.read_shp('test.shp') # doctest: +SKIP
    References
    ----------
    .. [1] http://en.wikipedia.org/wiki/Shapefile
    """
    try:
        from osgeo import ogr
    except ImportError:
        raise ImportError("read_shp requires OGR: http://www.gdal.org/")

    if not isinstance(path, str):
        return

    net = nx.DiGraph()
    shp = ogr.Open(path)
    for lyr in shp:
        fields = [x.GetName() for x in lyr.schema]
        for f in lyr:
            flddata = [f.GetField(f.GetFieldIndex(x)) for x in fields]
            g = f.geometry()
            attributes = dict(zip(fields, flddata))
            attributes["ShpName"] = lyr.GetName()
            # Note:  Using layer level geometry type
            if g.GetGeometryType() == ogr.wkbPoint:
                net.add_node((g.GetPoint_2D(0)), attributes)
            elif g.GetGeometryType() in (ogr.wkbLineString,
                                         ogr.wkbMultiLineString):
                for edge in edges_from_line(g, attributes, simplify,
                                            geom_attrs):
                    e1, e2, attr = edge
+                   if e1 != e2 and ignore_rings == True:
+                       net.add_edge(e1, e2)
+                       net[e1][e2].update(attr)
+                   else:
+                       net.add_edge(e1, e2)
+                       net[e1][e2].update(attr)
-                   net.add_edge(e1, e2)
-                   net[e1][e2].update(attr)
            else:
                raise ImportError("GeometryType {} not supported".
                                  format(g.GetGeometryType()))

    return net

[hasecbinusr]

Perhaps a better option would be to add in functionality to handle this without removing the rings? In def edges_from_lines():

    if geom.GetGeometryType() == ogr.wkbLineString:
        if simplify:
            edge_attrs = attrs.copy()
+           sp = geom.GetPoint_2D(0)
            last = geom.GetPointCount() - 1
+           ep = geom.GetPoint_2D(last)
            if geom_attrs:
                edge_attrs["Wkb"] = geom.ExportToWkb()
                edge_attrs["Wkt"] = geom.ExportToWkt()
                edge_attrs["Json"] = geom.ExportToJson()
+           if sp == ep:
+               middleish = int(geom.GetPointCount() / 2)
+               mp = geom.GetPoint_2D(middleish)
+               yield (sp, mp, edge_attr)
+               yield (mp, ep, edge_attr)
+           else:
+               yield (sp, ep, edge_attrs)
-           yield (geom.GetPoint_2D(0), geom.GetPoint_2D(last), edge_attrs)

Thoughts?

[hagberg]

I like adding the rings as in your last example.

[hasecbinusr]

One problem with the last option of adding the rings is that it obliterates the ring shape and creates two edges that are essentially the same in reverse. Would it be better to convert the ring into a node or into a triangle?

EDIT: Or even for that matter, just add the edge once with yield (sp, mp, edge_attrs) and ignore the backtrace from mp -> ep?

[hagberg]

I'm not familiar with the use of these data as graphs. Which do you think is the best (or least surprising) solution?

[hasecbinusr]

A ring at the end of a line typically represents a cul-de-sac, but not always. A line of short enough length would be a cul-de-sac, but a longer line would represent a circle drive.

In the following example rendered by openstreetmap.org, the example above is converted into a cul-de-sac point as a line cap, as seen at the end of Shire Horse Blvd. Another example where the ring is not a cul-de-sac is found at Heatherwood Dr, although this is not a great example because Heatherwood Dr is a 'C' shape closed by Wickersham.

Imagining this were an actual ring drive, it would be unexpected and inappropriate to completely remove the road, just as it would be wrong to turn it into a line. A triangle may be better and a diamond perhaps the best simplification.

But for the cul-de-sacs, removing them or replacing them with a node would be expected. Essentially, they only represent the driveable surface of a road network anyway, as shown in this illustration:

In a simplified network rendering, it would be fine to remove cul-de-sacs. This introduces the issue of determining how long a ring must be before it is no longer appropriate to remove it and instead to render it as a simplified ring.

Setting simplify=False would be a better option for road networks, but this error would still be thrown by a network with closed lines.

[hagberg]

Thanks for explaining this in detail. You are far more expert me and probably most others here. So let's fix it the way you think makes the most sense.

[dschult]

I'd love to see this turned into a PR so we could include the improvements.
The time is right for getting it into v2.2. Otherwise we'll bump it to v2.3.

[dschult]

The original posted error is casued by turning an edge into a set.
The comment line there says this is to speed up the "in" in case of large objects.
But edges have at most 2 nodes and most importantly a self-loop edge becomes a single element set when you turn it into a set. That's why the X.add_edges_from creates an error.

The solution is to leave the edge as a tuple.

The second issue discussed here is more related to drawing self-loops in networkx.
This is subtle and tricky and in some sense beyond the current state of the drawing package.
nxviz may have solved drawing self-loops. I'm not sure.

I'm going to close this issue.