1

I'll try to be as much detailed as possible but please let me know if anything is missing:

I'm trying to get the number of lines in layerA that intersect with each line in LayerB and stores that number of intersections in a field in LayerB for each line, or row in tabular format. Both layers are huge and complicated loops may be problematic, but share your opinion please

the following is my logic:

-For each row in LayerB's Attribute table : ->Select(currentRow)->makefeature(currentRow)->selectByLocation(all lines that intersect with currentRow)->getCount(selected)->CalculateField(getCount)

will that theoretically work, please help me as I'm new to both arcGIS and python.

2

It is clear that we want a count of intersecting features from layer A, written to a field in Layer B - I agree the approach you've shown here should theoretically work, in terms of the conceptual steps you have sequenced to arrive at your answer.

As an alternative though, maybe you could run a spatial join using the INTERSECT rule, and look at the JoinCount field in the output layer:

import arcpy
from arcpy import env

env.workspace = "Your\\Working\\Folder\\Path"

LayerA = env.workspace + "\\" + 'TheFirstLayer.shp'
LayerB = env.workspace +  "\\" + 'TheSecondLayer.shp'

arcpy.SpatialJoin_analysis(LayerB, LayerA, "OutPutLayer", match_option = 'INTERSECT')

We can then check the 'Join_Count' field in "OutPutLayer" to view the count of LayerA features INTERSECT-ing LayerB.

Best Luck Yo.

Also - you mentioned you're 'new to arcpy' - arcpy is cool. As an aside, the syntax match_option = 'INTERSECT' is called supplying a 'keyword argument' - it lets you skip around the order of optional arguments when calling a tool like arcpy.SpatialJoin_analysis()

  • that will do it!, I can certainly use that, but might do it within a distance instead of intersect to avoid creating buffers prior to this ( I haven't mentioned that part of my story in the question). thank you very much – Solaire Oct 7 '14 at 20:17
  • How can I get the Count_Join into LayerA? I have tried creating a join using a unique ID for each row in LayerA, However I cant calculate fields in a joined layer, so I cant get this number into the original LayerA. any Ideas? – Solaire Oct 8 '14 at 18:22
  • Sometimes ArcGIS gives warnings/errors on attempting to do field calculations outside of an edit session, on joined tables. Maybe try starting an edit session on the target layer, and attempt the Join-then-FieldCalculate task again. Lemme know how this goes. Good luck. – Jim Oct 8 '14 at 18:49
  • I think this problem deserves a new question, the problem is that joins aren't saved on the disk and therefor don't actually exist (unlike spatial joins, u cant choose an output shape file). – Solaire Oct 8 '14 at 19:13
2

Your logic is fairly close, here is what I would do:

  1. Make feature layer on layer A and B
  2. Use update cursor on layer B, (make sure you pull in a unique id field or FID column to query against AND the field to store the layer A feature count)
  3. Loop through table and use unique id value from step 2 to select layer by attribute (new selection method)
  4. While in loop, use select layer by location (new selection method) against step 3 layer to layer A feature layer
  5. Use get count method for layer A
  6. Update layer B selected row with step 5 count value (cursor.updateRow(row) method)
  • will that look something like nmpeterson's answer? @nmpeterson – Solaire Oct 7 '14 at 19:47
  • Close, @nmpeterson is combining my step 3 and 4 in one statement, but you will need to include my step 1 (make feature layer). – artwork21 Oct 7 '14 at 19:49
  • Yes, I omitted your step 1 in my answer because the wording of the question made it sound like the data already existed as layers. – nmpeterson Oct 7 '14 at 20:34
1

Your intuition is correct, and the steps you've laid out will work. Doing this with an arcpy.da.UpdateCursor (ArcGIS 10.1 or later) will make the logical flow pretty simple and straightforward. Something like this:

with arcpy.da.UpdateCursor(layerA, ['SHAPE@', 'COUNT_FIELD']) as cursor:
    for row in cursor:
        arcpy.SelectLayerByLocation_management(layerB, 'INTERSECT', row[0])
        row[1] = int(arcpy.GetCount_management(layerB).getOutput(0))
        cursor.updateRow(row)

Depending on the size of your feature classes, this may take a long time to run, but that's going to be pretty hard to avoid no matter how you approach the problem.

  • thanks for putting my thoughts into code as that would take me hours, not a bad thing but still time consuming XP.. Im going with spatial joins after all. – Solaire Oct 7 '14 at 20:18
  • I'm wondering what would happen if I remove shape@, or used shape@xy? @nmpeterson – Solaire Oct 9 '14 at 15:20
  • The 'SHAPE@' field is needed to provide you with the full feature geometry for the intersection, and is used directly in the Select By Location step. If you used 'SHAPE@XY', I believe the intersection would be performed with the feature's centroid point (which may not actually overlap the feature itself!) -- definitely stick with 'SHAPE@'. – nmpeterson Oct 9 '14 at 15:24
  • Thanks, I have tried running this code block just for curiosity, and it took 16 hours to loop through 8000+ rows. it has the advantage over spatial joins of getting me the count in the same layer. however joins take few minutes to run. – Solaire Oct 9 '14 at 15:59
  • Ouch! If you're interested in speeding it up, I've noticed that GetCount_management() can be a lot slower than counting rows in a SearchCursor. If interested, you could try row[1] = sum((1 for r in arcpy.da.SearchCursor(layerB, ['OID@']))) instead. But the slowest part is definitely the Select By Location... – nmpeterson Oct 9 '14 at 16:03
1

your workflow would theoretically work, except that select by attribute works with a layer, so you need to run makefeaturelayer before select by attribute.

However, note that what you need can be done with the built in spatial join function.

  • thanks for the info, I wouldn't have guessed that would be a problem. – Solaire Oct 7 '14 at 20:20

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.