Building feature classes from NumPy arrays

With ArcGIS 10.1, a NumPy array can be easily converted into a point feature class using the arcpy.da.NumPyArrayToFeatureClass function.

Notable is that other geometry types such as Polygon, Polyline or Multipoint are not supported by NumPyArrayToFeatureClass. However, all the tools needed to create other geometry types from NumPy are there. The numpy_array_to_features function below combines the new arcpy.da.InsertCursor and NumPy methods to turn a NumPy array into features.

import numpy
import arcpy


def numpy_array_to_features(in_fc, in_array, geom_fields, id_field):
    """
    Insert new features into an existing feature class (polygon,
    polyline or multipoint) based on a NumPy array.

    Parameters
    ----------
    in_fc : string
        An existing feature class to which new features will be added.

    in_array : structured NumPy array
        Array must include fields representing x and y coordinates, and
        an ID field.

    geom_fields: list of strings | string
        Field(s) representing x- and y-coordinates.
        If only a single numpy field is required (such as a field that
        has x,y coordinates included in a tuple) the field name can be
        passed in within a list or as a string.

    id_field: string
        The field that identifies how coordinates are grouped.  All
        coordinates with a common id value will be combined (in order
        of occurrence) into an output feature.
        The id_field is used in both the array and the feature class
        (i.e., the field name must exist in the feature class)

    """
    # Establish minimum number of x,y pairs to create proper geometry
    min_xy_dict = {'Polygon': 3, 'Polyline': 2, 'Multipoint': 1}
    min_xy_pairs = min_xy_dict[arcpy.Describe(in_fc).shapeType]

    if isinstance(geom_fields, list) and len(geom_fields) == 1:
        # Can't access a single field via a list later, extract the
        # only value
        geom_fields = geom_fields[0]

    with arcpy.da.InsertCursor(in_fc, ['SHAPE@', id_field]) as cursor:
        unique_array = numpy.unique(in_array[id_field])  # unique ids

        # Iterate through unique sets, get array that matches unique
        # value, convert coordinates to a list and insert via cursor.
        for unique_value in unique_array:
            a = in_array[in_array[id_field] == unique_value]
            if len(a) >= min_xy_pairs:  # skip if not enough x,y pairs
                cursor.insertRow([a[geom_fields].tolist(), unique_value])
            else:
                pass  # skip if not enough x,y pairs

    return

For example, you have an array that looks like this small one below:

>>> my_array
array([(21239.1, 51531.6, u'A'), (18648.5, 51551.1, u'A'),
       (18638.7, 54173.9, u'A'), (15989.5, 54153.5, u'B'),
       (13365.4, 51510.9, u'B'), (16004.5, 51524.3, u'B')], 
      dtype=[('SHAPE@X', '<f8'), ('SHAPE@Y', '<f8'), ('ID', '<U1')])

To convert this array into features, pass in an existing feature class, the array, a list of the x,y coordinate fields, and another field which controls how those coordinates are grouped. With the array listed above, 2 triangle-shaped polygon features will get added, one for the ID of ‘A’ and one for ‘B’.


numpy_array_to_features(polygon_fc, my_array, ['SHAPE@X', 'SHAPE@Y'], 'ID')

Also see Working with NumPy in ArcGIS.

About these ads

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s