topic Re: SelectLayerByLocation Alternatives? in Python Questions

SelectLayerByLocation Alternatives?

JamesCrandall — Fri, 29 Sep 2017 18:46:03 GMT

Being a bit greedy here, I know but has anyone come up with more efficient/faster alternatives to arcpy.SelectLayerByLocation_management?

I'd say for the vast majority (ie. 99%) implementing this is entirely acceptable. However, I've run into a couple of instances lately that require several of these calls to be made in order to build a list/set of results that ultimately return an expected response.

...or in another instance, there are a large amount of point features (85 million) that need to be identified as COMPLETELY_WITHIN a simple single polygon feature. In most cases it's sub 30 seconds to return the expected result.

However if these processes are to be published as synchronous GP services, the processing times can be excessive in that context.

numpy? dictionaries? pandas? Anything at all that you may have tried to exceed the performance of arcpy.SelectLayerByLocation_management???

Case / scenario: select by location on points within a simple polygon feature and return a yes or no if the count is > 0

Re: SelectLayerByLocation Alternatives?

JoshuaBixby — Fri, 29 Sep 2017 20:19:35 GMT

Even for 85 million points, 30 seconds seems long to select points that intersect a simply polygon feature. Can you elaborate a bit more on how these points are stored, file geodatabase, enterprise geodatabase??? Also, you say "simple polygon feature," how simple is it? Are we talking a square, or at least something convex? Is the polygon an ArcPy Polygon or a polygon in a feature class?

Re: SelectLayerByLocation Alternatives?

JamesCrandall — Sun, 12 Dec 2021 06:31:33 GMT

File Geodatabase (10.3)

On a network share

Spatial index applied

Registered with ArcGIS Server (10.4) site

Polygon: is quite simple maybe 25 vertices, an in_memory feature class that gets created from string input of rings, but clearly the SelectbyLocation is the culprit (timing each component):

jsonFrmt = geo_convert(feature_info)
jsonData = json.loads(jsonFrmt)
inSr = arcpy.SpatialReference(3857)

for key,value in jsonData.iteritems():
    if value == 'FeatureCollection':
        pass
    else:
        for i in value:
            for k,v in i.iteritems():
                if k == 'geometry':
                    features = []
                    features2 = []
                    feature = v['coordinates']
                    polygon = arcpy.Polygon(arcpy.Array([arcpy.Point(*coords) for coords in feature]), inSr)
                    features.append(polygon)
                else:
                    pass

fc2 = arcpy.CreateFeatureclass_management('in_memory', 'tmpFC', "POLYGON")
with arcpy.da.InsertCursor(fc2, ['SHAPE@']) as cur:
    cur.insertRow([polygon])

pnt_fc = r'\\pathtonetworkfolder\myPointFeatureClass'
fl = "PointFeatures"
arcpy.MakeFeatureLayer_management(pnt_fc, "PointFeatures")
arcpy.SelectLayerByLocation_management(fl, "COMPLETELY_WITHIN", fc2)

rowcount = int(arcpy.GetCount_management(fl).getOutput(0))
if rowcount>0:
    data.update({'HasPoints': 'yes'})
else:
    data.update({'HasPoints': 'no'})‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

Re: SelectLayerByLocation Alternatives?

DanPatterson_Retired — Fri, 29 Sep 2017 20:56:30 GMT

nothing is clearly the culprit until your profile your code. You can do this with actual code profiling tools or even simple decorators (I will look for my blog post on this)

Re: SelectLayerByLocation Alternatives?

JamesCrandall — Sun, 12 Dec 2021 06:31:36 GMT

I timed the:

arcpy.CreateFeatureclass_management (.02 secs)

arcpy.MakeFeatureLayer_management (2 secs)

arcpy.SelectLayerByLocation_management (30 secs)

I'm implementing a class:

class Timer:
    def __enter__(self):
        self.start = time.clock()
        return self
    def __exit__(self, *args):
        self.end = time.clock()
        self.interval = self.end - self.start‍‍‍‍‍‍‍

Then just messaging the print statement to see results for each arcpy method as such:

t1 = Timer()
with t1:
   fc2 = arcpy.CreateFeatureclass_management('in_memory', 'tmpFC', "POLYGON")
   with arcpy.da.InsertCursor(fc2, ['SHAPE@']) as cur:
      cur.insertRow([polygon])

msgt1 = "%s...CreateFeatureClass time %.02f seconds" % (bndy,t1.interval)‍‍‍‍‍‍‍

Re: SelectLayerByLocation Alternatives?

JamesCrandall — Fri, 29 Sep 2017 21:17:19 GMT

You made an interesting comment about just using the polygon feature in the SelectbyLocation rather than creating the in_memory feature class. This works too and I can take out that extra step of CreateFeatureClass_management and shave processing time right there.

Thanks!

Still interested in what might be wrong with the SelectbyLocation

Re: SelectLayerByLocation Alternatives?

JoshuaBixby — Fri, 29 Sep 2017 21:17:27 GMT

We have had some situations where network-share file geodatabases perform terribly with selections or geoprocessing operations. Just to test the impact of storing the data on a network share, can you copy the data local and see if that speeds it up much? Also, what about a different spatial operators, some perform better than others. If you can live with "WITHIN" instead of "COMPLETELY_WITHIN", then it opens the door to try "INTERSECT."

Re: SelectLayerByLocation Alternatives?

DanPatterson_Retired — Sun, 12 Dec 2021 06:31:39 GMT

testing time....

"""
Script:   point_in_rect
Author:   Dan.Patterson@carleton.ca
Purpose:
  To determine whether points are within the extents of polygons.
"""
# ----10| ------20| ------30| ------40| ------50| ------60| ------70| ------80|
import numpy as np
np.set_printoptions(edgeitems=5,linewidth=80,precision=4,
                    suppress=True,threshold=20)


def time_deco(func):  # timing originally
    """timing decorator function
    :print("\n  print results inside wrapper or use <return> ... ")
    """
    import time
    from functools import wraps

    @wraps(func)
    def wrapper(*args, **kwargs):
        t0 = time.perf_counter()        # start time
        result = func(*args, **kwargs)  # ... run the function ...
        t1 = time.perf_counter()        # end time
        dt = t1 - t0
        print("\nTiming function for... {}".format(func.__name__))
        print("  Time: {: <8.2e}s for {:,} objects".format(dt, len(result)))
        return result                   # return the result of the function
        return dt                       # return delta time
    return wrapper


@time_deco
def contains(pnts, ext):
    """Performs a logical_and to find points within a box/ext
    comp = np.logical_and( great-eq LB, less RT) # logic
    inside = np.where(np.prod(comp, axis=1) == 1) # logic
    """
    L, B, R, T = ext.flatten()
    comp = np.logical_and(([L, B] <= pnts), (pnts <= [R, T]))  # < or <=
    idx = np.where(comp[:, 0] * comp[:, 1] == 1)[0]   # faster incarnation
    inside = pnts[idx]
    return inside

if __name__=="__main__":
    """make some points for testing, create and extent,
    time each option and optionally graph"""
    N = 10000000
    pnts = np.random.random_sample((N, 2)) * 5
    ext = np.array([[2., 2.],[3., 3.]])
    contains(pnts, ext)‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

N = 10,000,000 points and the number of points returned within the extent.

Timing function for... contains
Time: 2.44e-01s for 399,628 objects

N = 100,000,000

Timing function for... contains
Time: 3.02e+00s for 4,001,834 objects

N = 1, 000, 000, 000

damn... ran out of memory on my toy computer

So even though this is a point in simple rectangle test, you could use it to prune out the candidates before passing on to the SelectByLocation. Some fodder for thought

Re: SelectLayerByLocation Alternatives?

DanPatterson_Retired — Sun, 01 Oct 2017 17:35:39 GMT

Recent document I have started and will add to Point in .... Analysis

Re: SelectLayerByLocation Alternatives?

JamesCrandall — Mon, 02 Oct 2017 13:14:36 GMT

Great work Dan -- I'll take what I can and implement something now. Following the document too.

Re: SelectLayerByLocation Alternatives?

DanPatterson_Retired — Mon, 02 Oct 2017 21:56:58 GMT

Not going to bother with the Pandas option since the 'arcgis' module, which uses pandas in its spatialdataframe construct, essentially uses searchcursors (at some level). I can't imagine geopandas being any faster.

I think the problem with working with coordinates is the ability to vectorize the read/write process. cursors don't as your know.

As a side note, I tried to speed up the point-in-extent using numba and found no meaningful speed improvement until you got into the millions of points.... but a percentage of something fast isn't worth pursuing either

I have played around with full point in polygon searches (several standard approaches exist), just got to find my notes. In some cases (not many), the point-in-extent seemed to simply/prune the candidate tree before doing a full check on any particular point in polygon. But in your case, it may be useful to cut out the useless locations first, then do Sel By Att on the remainders. Report if you try it. However, going by my rule... 1 million of anything done in under a second ... is good enough, that's why there is coffee

Re: SelectLayerByLocation Alternatives?

ChrisSnyder — Mon, 02 Oct 2017 23:10:26 GMT

Depending on what you are going to do with the results, 'Generate Near Table' might be a faster solution or even 'Intersect'. Both of these of course write an output dataset, so....

Another approach might be to break the geometry of your polygon into several simpler geometries. For example, a large rectangle and then several smaller triangles.

As I recall, the ESRI algorithm uses this search hierarchy:

1. Bounding rectangle

2. Convex hull

3. Full geometry

Re: SelectLayerByLocation Alternatives?

JamesCrandall — Tue, 03 Oct 2017 12:52:15 GMT

Thank you for the idea!

The Generate Near Table was not as efficient compared to the SelectByLocation method.

Re: SelectLayerByLocation Alternatives?

JamesCrandall — Tue, 03 Oct 2017 17:32:44 GMT

It's not clear to me how I will efficiently set the "pnts" array from the point feature class containing the 85million features. arcpy.da.FeatureClassToNumPyArray is very slow.

Re: SelectLayerByLocation Alternatives?

ChrisSnyder — Tue, 03 Oct 2017 17:52:08 GMT

Like bixb0012 said, anything that you can do to get all the base datasets on the local machine (and off a network connection) would certainly be a big boost. Might there be z values in your point layer? If so, I bet dropping the z values would speed the read performance a bit. Hands down, all these web services are way faster with local data that is stripped down to the bare minimum.

A casual armchair observation that .shp format generally remains faster to read/write that FGDB...

I hate to say it, by most of ESRI's spatial search (select by location, etc. ) and overlay algorithms (union, etc.) are in fact "the best in the biz". Not all though...

Re: SelectLayerByLocation Alternatives?

JamesCrandall — Tue, 03 Oct 2017 18:01:05 GMT

Ultimately, the FGDB is local to the ArcGIS Server that it is registered to. I will check on the z-values.

Thanks for your input!

Re: SelectLayerByLocation Alternatives?

DanPatterson_Retired — Tue, 03 Oct 2017 18:31:28 GMT

Ahhh two tricks

dump attributes when working with geometry.
If you work a lot with the geometry at least save the featureclass geometry to a *.npy file if you are going to use it time and again. Loading native numpy binary data is zippy zippy.

If you examine how the 'arcgis' module handles getting fc's out to arrays (panda, spatialdataframe or whatever) it is basically a cursor... so until that is speeded up files that large are just going to be slow. I haven't been able to understand the difference since arcpy (and its cloaked dll's) are a black box.

If you dump the attributes, things get faster (as long as you keep track of an index like pandas), so read only the geometry... do the stuff... then reconnect with the attributes after.

Re: SelectLayerByLocation Alternatives?

JamesCrandall — Sun, 12 Dec 2021 06:31:41 GMT

Edit: added clip_analysis to generate the final polygons to use in SelectLayerByLocation process.

I think the best solution I've tested so far is to break up the polygon into individual polygons and then iterating over each polygon to perform the SelectbyLocation on the point features, checking to see if any features fall inside, then break out of that loop if that's the case.

Pretty straight forward:

1. CreateFishnet_management of the polygon feature.

2. FeatureToPolygon_management of the output Fishnet

3. Search cursor loop over each polygon in the fishnet

4. SelectLayerByLocation and perform test to count result.

5. If count > 0 then break and return with Boolean TRUE

What's interesting and the challenge is that the combination of rows/cols for the fishnet vs. the overall size of the input polygon geometry will determine the efficiency of the overall process. That is, if it's a large input polygon and too many rows/cols are created for the fishnet, then it will take longer to loop, but it also takes less time to SelectLayerByLocation.

There's a direct correlation between the content of the fishnet vs. the processing time of the SelectLayerByLocation. That is, if the faster it finds points within one of the fishnet polygons, the faster the overall process.

This beats any straight-up SelectLayerByLocation call against the full extent of the original point / polygon feature classes!

        arcpy.MakeFeatureLayer_management(pnt_fc, "PointFeatures")
        extent = polygon.extent

        originXY = extent.lowerLeft
        yaxisXY = extent.upperLeft
        opp_corner =extent.upperRight

        print "origin: " + str(originXY)
        print "yaxis: " + str(yaxisXY)
        print "opposite corner: " + str(opp_corner)

        #keep the cell width/height set to zero as it will use the opp_corner to figure it out
        cellWidth = 0
        cellHeight = 0

        #plug in the #rows/#cols you want the output grids to contain
        nrows = 4
        ncols = 4
        outFC = r'in_memory\outsections_line'
        outFishnet = r'in_memory\outsections' 
        outClip = r'in_memory\clipfeatures'
        arcpy.CreateFishnet_management(outFC, str(originXY), str(yaxisXY), cellWidth, cellHeight, nrows, ncols, str(opp_corner))
        arcpy.FeatureToPolygon_management(outFC, outFishnet)
        arcpy.Clip_analysis(outFishnet, polygon, outClip)

        with arcpy.da.SearchCursor(outClip, ['SHAPE@']) as outcur:
            for outrow in outcur:

                tshape = outrow[0]

                arcpy.SelectLayerByLocation_management("PointFeatures", "INTERSECT", tshape)
                rowcounter = int(arcpy.GetCount_management("PointFeatures").getOutput(0))

                if rowcounter > 0:
                    print "rowcounter: " + str(rowcounter)
                    return True
                    break

                else:
                    print "rowcounter: " + str(rowcounter)‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

Re: SelectLayerByLocation Alternatives?

ChrisSnyder — Tue, 03 Oct 2017 20:49:27 GMT

Very nice!

Maybe the 'Dice' tool would provide a more "optimal" way to break the polygon geometry up?

Re: SelectLayerByLocation Alternatives?

JamesCrandall — Wed, 04 Oct 2017 17:18:52 GMT

Thanks for the idea! I implemented it and it performed well, very close to the Fishnet approach I laid out. Unfortunately, the Dice_management requires an Advanced license level, which is not available on the ArcGIS Server this GP service will be published too.

This also has me back at the drawing board to get the output Fishnet line features converted to polygons too because the FeatureToPolygon_management I implemented in my example solution requires that same Advanced license.