Dan_Patterson

Geometry ... ArcPy and NumPy... # 2

Blog Post created by Dan_Patterson Champion on Apr 10, 2019

Part 2... Deconstructing Geometry

---- On to the examples ----

The polygons that I will be using are shown to the right.

  1. A square, 5 points, first and last duplicates
  2. Donut with a Timbit inside
  3. A multipart with a donut hole in each
  4. The letter 'C'
  5. The letter 'V'

Each part is labelled at the labelpoint rather than the centroid, hence each part gets labelled.

 

 

 

I will be omitting examples that rely on json representation or the __geo_interface__ method since they don't add much to the functionality of constructing and deconstructing poly* type features.

 

---- Getting Geometry from FeatureClasses ----

---- The SearchCursor Approach ----

The functions getSR and _view_ are described at the end.  They are helper functions used to derive the spatial reference and to reshape the coordinates.  The key operative in this approach is on line 8, _as_narray, which does the conversion behind the scenes.

def cur_xy(in_fc, to_pnts=True):
    """Convert featureclass geometry (in_fc) to a simple 2D structured array
    with ID, X, Y values. Optionally convert to points, otherwise centroid.
    """

    SR = getSR(in_fc)
    flds = ['SHAPE@X', 'SHAPE@Y']
    cur = arcpy.da.SearchCursor(in_fc, flds, spatial_reference=SR,
                                explode_to_points=to_pnts)
    a = cur._as_narray()
    a = _view_(a)
    return a

---- Results ....

Essentially you have a simplified list of X, Y coordinates since the 'shp' was defined as 'SHAPE@X' and 'SHAPE@Y' with 'explode_to_points' set to True (False, returns centroids). Sadly you can't reconstruct the polygons unless you deal with the which points belong to what polygon.

array([[ 300020., 5000000.],
       [ 300010., 5000000.],
       [ 300010., 5000010.],
       ...,
       [ 300002., 5000002.],
       [ 300008., 5000002.],
       [ 300005., 5000008.]])

---- The FeatureClassToNumPyArray Approach ----

This produces the same results above, requires the same sort of inputs.  Timing shows that they are strongly related, especially since _as_narray has a fields and dtype property.

def fc_xy(in_fc):
    """Return the x,y coordinates for points in a featureclass (in_fc) using a
    data access searchcursor.
    """

    SR = getSR(in_fc)
    a = arcpy.da.FeatureClassToNumPyArray(in_fc, ['SHAPE@X', 'SHAPE@Y'],
                                          spatial_reference=SR,
                                          explode_to_points=True)
    a = _view_(a)
    return a

---- Retrieving shape objects ----

If you want to use arcpy directly because of the builtin methods, you need the contents of the 'shape' fields using SHAPE@' rather than just extracting the X and Y coordinates as in the previous example

def fc_shapes(in_fc, as_array=True):
    """Derive, arcpy geometry objects from a featureClass searchcursor.

    Parameters
    ----------
    in_fc : text
        Path to the input featureclass
    as_array: boolean
        True, return an object array of arcpy polygon objects.  False, returns
        a list.
    """

    SR = getSR(in_fc)
    with arcpy.da.SearchCursor(in_fc, 'SHAPE@', None, SR) as cursor:
        a = [row[0] for row in cursor]
    if as_array:
        return np.asarray(a)
    return a

---- Put it to work ----

polys = fc_shapes(in_fc, as_array=True)

polys
array([<Polygon object at 0x197f1bcebe0[0x197f0199968]>,
       <Polygon object at 0x197f1bcec18[0x197e9b48da0]>,
       <Polygon object at 0x197f1bceb70[0x197f005b738]>,
       <Polygon object at 0x197f1bceb38[0x197f005b5a8]>,
       <Polygon object at 0x197f1bceac8[0x197f005b760]>], dtype=object)

---- Geometry as a structured array ----

Nothing fancy, but there is an integer ID field indicating which feature a point belongs to and the coordinates.

A simpler version of the above... just an ID field and coordinates.  

def fc_xyID(in_fc, to_pnts=True):
    """Convert featureclass geometry (in_fc) to a simple 2D structured array
    with ID, X, Y values. Optionally convert to points, otherwise centroid.
    """

    SR = getSR(in_fc)
    flds = ['OID@', 'SHAPE@X', 'SHAPE@Y']
    cur = arcpy.da.SearchCursor(in_fc, flds, spatial_reference=SR,
                                explode_to_points=to_pnts)
    a = cur._as_narray()
    a.dtype = [('IDs', '<i4'), ('X_s', '<f8'), ('Y_s', '<f8')]
    return a

---- The results ----

a = fc_xyID(in_fc, to_pnts=True)

a
array([(1, 300010., 5000000.), (1, 300000., 5000000.), (1, 300000., 5000010.),
       (1, 300010., 5000010.),(1, 300010., 5000000.),
... snip
       (5, 300020., 5000010.),(5, 300022., 5000010.), (5, 300025., 5000002.),
       (5, 300028., 5000010.), (5, 300030., 5000010.),(5, 300026., 5000000.),
       (5, 300024., 5000000.), (5, 300020., 5000010.)],
      dtype=[('IDs', '<i4'), ('X_s', '<f8'), ('Y_s', '<f8')])

The polygon ID that each point belongs to is retained, however, it is replicated many times and the null points separating polygon parts is removed.

 

--------------------------------------------------------------------------

Helper functions

 

Helper Functions 
_view_   .... view structured arrays as an unstructured array

_getSR  .... spatial reference for a featureclass

def _view_(a):
    """Return a view of the array using the dtype and length
    Notes
    -----
    The is a quick function.  The expectation is that the array contains a
    uniform dtype (e.g 'f8').  For example, coordinate values in the form
    ``dtype([('X', '<f8'), ('Y', '<f8')])`` maybe with a Z.
    References
    ----------
    ``structured_to_unstructured`` in np.lib.recfunctions and its imports.
    `<https://github.com/numpy/numpy/blob/master/numpy/lib/recfunctions.py>`_.
    """

    v =  np.version.version.split('.')[1]  # version check
    if int(v) >= 16:
        from numpy.lib.recfunctions import structured_to_unstructured as stu
        return stu(a)
    else:
        names = a.dtype.names
        z = np.zeros((a.shape[0], 2), dtype=np.float)
        z[:,0] = a[names[0]]
        z[:,1] = a[names[1]]
        return z

Result...

a = np.array([(300015.  , 5000005.  ),
              (300005.  , 5000015.  ),
              (300010.49, 5000010.59)],
              dtype=[('SHAPE@X', '<f8'), ('SHAPE@Y', '<f8')])

_view_(a)

array([[ 300015.  , 5000005.  ],
       [ 300005.  , 5000015.  ],
       [ 300010.49, 5000010.59]])

getSR ....

def getSR(in_fc):
    """Return the spatial reference of a featureclass"""
    desc = arcpy.da.Describe(in_fc)
    SR = desc['spatialReference']
    return SR

 

Which you can use for basic spatial reference facts ( dir(SR) for a full list )

SR.type, SR.name, SR.factoryCode, SR.centralMeridian, SR.falseEasting

('Projected', 'NAD_1983_CSRS_MTM_9', 2951, -76.5, 304800.0)

 

Geometry in NumPy... # 1 

Outcomes