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

DanPatterson_Retired · ‎04-10-2019

Part 2... Deconstructing Geometry

target="_self">---- Getting Geometry from FeatureClasses ----
---- The SearchCursor Approach ----
---- The FeatureClassToNumPyArray Approach ----
---- Retrieving shape objects ----
---- Geometry as a structured array ----
Helper functions
_getSR .... spatial reference for a featureclass
Other links

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

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

A square, 5 points, first and last duplicates
Donut with a Timbit inside
A multipart with a donut hole in each
The letter 'C'
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)‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍`