Dissolve with arcpy.da.UpdateCursor

5998
19
Jump to solution
09-13-2015 03:50 PM
PeterShoemark
New Contributor II

Hi, can anyone give me some guidance as to how (or if it is possible) to use arcpy.da.updatecursor to loop through multiple polylines and dissolve them one by one ? (Im on Win7, arcgis 10.2.1, Python 2.7 and the data is in an Oracle 11g sde db)

The problem I have is numerous layers that I am trying to automate removal of self overlaps but that needs to occur by objectid not by layer and I also dont want to create another dataset from the one being updated so maybe need to start an edit session first??.

While Im not a wizard at this (or I may have been able to resolve it by now), I do understand enough to be able to adjust the code below to do some other maintenance functions.

        # The "if (row[1] or 0) == 0: " code: "row[1] or 0" will return either the

        # value of the outputField, or if that value is None (null), then it will return 0...

        inputField = "SEVENTH_ID"

        outputField = "ELEVENTH_ID"

        fc = "REGION_46"

        fields = [inputField, outputField]

        with arcpy.da.UpdateCursor(fc, fields) as cursor:

            for row in cursor:

                if (row[1] or 0) == 0:

                    row[1] = row[0]

                    cursor.updateRow(row)

        del cursor, row    

I think what I am trying to do would look something like;

        # Dissolve the self overlaps by row for REGION_46

        # If the current row has no self overlaps, move onto the next row

        inputField = "OBJECTID"

        fc = "REGION_46"

        fields = [inputField]

        with arcpy.da.UpdateCursor(fc, fields) as cursor:

            for row in cursor:

                if (row[0] has self overlaps, use "DISSOLVE" to eliminate them

                    cursor.updateRow(row) then move onto the next row

        del cursor, row    

Anyone successfully acjieved this previously or can give me some useful guidance?

(Noting I have read the help menu and other resources numerous times and not been successful)

Cheers, Peter

0 Kudos
1 Solution

Accepted Solutions
JoshuaBixby
MVP Esteemed Contributor

Try this:

with arcpy.da.UpdateCursor(fc, ["SHAPE@"]) as cur: 
    for shape, in cur: 
        if shape: 
            shape = shape.union(shape) 
            cur.updateRow([shape])

The error you get is likely from Line #5 or my earlier code.  What is odd is that my earlier code worked when I tested on a feature class earlier today, but it generated the same error when I tried it tonight.  I must have tinkered enough to change something.  The code above here should work and not require removing NULLs or empty geometries.

View solution in original post

19 Replies
WesMiller
Regular Contributor III

Is this something you'll need to do more than once? A topology may be the better route.

PeterShoemark
New Contributor II

Thanks @ Wes Miller, Coincidentally the fc is a line topolgy.

I have run 'export topology errors' and can see which lines have overlap by joining the original and export tables so now I want to loop through the original lines by objectid to dissolve each lines errors.

0 Kudos
JoshuaBixby
MVP Esteemed Contributor

Can you provide an real-world example, or illustrative example, of a polyline and what the expected outcome would be?

Unfortunately, this is a case where the weaknesses of the ArcPy Geometry classes really show throw.  Other Python geometry libraries/packages like GeoDjango and Shapely offer more functionality at the geometry level.  For example, ArcPy Geometry classes don't offer a simple/is_simple functionality to tell whether a linestring is simple or complex, which would be helpful in this case at determining whether a line either crosses or overlaps itself.

I may be able to offer additional feedback if I have a specific example to work with.

0 Kudos
PeterShoemark
New Contributor II

Thnaks @ Joshua Bixby, Im not sure how to provide an example so let me expalin long hand.

I have a line which has several self overlaps on it where the previous operator has snapped 'back and forth' at times.

I want to take the line and dissolve all the self overlaps of that line without either doing the dissolve by layer or creating a new fc or layer.

Each of my dat sets has many lines within and objectid is the field i want to individually dissolve by. hth.

0 Kudos
JoshuaBixby
MVP Esteemed Contributor

Load the feature class in question, select 1 of the problem features using the ArcMap GUI, and then run the following code in the interactive Python window:

print next(arcpy.da.SearchCursor(fc, "SHAPE@"))[0].WKT

Replace fc with the name of the feature class/layer with the selected feature.  This should print out the WKT representation of the feature, which can be used to recreate the feature for testing.

PeterShoemark
New Contributor II

Thanks Joshua Bixby some sample lines below to demonstrate.

OBJECTID 1 has self ovelap and I want to dissolve it (and go to next OBJECTID)

OBJECTID 1 MULTILINESTRING ((9700614.7831249982 4432945.8557500001, 9700614.7831249982 4427945.8557500001, 9700614.7831249982 4428945.8557500001))

OBJECTID 2 has self ovelap and I want to dissolve it (and go to next OBJECTID)

OBJECTID 2 MULTILINESTRING ((9702086.5310000032 4431846.6370000001, 9702086.5308749974 4427846.6367499996, 9702086.5308749974 4428846.6367499996))

OBJECTID 3 has NO self ovelap so no requirement to dissolve it (skip and go to next OBJECTID)

OBJECTID 3 MULTILINESTRING ((9703086.5310000032 4431846.6370000001, 9703086.5308749974 4427846.6367499996, 9704269.3477500007 4427852.6163749993))

One of my fcs has more than 2000 OBJECTIDs that have self overlap and I want to dissolve them individually by OBJECTID. (and from my research perhaps arcpy.da.updatecursor is the way to go as it could cycle through the OBJECTIDs ??)  Cheers, Peter

0 Kudos
JoshuaBixby
MVP Esteemed Contributor

Thanks for providing the WKT representations, quite helpful.

A quick, or not so quick, aside before moving on.  When talking about valid geometries and spatial relationships, the nuance and complexity of the topics are often lost or oversimplified in conversation.  There are very specific definitions, and sometimes there are disagreements or discrepancies between different groups/companies/organizations when it comes to the meaning and implementation of those definitions.

Looking at the examples provided in the context of ArcPy Geometry classes:

>>> #OBJECTID 1
>>> l1 = arcpy.FromWKT('MULTILINESTRING((9700614.7831249982 4432945.8557500001,'
...                                     '9700614.7831249982 4427945.8557500001,'
...                                     '9700614.7831249982 4428945.8557500001))')
...                               
>>> l1.overlaps(l1)
False
>>> l1.crosses(l1)
False
>>> l1.touches(l1)
False
>>>
>>> #OBJECTID 2
>>> l2 = arcpy.FromWKT('MULTILINESTRING((9702086.5310000032 4431846.6370000001,'
...                                     '9702086.5308749974 4428846.6367499996,'
...                                     '9702086.5308749974 4428846.6367499996))')
...                               
>>> l2.overlaps(l2)
False
>>> l2.crosses(l2)
False
>>> l2.touches(l2)
False
>>>
>>> #OBJECTID 3
>>> l3 = arcpy.FromWKT('MULTILINESTRING((9703086.5310000032 4431846.6370000001,'
...                                     '9703086.5308749974 4427846.6367499996,'
...                                     '9704269.3477500007 4427852.6163749993))')
...                               
>>> l3.overlaps(l2)
False
>>> l3.crosses(l3)
False
>>> l3.touches(l3)
False
>>>

From one set of definitions, the examples don't overlap, cross, or even touch themselves, which makes it fairly difficult to identify the geometries that "self overlap."  The ArcPy Geometry classes, unfortunately, do not have a self-overlap method or anything similar.  For the Python implementations (Shapely, GeoDjango) of JTS/GEOS, the closest approximation to a self-overlap method for lines is is_simple:

>>> from shapely import wkt
>>> 
>>> #OBJECTID 1
>>> _l1 = wkt.loads('MULTILINESTRING((9700614.7831249982 4432945.8557500001,'
...                                  '9700614.7831249982 4427945.8557500001,'
...                                  '9700614.7831249982 4428945.8557500001))')
...                                 
>>> _l1.is_simple
False
>>> 
>>> #OBJECTID 2
>>> _l2 = wkt.loads('MULTILINESTRING((9702086.5310000032 4431846.6370000001,'
...                                 '9702086.5308749974 4428846.6367499996,'
...                                 '9702086.5308749974 4428846.6367499996))')
...                                 
>>> _l2.is_simple
True
>>> 
>>> #OBJECTID 3
>>> l3 = wkt.loads('MULTILINESTRING((9703086.5310000032 4431846.6370000001,'
...                                 '9703086.5308749974 4427846.6367499996,'
...                                 '9704269.3477500007 4427852.6163749993))')
...                                 
>>> _l3.is_simple
True
>>>

While it is clear that OBJECTID 1 is not simple because it "self overlaps," or folds back onto itself, it is less clear why OBJECTID 2 is simple.  Since JTS/GEOS libraries are OGC compliant, consecutive equal vertices are allowed in a LineString.  A basic simplify operation will remove not only redundant vertices but also some intermediate vertices if not necessary to the topology of the line.

>>> _lc = wkt.loads('LINESTRING(0 0, 1 0, 1 0, 2 0, 2 0, 3 0, 3 0, 3 0)')
>>> _lc.is_simple
True
>>> _lc.simplify(0).wkt
'LINESTRING (0 0, 3 0)'
>>>

Instead of continuing this already-too-long comment, I will wrap us this aside and post another comment with additional thoughts.

UPDATE:  OBJECTID 2 in my code above is not the same as OBJECTID 2 provided by OP.  I must have miscopied the WKT representation of OBJECTID 2 and created a line with duplicate end points.  Instead of correcting OBJECTID 2, I will leave it as is to discuss the situation when a line has consecutive vertices that are equal.

0 Kudos
JoshuaBixby
MVP Esteemed Contributor

As I alluded to in my previous comment, identifying the lines with "self overlaps" is not easy using the ArcPy Geometry classes.  One could, if so inclined, write his/her own Shamos-Hoey or Bentley–Ottmann algorithm to identify self-intersections with lines, but it might just be easier at that point to convert ArcPy geometry objects over to Shapely or GeoDjango geometry objects and run an is_simple test on them.  Another option could be generating a list of problem lines first, using Esri topology tools for example, and then using that list to only process lines that you know have issues.

I can't say it will work in all situations, but one approach worth looking into is using the union() method of the ArcPy Geometry class.

>>> #OBJECTID 1
>>> l1.length
6000.0
>>> l1u = l1.union(l1)
>>> l1u.length
5000.0
>>> l1u.WKT
u'MULTILINESTRING ((9700614.7830810547 4427945.8558959961, '
                   '9700614.7830810547 4428945.8558959961,'
                   '9700614.7830810547 4432945.8558959961))'
>>>
>>> #OBJECTID 2
>>> l2.length
3000.0002500005094
>>> l2u = l2.union(l2)
>>> l2u.length
3000.0001831054788
>>> l2u.WKT
u'MULTILINESTRING ((9702086.5311279297 4431846.6370849609,'
                   '9702086.5308837891 4428846.6369018555))'
>>>
>>> #OBJECTID 3
>>> l3.length
5182.8322396370995
>>> l3u = l3.union(l3)
>>> l3u.length
5182.832314284238
>>> l3u.WKT
u'MULTILINESTRING ((9703086.5311279297 4431846.6370849609,'
                   '9703086.5308837891 4427846.6369018555,'
                   '9704269.3479003906 4427852.6165161133))'
>>> 

As you can see, the union() method of the ArcPy Geometry class does some amount of geometry validating when it creates a new geometry.

  • For OBJECTID 1, the line was cleaned up by reordering the vertices so the line no longer folded back onto itself.
  • For OBJECTID 2, the duplicate vertex was simply removed.
  • For OBJECTID 3, the union operation effectively recreated the original line.
PeterShoemark
New Contributor II

Thank you for the comprehensive and concise response Joshua. Joshua Bixby

While I hear what you are saying, let's not get hung up on "valid geometries and spatial relationships etc etc"

Suffice to say I can open the fc attribute table and manually select each row and run the stock Dissolve tool and achieve the precise desired outcome. What I dont want to do is spend the last 6 years of my working life doing that when I suspect there are more efficient methods out there that users can point me to.

So to my original question, and in light of your findings with Union, "can anyone give me some guidance as to how (or if it is possible) to use arcpy.da.updatecursor to loop through multiple polylines and dissolve (or UNION) them one by one? (addendum, I want to use this in a stand alone Python script) Cheers, Peter

0 Kudos