Are searchcursors brutally slow? They need not be

Question

Background... exploring the ArcGIS module and looking at how it worked which lead me to SpatialDataFrames, the cloaked version of a Pandas df with geometry stuff, which lead me to some of its inner workings etc etc. I use Numpy and lot and the FeatureClassToNumPy array with some fancy array splitting generally satisfies all my needs to get the geometry yanked out of a featureclass... do the magic... then send it back to a featureclass. I tend not to use searchcursors, but the ArcGIS code uses it in several places, so, I thought... give it a shot. So... helper function ... import sys from textwrap import dedent import numpy as np import arcpy from arcgis.geometry import _types from arcgis.features._data.geodataset import SpatialDataFrame as SDF import pandas as pd def fc_info(in_fc): """Return basic featureclass information, including... : SR - spatial reference object (use SR.name to get the name) : shp_fld - field name which contains the geometry object : oid_fld - the object index/id field name : - others: 'areaFieldName', 'baseName', 'catalogPath','featureType', : 'fields', 'hasOID', 'hasM', 'hasZ', 'path' : - all_flds =[i.name for i in desc['fields']] """ desc = arcpy.da.Describe(in_fc) args = ['shapeFieldName', 'OIDFieldName', 'spatialReference', 'shapeType'] shp_fld, oid_fld, SR, shp_type = [desc[i] for i in args] return shp_fld, oid_fld, SR, shp_type‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍ @time_deco def to_arr(in_fc, use_geo=False): """Convert a featureclass to a structured or recarray using a searchcursor. : :Requires: import arcpy, numpy as np :-------- : in_fc - featureclass : use_geo - True .__geo_interface__ : - list comprehension :get the row information : cycle through all geometries and get xy pairs : :References: :---------- : - see the polygon, polyline etc classes in : C:\ArcPro\Resources\ArcPy\arcpy\arcobjects\geometry.py """ shp_fld, oid_fld, SR, shp_type = fc_info(in_fc) # get the base information flds = arcpy.ListFields(in_fc) fnames = [f.name for f in flds if f.type not in ['OID', 'Geometry']] geom_flds = ['SHAPE@', oid_fld] + fnames flds = [shp_fld, oid_fld] + fnames vals = [] geoms = [] coords = [] idx = flds.index(shp_fld) with arcpy.da.SearchCursor(in_fc, field_names=geom_flds, where_clause=None, spatial_reference=SR, explode_to_points=False, sql_clause=(None, None)) as rows: for row in rows: row = list(row) geom = row.pop(idx) vals.append(row) geoms.append(geom) # pop the geometry out if use_geo: xy = geom.__geo_interface__['coordinates'] else: xy = [np.array([(pt.X, pt.Y) for pt in arr if pt]) for arr in geom] # if pt else None coords.append(np.asarray(xy)) # maybe dump the last as np.asarray del row, geom, xy del rows return vals, coords, geoms‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍ Now on to the next '2' #@time_deco def two_arrays(in_fc, keep_geom=True, dtype0=True, as_recarray=False): """Send to a numpy structured/array and split it into a geometry array : and an attribute array. : :Requires: :-------- : fc_info(in_fc) - function needed to return fc properties : keep_geom : - True, split the points into their geometry groups as an object array : - False, a sequential array of shape = (N,) : dtype0 : - True dt = [('Idx', ' shp_fld, oid_fld, SR, shp_type = fc_info(in_fc) all_flds = arcpy.ListFields(in_fc) b_flds = [i.name for i in all_flds] a = arcpy.da.FeatureClassToNumPyArray(in_fc, field_names=[oid_fld, shp_fld], explode_to_points=True, spatial_reference=SR) if dtype0: dt = [('Idx', '), ('Xs', '), ('Ys', ')] # option 1 dtb = [('Idx', '), ('Xc', '), ('Yc', ')] else: dt = [('Idx', '), ('XY', ', (2,))] dtb = [('Idx', '), ('XYc', ', (2,))] a.dtype = dt oid_fld = 'Idx' if keep_geom: a = np.split(a, np.where(np.diff(a[oid_fld]))[0] + 1) # oid_fld a = np.r_[a] if as_recarray: a = a.view(object, np.recarray) # # b = arcpy.da.FeatureClassToNumPyArray(in_fc, field_names=b_flds, explode_to_points=False, spatial_reference=SR) dtb.extend(b.dtype.descr[2:]) b.dtype = dtb return a, b‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍ Now testing on trivial samples is trivial. Before I moved on to big samples (geometries with millions of points), I tried it on a suitable one that was readily available... # ---- A multipart polygon file, with 3 parts and 0, 1 or 2 subparts a.size # = 3 a.shape # = (3,) sum([i.size for i in a]) # = 1,814,562 points‍‍‍‍‍ # ---- Big snips because of array print options ----- array([ array([(1, -65.55527496337874, 49.25722122192411), (1, -65.55305480957014, 49.25694274902389), ..., # big snip (1, -55.902732849121094, 51.63003158569353), (1, -55.90134429931635, 51.63085556030313)], dtype=[('Idx', '), ('Xs', '), ('Ys', ')]), array([(2, -81.8187484741211, 68.89791870117205), (2, -81.81718444824207, 68.89531707763712), ..., # other snip (2, -99.49025726318348, 51.53623199462936), (2, -99.49021148681629, 51.536178588867415)], dtype=[('Idx', '), ('Xs', '), ('Ys', ')]), array([(3, -114.51197814941384, 73.37343597412138), (3, -114.5083312988279, 73.37291717529314), ..., # and another (3, -69.89323425292969, 83.10884857177757), (3, -69.89251708984364, 83.10861206054705)], dtype=[('Idx', '), ('Xs', '), ('Ys', ')])], dtype=object) # ---- again, 1.8 ish million points ----‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍ Timing results Timing function for... to_arr Time: 4.63e+01s for 3 objects # 46 seconds_ = two_arrays(in_fc)Timing function for... two_arrays Time: 4.45e+00s for 2 objects # 4.5 seconds I kept the scripts intact for people to try out on their data. I even went as far as to strip out anything unnecessary in other versions but nothing reduced the time of the searchcursor approach. I am trying to get a handle on 'why'. I would have guessed that a similar base code would have been recycled in the searchcursor and the FeatureClassToNumPyArray functionality. I guess I am wrong OR there is some fundamental flaw in my implementation of 'to_arr' vs 'two_arrays'... so if you have a moment or two... have a look-seeDan

DanPatterson_Retired · Answer

I should re-read my notes more often, 'subclassing'? ... _as_narray def _shps_recarray(in_fc, to_pnts=True): """Convert shapes to a recarray quickly """ cur = arcpy.da.SearchCursor(in_fc, "*", explode_to_points=to_pnts) flds = cur.fields dt = cur._dtype a = cur._as_narray() return a, flds, dt‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍ Canada level 0 gdb tile 4ish million points with a bunch of attributes a, flds, dt = _shps_recarray(in_fc) Timing function for... _shps_recarray Time: 8.99e+00s for 3 objects a.shape(3890048,) flds ('OBJECTID', 'Shape', 'Shape_Length', 'Shape_Area', 'AREA_GEO', 'PERIM_GEO', 'PNT_COUNT') dt dtype([('OBJECTID', '), ('Shape', ', (2,)), ('Shape_Length', '), ('Shape_Area', '), ('AREA_GEO', '), ('PERIM_GEO', '), ('PNT_COUNT', ')])‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍ Certainly cuts out all the with arcpy.da.SearchCursor(in_fc, "*") as cursor:‍‍ stuff and you can then pull out what you either before or after. I knew there had to be direct support of numpy arrays beyond SpatialDataFrames. The observational gallery to this thread has grown but no comments, I guess I will turn into a discussion if there is no further contributions. References (so I don't forget again) UpdateCursor link Joshua Bixby's ... comment ... and my supplemental

DanPatterson_Retired · Answer

Update Seems that when they say 'SHAPE@' is an 'expensive' operation in the help files, they are being generous/cautious. I fiddled some and produced two variants of getting just the id and geometry data using 'SHAPE@X' and 'SHAPE@Y' to get the geometry. The first uses FeatureClassToNumPyArray and the second uses arcpy.da.SearchCursor #@time_deco def to_arr0(in_fc): """Just get the geometry and id field """ shp_fld, oid_fld, SR, shp_type = fc_info(in_fc) flds = ['OID@', 'SHAPE@X', 'SHAPE@Y'] a = arcpy.da.FeatureClassToNumPyArray(in_fc, field_names=flds, spatial_reference=SR, explode_to_points=True) dt = [('Idx', '), ('Xs', '), ('Ys', ')] a.dtype = dt return a‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍ And now for the searchcursor approach #@time_deco def to_arr1(in_fc, split_geom=True): """Pulls out the geometry and index values from a featureclass. Optionally : the full array can be split into an object array containing an array : of individual geometries. Note, this is only useful for poly* objects : :Requires: numpy, arcpy and fc_info :-------- : in_fc - featureclass : split-geom - True, separate arrays are created for each geometry object : :Notes: : - arcpy.da.SearchCursor( : in_fc, field_names, where_clause, spatial_reference, : explode_to_points, sql_clause=(None, None)) :------------------------------------------------------------------------- """ shp_fld, oid_fld, SR, shp_type = fc_info(in_fc) # flds = arcpy.ListFields(in_fc) # fnames = [f.name for f in flds if f.type not in ['OID', 'Geometry']] # flds = [oid_fld, shp_fld] g_flds = ['OID@', 'SHAPE@X', 'SHAPE@Y'] # 'SHAPE@', SHAPE@Z etc vals = [] with arcpy.da.SearchCursor(in_fc, g_flds, None, SR, True) as rows: for row in rows: vals.append(row) del row, rows # ---- construct the array ---- dt = [('Idx', '), ('Xs', '), ('Ys', ')] a = np.array(vals, dtype=dt) # ---- split out into an object array containing arrays of dt ---- if split_geom: a = np.split(a, np.where(np.diff(a['Idx']))[0] + 1) # oid_fld a = np.r_[a] return a‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍ You can ignore the commented out lines and the split_geom section. In theory, the functions do the same if split_geom=False The results are comparable within a coffee sip for a medium sized featureclass (too big to attach) a = to_arr0(in_fc) Timing function for... to_arr0 Time: 1.75e+00s for 1,814,562 objects a1 = to_arr1(in_fc, split_geom=False) Timing function for... to_arr1 Time: 2.18e+00s for 1,814,562 objects # ---- comparison of sequential Xs and Ys ---- np.allclose(a0['Xs'], a0['Xs']) True np.allclose(a0['Ys'], a1['Ys']) True‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍ The original two_arrays function can be modified to use either of the above, but to_arr1 allows the user to keep the geometries separate with a known dtype. This could be reshaped to comply with a SpatialDataFrame if desired. Sadly, Pandas is used as the middle glue. I will leave this open to see if anyone has any useful information to contribute. Dan