Split into equal length features

Request came in last week for a way to split a line feature into 10 equal length line features.

The input looked like this
stream1

The accepted solution came from Dave on the team who sent this elegant and efficient solution.

in_fc = r'c:\projects\waterway.gdb\stream'
out_fc = r'c:\projects\waterway.gdb\stream10'
out_count = 10 # how many features desired

line = arcpy.da.SearchCursor(in_fc, ("SHAPE@",)).next()[0]
arcpy.CopyFeatures_management([line.segmentAlongLine(i/float(out_count), ((i+1)/float(out_count)), True) for i in range(0, out_count)], out_fc)

Which outpus a new feature class containing 10 line features like so
stream10

How does it work? The first 3 lines are self explanatory. So we will skip.

The following line is also fairly simple, what it does is get the geometry (shape@) of the first record (we only ask for next() once). The [0] is needed to get just the first value in the record (remember cursors return lists of values).


polyline = arcpy.da.SearchCursor(in_fc, ("SHAPE@",)).next()[0]

Next is where the magic happens, python list comprehension is used to turn the polyline object into a list of 10 (as per the out_count variable) equal length segments generated by the segmentAlongLine function. This list of polyline is then used as input to CopyFeatures (as per Using geometry objects with geoprocessing tools) which writes out the 10 polyline as individual features into the output feature class (out_fc).


arcpy.CopyFeatures_management([polyline.segmentAlongLine(i/float(out_count), ((i+1)/float(out_count)), True) for i in range(0, out_count)], out_fc)

EDIT: As was pointed out in the comments, the segmentAlongLine is new at 10.3.

Cheers

Advertisements

python/arcpy code from the 2014 dev summit

Hey Everyone,

it was great meeting and talking to everybody at the dev summit.

A number of python/arcpy presenters already shared their code to github repos below.

Happy coding.

merge polyline paths

Somebody on the team was having trouble tracing streams due to breaks in the streams introduced when converting raster to feature.

So here is slick little solution which takes n paths in a line features and makes a 1 path (aka 1 part) out of it.
None of the vertices are moved. Where there were two paths separated by a gap you will have one path with no gap.

To work property the paths have to be pointed in the same direction and not be converging. To deal with that would require quite a bit more logic.

As is script modified data in place, so back up your data before using.

Ranking field values

At the Esri International User Conference this year, an attendee came to the analysis island to ask “how do I create a rank field?”. They had run the Generate Near Table geoprocessing tool (see illustration of table below) and were looking for a way to further rank the distances. Ideally, the table would be updated to include a RANK field starting at ‘1’ for the smallest distance for each PATIENT and increasing sequentially based on the distance. The rank could then be used to facilitate further review and reporting. We were able to come up with the addRanks function below in a few minutes automating a key missing piece of the user’s workflow.

Original table
pre_ranking

Table after running addRanks function.
post_ranking

import arcpy

def addRanks(table, sort_fields, category_field, rank_field='RANK'):
    """Use sort_fields and category_field to apply a ranking to the table.

    Parameters:
        table: string
        sort_fields: list | tuple of strings
            The field(s) on which the table will be sorted.
        category_field: string
            All records with a common value in the category_field
            will be ranked independently.
        rank_field: string
            The new rank field name to be added.
    """

    # add rank field if it does not already exist
    if not arcpy.ListFields(table, rank_field):
        arcpy.AddField_management(table, rank_field, "SHORT")

    sort_sql = ', '.join(['ORDER BY ' + category_field] + sort_fields)
    query_fields = [category_field, rank_field] + sort_fields

    with arcpy.da.UpdateCursor(table, query_fields,
                               sql_clause=(None, sort_sql)) as cur:
        category_field_val = None
        i = 0
        for row in cur:
            if category_field_val == row[0]:
                i += 1
            else:
                category_field_val = row[0]
                i = 1
            row[1] = i
            cur.updateRow(row)

if __name__ == '__main__':
    addRanks(r'C:\Data\f.gdb\gen_near_table_patients2hosp',
             ['distance'], 'patient', 'rank')

Note: dBase (and shapefiles) does not support ORDER BY as used above by arcpy.da.UpdateCursor’s sql_clause argument.

Altering spatial references in Python

With ArcGIS 10.1, a spatial reference object can be created using a name or well-known ID (WKID).

# By name
sr = arcpy.SpatialReference('WGS 1984 UTM Zone 11N')

# By WKID
sr = arcpy.SpatialReference(32611)

However, once a spatial reference is created, many of the properties cannot be altered because they are read-only.

# Not possible
sr.centralMeridian = -110.0

Instead, if you need to change a property, you will need to take advantage of Python’s string manipulation capabilities. Since spatial reference properties can be expressed as well known strings, one solution is to export the spatial reference to a string, modify that string, and then use the altered string to create a new spatial reference.

import arcpy
import re
sr = arcpy.SpatialReference('WGS 1984 UTM Zone 11N')

# Change the central meridian.
sr.loadFromString(re.sub('PARAMETER\[\'Central_Meridian\'\,.+?]',
                         'PARAMETER\[\'Central_Meridian\',-120.0]',
                         sr.exportToString()))

References:
Well-known text representation of spatial reference systems

In addition to the documentation above, storage parameters like the coordinate domains and resolution, as well as tolerances are included in the spatial reference string.

List of Projected Coordinate system well-known IDs
List of Geographic Coordinate System well-known IDs

Splitting a line with a point

Recently I got an email where a colleague was trying to determine the distance a point on a line was from the beginning.  They were spinning their wheels looking at geometry operators and not seeing a solution.

First and foremost arcpy offers a broad broad collection of geoprocessing tools. What is perhaps less obvious is how those tools can be used directly with geometry objects. In the above scenario, using geometry as both input and output offers the solution to the problem. At a functional level, the Split Line at Points tool offers the solution. Being able to use it directly with geometry saves the hassle (and mess) of creating temporary datasets to work with.

import arcpy

def split_line_with_point(in_line, in_point, search_radius=0):
    """Splits a Polyline object with a PointGeometry object. Returns
    a list of two Polyline geometries. First line will include the
    starting point of the in_line. Curves are supported.

    Parameters:
    in_line: arcpy.Polyline
    in_point: arcpy.PointGeometry
    search_radius: string | float | int
        Use to split the lines based on its proximity to the point.
        Line will be split at the nearest location to the point.
        If search_radius is 0, the point must be on the line.
        Can be expressed as a number or linear unit.
    """

    split_lines = arcpy.management.SplitLineAtPoint(
        in_line,
        in_point,
        arcpy.Geometry(),
        search_radius)

    if len(split_lines) == 1:
        raise Exception('Line could not be split')
    else:
        return split_lines

So to answer the original question, to get the distance a point on a line is from the beginning, you could use the above function and then ask for the length property of the first Polyline returned in the list.

distance = split_line_with_point(line, point)[0].length

Inventorying data: a new approach

ArcPy list functions give you the options to list out a particular data type for a given workspace, but expanding that out to a directory tree meant cobbling together those list functions with Python’s os.walk and lots of updates to the workspace environment. It can be done (as shown here), but in my experience it is a process which is easy to get wrong.

Python’s os.walk is noteworthy and useful, because it does all of this, but is limited to file types. It can’t peer into a geodatabase to identify feature classes for example.

At 10.1 service pack 1, we added arcpy.da.Walk, Walk takes care of all that workspace handling for you and mimics os.walk in arguments and behaviors.

The below code wraps arcpy.da.Walk in a generator function to return a full path to all appropriate datatypes under a given workspace.

import os
import arcpy

def inventory_data(workspace, datatypes):
    """
    Generates full path names under a catalog tree for all requested
    datatype(s).

    Parameters:
    workspace: string
        The top-level workspace that will be used.
    datatypes: string | list | tuple
        Keyword(s) representing the desired datatypes. A single
        datatype can be expressed as a string, otherwise use
        a list or tuple. See arcpy.da.Walk documentation 
        for a full list.
    """
    for path, path_names, data_names in arcpy.da.Walk(
            workspace, datatype=datatypes):
        for data_name in data_names:
            yield os.path.join(path, data_name)


for feature_class in inventory_data(r"c:\data", "FeatureClass"):
    do_something(feature_class)