Identifying and deleting overlapping or near lines with similar running?

Question

I have shapefiles with lines that shows surface water runoff. I need to clear and simplify these shapefiles by condition that if some lines are near by other and has similar running, this lines has to be deleted.

I'm using Python 2.7, ArcMap 10.2 and Windows 7.

I have tried two solutions like:

Solution 1: Iterate through lines using ESRI arcpy module, make buffer with 10 meters distance on each line separately and then use Select By Location with Completely Contains condition (without the original line from which it was buffer made) and delete this lines. Then iterate next line. This works fine, but Select By Location was TOO slow, and because I have 244 shapefiles that needs to be cleaned and together contains 48 milions of lines, so it approx. 200 000 line per shapefile and I need script to do it quicky.

Solution 2: Like in previous solution, I use iteration through lines and separately use buffer, but then I want to avoid to use Select By Location, so I use Clip the lines with buffer and then compare the lengths of unclipped lines and clipped lines. If the lengths was the same - so it means that lines completely belong into buffer - I sum the IDs of this lines to list and delete from original shapefile (no from the clipped, it was only temp). But in this solution, there was Clip function so damn slow, so it is not usable for me.

Now I want to leave using arcpy module, because it not suitable for using in such a large data. I want to use GDAL/OGR to reconsturct a procedure of one of previous solution, but I'm really loosing in OGR algorithms. All that I founf out till now is load shapefile and iterate through lines, but I can't figureed out how to make buffer, identify nearby lines with similar running and delete them.

Pictures showing iteration on one line:

And pictures showing final result:

Answer 1

I think the speed problem you are having is looping around each feature and using geoprocessing tools inside the loop. They are not designed for that. They expect to process the whole dataset at once. So restructure your script to avoid the loop and it should complete "within the time to drink a cup of coffee". The data is in my opinion so small that it should not take any longer.

Step One: Forget shapefiles! They are so old school, are nowhere near as fast, lack a better spatial index and are just obsolete. Use a single filegeodatabase.

Make sure your database is on a local fast disk, preferably a SSD.

you might partition the data during the process to avoid running out of memory, but I find that using geometry objects in arcpy has all the spatial operators and is extremely fast.

There may be better approaches than weeding vectors. Surely they were generated from some sort of raster surface? Could you not look at the source and regenerate the lines again?

I can't think of an easy way to generate a sample vector set. Can you post some somewhere for me to write a demo?

Here is a demo script that runs in two seconds on 1000 lines:

#-------------------------------------------------------------------------------
# Name:        collapse_line3.py
# Purpose:     process to collapse nearby flow lines
#
# Author:      kimo
#
# Created:     14/02/2016
# Copyright:   (c) kimo 2016
# Licence:     Creative Commons 3.0 New Zealand
#-------------------------------------------------------------------------------
# Method:
#        Sort lines in study area by length for priority
#        Create a partition geometry of equal feature counts
#          For each partition, select features
#            Create a Near_partition OD table
#            Append to Near table
#          For each Near OD record:
#             if not already in remove set:
#                add FID to the keep set
#                and NEAR_FID lines, to remove set
#
#        Create SQL query layer of original ids and export

import arcpy
import sys
import os
import datetime
import logging

# global variables
start = datetime.datetime.now()

ws = 'd:/project/stackExchange/sample.gdb'
arcpy.env.workspace = ws
arcpy.env.overwriteOutput = True
in_fc = "contour"
fc = in_fc + '_sort'
buf_distance = 250 # to suit data
max_near_count = 5
keep_fc = in_fc + '_keep'
rem_fc = in_fc + '_remove'
near_tab = 'near'

logging.basicConfig(filename=os.path.dirname(ws)+"/contour.log",
# format='%(asctime)s %(message)s', datefmt='%m/%d/%Y %I:%M:%S %p',
filemode='w',
level=logging.DEBUG)
logging.debug(vars())

def create_near(layer):
    """ create a near table for each partition
        accumulate and discard"""
    arcpy.GenerateNearTable_analysis(
        in_features=layer,
        near_features=layer,
        out_table="in_memory/"+layer+'_tab',
        search_radius=buf_distance,
        location="NO_LOCATION",
        angle="NO_ANGLE",
        closest="ALL",
        closest_count=max_near_count,
        method="PLANAR")
    arcpy.management.Append("in_memory/"+layer+'_tab',near_tab)
    count = int(arcpy.management.GetCount("in_memory/"+layer+'_tab').getOutput(0))
    return count

def split_near_tab(near_tab,fc):
    # Iterate over the Near table to collect main and surrounding line IDs,
    # [could refine better still using a buffer on near line geometry to get shorter lines
    # contained as original method, but this limits the candidates]
    keep_set = set()
    remove_set = set()
    last_id = None
    with arcpy.da.SearchCursor(near_tab,['IN_FID','NEAR_FID','NEAR_DIST','NEAR_RANK']) as cur:
        for row in cur:
            # keep longest FID and remove all related lines
            fid = row[0]
            nearid = row[1]

            if fid not in remove_set:
                keep_set.add(fid) # only adds once to the set
                if nearid not in keep_set:
                    remove_set.add(nearid)

    # count what we got
    print "keep {} remove {} total {}".format(len(keep_set), len(remove_set),len(keep_set) + len(remove_set))

    # now just make a selection set and export
    # key must be indexed, but OBJECTID is already
    keep_SQL = """OBJECTID in {}""".format(str(tuple(keep_set)))
    print "{} ... {}".format(keep_SQL[0:20],keep_SQL[-20:])
    arcpy.management.MakeFeatureLayer(fc,"keep_lay",keep_SQL)
    arcpy.management.CopyFeatures("keep_lay",keep_fc)
    # deleted for the record
    remove_SQL = """OBJECTID in {}""".format(str(tuple(remove_set)))
    print "{} ... {}".format(remove_SQL[0:20],remove_SQL[-20:])
    arcpy.management.MakeFeatureLayer(fc,"remove_lay",remove_SQL)
    arcpy.management.CopyFeatures("remove_lay",rem_fc)
    return

# -------------- main -------------------------

# sort featureclass objects by length for priority by OID
if not arcpy.Exists(fc):
    arcpy.management.Sort(in_fc,fc,[['Shape_Length','DESCENDING']])
print "source count {}".format(arcpy.management.GetCount(fc))
# generalize temp copy down to 3 points for faster proximity processing

# Process all search distances in chunks (dear old Esri does not use partitioning
# for GenerateNear which has problems with featureclasses other than points)
# create a OD matrix of distances to each adjacent feature with the tolerance limit,
# and the nearest maxcount to ensure table does not explode in size

# tile count is critical for performance, try 100K for points, 10K for polylines
arcpy.CreateCartographicPartitions_cartography(in_features=fc,
    out_features="partition",
    feature_count="10000")
g = arcpy.Geometry()
total_near = 0
tiles = arcpy.management.CopyFeatures("partition",g)
msg = "Tiles {}".format(len(tiles))
logging.debug(msg)
print msg

# create an empty file for append so first tile has accurate time
if arcpy.Exists(near_tab):
    arcpy.management.Delete(near_tab)
arcpy.management.CopyRows('near_template',near_tab)
# main loop over a geometry list, not a cursor
for tile in tiles:
    arcpy.management.MakeFeatureLayer(fc,'fc_lay')
    arcpy.management.SelectLayerByLocation('fc_lay','INTERSECT',tile)
    feature_count = int(arcpy.management.GetCount('fc_lay').getOutput(0))
    if feature_count > 0:
        start_tile = datetime.datetime.now()
        n = create_near('fc_lay')
        total_near+=n
        # time log essential to identify tuning and bottlenecks
        elapsed = datetime.datetime.now() - start_tile
        rate = elapsed.seconds / ( n / 1000000.0 )
        msg = "tile {:3} count {:5} near {:5} elapsed {} rate {:8.1f} secs per M".format(tiles.index(tile),feature_count,n,elapsed,rate)
        logging.debug(msg)
        print msg
    else:
        msg = "skip {:3}".format(tiles.index(tile))
        print msg
        logging.debug(msg)
msg = "total_near {}".format(total_near)
print msg
logging.debug(msg)

# now process combined near table in one step
split_near_tab(near_tab,fc)
#-------------------------------
msg = "Well Done {}".format(datetime.datetime.now() - start)
logging.debug(msg)
print msg
logging.shutdown()

I now have scaled up the data by using 1:50K contours for a whole country and trying to thin it to match 1:250K contours. That is a bit closer to flowlines than the rivers that I first used. The GenerateNearTable fails on a larger set so I need to use Cartographic Partitions. But the tool to partition SelectByLocation also fails! Not surprising really, nothing in the Toolbox scales well.

So the next step is to replace the partitioning with something faster... However putting in a logging module and some timing, together with switching off the GenerateNear shows that the partitioning is not the slow part. If I replace the polylines with midpoints the whole process flys along at average of 58 seconds per million records with large partitions of 100,000 features, but with polylines it slows to 28,000 seconds with partitions of 5000 records. So why not just use midpoints? Well that is a bit too simplified, but suppose the polylines were generalised (for selection only) to the ends and a midpoint? That might speed up the selection, and the final lines are still selected from the source. I have replaced the test script with the latest iteration that includes partitions.

My sample of contours are very long sinuous polylines but I notice that the flowlines are almost straight but nicely curved. Maybe they have far too many vertices for analysis. If they were heavily generalised then it may process faster. The original line can be substituted by ID at the end. If I reduced the contours to have a maximum of 500 verticies it improved the speed by 10 times, and if they could be generalised further I predict another factor of 10.

I did try using Python 64 bit and it did reduce the time by 50% but I am looking for an improvement of 1000%.

Answer 2

I believe your second solution is way slower that the first one, since Clip has an inner Selection By Location plus other geometrical operations. So I suggest some improvements to your first algorithm:

1. As per @artwork21 comment, merge the shapefiles if possible. e.g. they share the same schema. This procedure removes an extra loop over the shapefiles, hence diminishes the calls to geoprocessing functions.
2. Instead of buffering in the loop, buffer the whole shapefile outside of the loop and then loop over the buffered shapefile features.
3. While buffering, make sure to use in_memory keyword to avoid writing the temporary buffer results to disk which is I/O expensive.

Answer 3

Some things to think about.

1. If your example is any indicator, you are going to have over 5 million records when you are finished. Neither shapefiles or personal geodatabases will handle that very well, so you should use either a file geodatabase or SDE.

2. Processing 48 million records won't be quick, especially when it involves a spatial selection. Perhaps think instead about building a stable script that can run overnight and in the background, so it won't be in the way. No matter how you cut it, this is going to make your computer do lots of intensive number crunching.

3. Doing a selection on 48 million records would be very resource intensive. At the very least, this will be very slow. You do NOT want to merge all the shapefiles.

4. You can't compare a script test on 1,000 records with 48 million records. The more resources a computer has to use, the slower it is going to be.

5. If you use in_memory space, think about how much memory you are using. If you are going to be working with large datasets in memory, you will need to keep track of iterations, write to disk and clear the memory periodically.

6. Make sure your layers have spatial indexes, of course.

Getting to my actual answer, I think you are on the right track.

How are you looping through the lines, though, when you are deleting them at the same time? That might be causing some of your slowness.

Here is what I would suggest.

1. You'll need to keep the original shapefiles, so convert them into a file geodatabase to create a working copy.
2. Create an output feature class for your filtered lines.
3. Loop through the shapefiles (which are now feature classes)...
4. Set a variable for the record_count (and get count of shapefile records--now feature classes--just a figure of speech).
5. Start a while loop for record_count > 0
6. Grab the first line in the "shapefile" and buffer it.
7. Do your select by location.
8. Grab the first selected record and append it to the output feature class (use a cursor on the selected set).
9. Delete all the selected records.
10. Get a record count and re-set your record_count variable
11. Make sure you clean up everything you can at the bottom of your while loop (and end your while loop).
12. And make sure you clean up everything you can at the bottom of your "shapefile" loop (and end your shapefile loop).

You'll need to repeat the script on the output feature class if the original shapefiles were not spatially separate.

This may not be perfect, since I haven't actually written the script, but hopefully it will get you started.