1

I have a script attempting to identify Many-to-Many matches between two different parcel layers. It is an iterative process, during which I dissolve the selected polygons in layers A and B, intersect the dissolves, and compare the intersected area to the areas of both dissolves. If no match is found, I add features in B which intersect the A-dissolve, then add features in A which match the B-dissolve, etc... iteratively, until either a match is found, or no further parcels meet the criteria (currently 40% intersection with the dissolve) to be included.

This function currently takes about 8 hrs to run and, ideally, I would like to run it multiple times on decreasing match tolerance (alternating with other routines attempting 1-1, 1-Many, Many-1 matches).

Update: I've shifted from using the intersect_analysis and dissolve_management, to doing much of the work in SQL, via ArcSDESQLExecute. This seems to have sped things up significantly, but it is still too slow.

I just discovered the polygon geometry object, and its intersect method and union methods. I'm wondering if these are likely to be quicker.

I can, of course, rewrite my code, and test it. (And am working on that, as I await a response). But this looks like a major re-write, and before I open this Pandora's Box, I was wondering if anyone, with under-the-hood knowledge could tell me whether this is likely to be worth the effort?

Here is my function:

def Match_MM(tolerance):
    # get list of OIDs that have not been matched
    sql='''
    select distinct a.DOR_OBJECTID from tmp_matches a left join (select DOR_OBJECTID from tmp_matches where matchType is not null) b
    on a.DOR_OBJECTID=b.DOR_OBJECTID where b.DOR_OBJECTID is null order by a.DOR_OBJECTID
    '''  
    dorOIDs = db_conn.execute(sql)
    if isinstance(dorOIDs,int): 
        dorOIDs=[dorOIDs]
    else:
        dorOIDs=[a[0] for a in dorOIDs]

    for dorOID in dorOIDs:    
        arcpy.AddMessage(dorOID)

        dorList=[dorOID]
        pwdList=[]

        found=False
        i=1

        db_conn.execute("IF OBJECT_ID('dorDissolve', 'U') IS NOT NULL drop table dorDissolve")
        db_conn.execute('select shape into dorDissolve from MATCHING_DOR where DOR_OBJECTID={}'.format(dorOID))
        while True:
            if (i%2)==1:
                id1='DOR_OBJECTID'
                id2='PARCELID'
                fc1='MATCHING_DOR'
                fc2='MATCHING_PWD'
                dissolve1='dorDissolve'
                dissolve2='pwdDissolve'
            else:
                id1='PARCLEID'
                id2='DOR_OBJECTID'
                fc1='MATCHING_PWD'
                fc2='MATCHING_DOR'
                dissolve1='pwdDissolve'
                dissolve2='dorDissolve'

            # get features in lyr2 that intersect dissolve1
            sql="""
            select b.{}
            from {} a join {} b
            on a.shape.STIntersects(b.shape)=1
            where (a.shape.STIntersection(b.shape).STArea()>(.4*a.shape.STArea())) or
                (a.shape.STIntersection(b.shape).STArea()>(.4*b.shape.STArea()))
            """.format(id2,dissolve1,fc2)
            newList = db_conn.execute(sql)

            # put the returned OIDs into a list. if nothing was returned, exit loop
            if isinstance(newList,bool): 
                break
            elif isinstance(newList,int): 
                newList=[newList]
            else:
                newList=[a[0] for a in newList]

            # if nothing was added, exit loop
            if (i%2==1):
                if len(newList)<=len(pwdList): break
                pwdList=newList
            else:
                if len(newList)<=len(dorList): break
                dorList=newList

            # dissolve lyr2
            db_conn.execute("IF OBJECT_ID('{0}', 'U') IS NOT NULL drop table {0}".format(dissolve2))
            sql='select geometry::UnionAggregate(shape) as shape into {} from {} where {} in ({})'.format(dissolve2,fc2,id2,str(newList)[1:-1])
            db_conn.execute(sql)

            # intersect dissolve1 and dissolve2
            sql='''
            select a.shape.STArea(), b.shape.STArea(), a.shape.STIntersection(b.shape).STArea()
            from dorDissolve a join pwdDissolve b 
            on a.shape.STIntersects(b.shape)=1'''
            val = db_conn.execute(sql)

            # if the intersection covers both dissolves, greater than the tolerance, success
            if ((val[0][2]/val[0][0])>tolerance) and ((val[0][2]/val[0][1])>tolerance): 
                found=True
                break   
            i += 1

        if found:
            arcpy.AddMessage('found')
            # update the tmp_matches list
            sql="""
            update tmp_matches
            set matchType='M-M', MatchID={}
            where dor_objectid in ({}) and parcelid in ({})
            """.format(min(dorList),str(dorList)[1:-1],str(pwdList)[1:-1])
            db_conn.execute(sql)

            # remove these OIDs from further consideration 
            for dor in dorList:
                if dor in dorOIDs: dorOIDs.remove(dor)
  • As per the Tour there should be one question per question and since you are asking about ArcPy a code snippet should be included. – PolyGeo Feb 13 at 20:42
1

I am assuming using Geomerty objects would be much faster than using feature classes to do what you want.

First of all, Geomerty objects do not hold attributes, just the geometry you need to work with. Of course, if you need to work with geometries and their attributes, I suggest not using Geomerty objects.

From this page of the official Esri documentation it reads:

In many geoprocessing workflows, you may need to run a specific operation using coordinate and geometry information but don't necessarily want to go through the process of creating a new (temporary) feature class, populating the feature class with cursors, using the feature class, then deleting the temporary feature class. Geometry objects can be used instead for both input and output to make geoprocessing easier.

So I am assuming the existance of the Geometry class is meant (also) to be used whenever you only need the geometries of your features, and this implicitly (I guess) should get your job done much faster.

They might also be faster because everything is done in memory and you don't have to access any file in the file system (just a guess, I am not an IT expert).

Also, in my experience, I do see an improve in my script times by using them, so, might be worth the effort.

Finally, consider using a 64 bit Python with arcpy if you have access to it (e.g. if you have ArcGIS Server or ArcGIS Pro installed). This would improve the performance in any case respect to the 32 bit!

  • 1
    Another interesting obseration I add as a separate comment is that there exist different alternatives to arcpy.Geometry which seem faster. Refere to these two links: 1) community.esri.com/blogs/tilting/2018/12/06/… and 2) gis.stackexchange.com/a/220277/9518 – umbe1987 Feb 13 at 15:47
  • I think Intersect_analysis will often/sometimes be faster but it will depend on precisely what the user is trying to do. Without code that's hard to tell. – PolyGeo Feb 13 at 20:45
  • @umbe. thank you. i switched to using ArcSDESQLExecute, as this seemed to be quicker, and made for simpler code. but its still too slow (>8hrs per cycle). i will try the geometry object. it seems like the coding is gonna be more complicated, but i gotta get these run-times down – ericwerfel Feb 15 at 15:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.