# Calculate line density in multiple overlapping polygons

I am working with two layers. One is polygons (animal locations buffered by 5 km), and the other is polylines (for example, roads). I would like to calculate the line density for each polygon, separately.

I did this for multiple other buffer sizes (ex 2km) and was able to do it easily by creating a line density raster, and then using Zonal Statistics as Table to extract the mean pixel value.

However, when I run this with my 5km buffers, I get a number of errors such as "number of unique values exceeds 500 .." and "Error in executing grid expression". I'm afraid my raster cell size might be too small (ie many unique values), but I would like to continue using this raster layer, if I can.

Anyone have any other ideas as to how I can extract this information? I've also tried Extract my Mask, playing around with Line Density, and Isectpolyrast in GME, to no avail. I'd also like to avoid running this for each polygon separately, as there are hundreds of thousands, unless I can write a python script (my coding isn't up to par...).

I would probably approach this a little differently: It might be easiest to intersect your line data with the polygon features. This would assign your unique polygon ID to each line feature that falls within the polygon. You could then divide the total length of the polylines within each polygon by the area of the polygon to get a density (for example 'linear feet per acre').

• I don't think this will work for me because the polygons overlap substantially - there are multiple lines in each polygon, and lines can be in multiple polygons...I haven't gotten this approach to work yet. I could do it in a loop of sorts through python, but again my coding isn't up to it yet.
– Mel
Jan 12, 2016 at 2:05

# Approach 1

Having read your question and seen @CStardbird's response then your reply there is a simple but probably not very efficient way of doing this. Using model builder you could step through each polygon and then run the intersect tool.

Be aware that a line could be clipped out by several polygons if you have overlapping polygons. Thus a segment of line may be duplicated. Is that an issue? The line is in all the overlapping polygons so why not report X number of times?

The model will intersect the line layer with each individual polygon, so if your polygon dataset had 999 polygons you get 999 shapefiles.

• You would then run the merge tool to bring all the datasets back into a single dataset.
• Recompute the polyline lengths if your data is a shapefile.
• Run a summary statistics grouping by polygon ID and summing the polyline length
• Join the merged table to your summary table by polygon ID
• Now you can calculate line density.

# Approach 2

This approach makes use of the fact you can Union a single dataset.

1. Simplify your polygon data to just a field which is it ID, i.e. delete other fields so when you carry out the following processing it will just be easier to understand.
2. Run a Union on just your polygon data. This will cut up your polygon data into ever smaller polygons and will resolve the issue of overlaps. Polygons within overlaps will be duplicated and retain the source ID.
3. Run this chopped up polygon data and your polyline data through the Intersect tool.
4. You then have all the necessary data, polylines per polygon ID. So a polyline that occurs in multiple overlapping polygons will get duplicated as many times. Then it's just a matter of running the Summary stats tool too aggregate up your data by ID's

This second approach will be very fast compared to approach 1

• You have the right idea for the outcome I would like, and I could do this, or more efficiently I could run it through python. However there are hundreds of thousands of polygons, and this would take a really long time. I was hoping someone might have a more efficient suggestion.
– Mel
Jan 12, 2016 at 12:52
• I have suggested an alternative and significantly faster approach. Jan 12, 2016 at 14:42

I ended up using python. If anyone else comes across this, here is some example code (you just have to merge them all together at the end). Its not the most efficient, but it works.

``````print ""
print "###############################################"
print "Starting loop"
for i in inds:
i = str(i)
print ""
print "###############################################"
print "Exporting polygon" + i
whereclause = "\""+ind+"\"="+i
arcpy.MakeFeatureLayer_management(inPoly,"lyr"+i,whereclause)
Individual = arcpy.CopyFeatures_management("lyr"+i, FileOut+"Ind_"+str(i))
# Intersect each polygon with line file
print "###############################################"
print "Intersecting polygon" + i
IndIntersect = arcpy.Intersect_analysis ([Individual, inLines], FileOut+"Ind_Int_"+str(i))
#add column for the length of the lines
print "###############################################"
print "Adding length column to polygon" + i