2019

April 2019

# Geometry ... Reconstructing Poly Features # 4

Posted by Dan_Patterson Apr 17, 2019
Geometry

### As before, the inputs

----The polygons that I will be using are shown to the right.
1. A square, 5 points, first and last duplicates
2. A lake on an island in a lake...
3. A multipart with a two different shaped donut holes
4. The letter 'C'
5. The letter 'V'
Each part is labelled at the labelpoint rather than the centroid, hence each part gets labelled.

---- An Alternate Geometry Reconstructor ----

#### ---- Poly* to arrays revisited ----

These two scripts help blow poly features into their constituent parts.  Initially a null point (np.nan, np.nan) is used to separate the poly features.  The location of these insertions is retained and a fr_to index list is maintained so that the resultant 2D points array can be reconstructed or used in other calculations.

The outputs are: a_2d, ids, fr_to, id_fr_to which represent the 2D array, the id values of where one poly feature ends, a from-to list is constructed and a final list with the ids prepended is also included.

``def _pp_(poly):    """See poly_pnts for details.  This is for single poly feature conversion    requires         null_pnt = (np.nan, np.nan)    """    sub = []#    for i, arr in enumerate(poly):    for arr in poly:        pnts = [[pt.X, pt.Y] if pt else null_pnt for pt in arr]        sub.append(np.asarray(pnts)) #append(pnts)    return np.asarray(sub)def poly_pnts(in_fc, as_obj=False):    """Poly features to points while maintaining the location of the points in    the input features. null points are replaced with their numpy equivalent.    Parameters    ----------    in_fc : text        Full path to the featureclass.    as_obj : boolean        True returns a list of arrays of the shapes.  The array types may be        mixed.  False, returns a 2D array of points and an array of indices.    Returns    -------    False, returns a 2D array of points and the location to split    those points should you need to reconstruct the poly feature.  A `nan`    point separates the parts of the features within the poly feature, so split    first. If `as_obj` is True, a list of arrays is returned.      Notes    -----    Use `polys = fc_shapes(in_fc)` to obtain the poly features separately.    """    SR = getSR(in_fc)        polys = []    with arcpy.da.SearchCursor(in_fc, 'SHAPE@', spatial_reference=SR) as cur:        polys = [row[0] for row in cur]    id_too = []    a_2d = []    for i, p in enumerate(polys):        r = _pp_(p)                              # calls to _pp_         id_too.extend([(i, len(k)) for k in r])        a_2d.extend([j for i in r for j in i])    a_2d = np.asarray(a_2d)    id_too = np.array(id_too)    ids = id_too[:, 0]    too = np.cumsum(id_too[:, 1])    frum = np.concatenate(([0], too))    fr_to = np.array(list(zip(frum, too)))    id_fr_to = np.array(list(zip(ids, frum, too)))    if as_obj:        a_obj = [a_2d[f:t] for f,t in fr_to]  # alternate constructor        return a_obj, ids, fr_to, id_fr_to    return a_2d, ids, fr_to, id_fr_to                ``

#### A sample run

• Line 1 reads the featureclass, produces the 2D point array and the various indices
• In Line 2 and 3, I usually calculate the mean of the point cloud accounting for the null points separating the individual poly features.
• A separate object array can be constructed if you have a need to work with the individual poly features at once, otherwise, you can recreate this using the a array and the fr_to indices as shown in line 4.
• Line 5 reconstructs the original polygons so that they can be moved back into ArcGIS Pro  ... more later

Normally you would 'do stuff' between line 4 and 5 and send back the resultant array, but this is just for illustration

``a, ids, fr_to, id_fr_to = poly_pnts(in_fc)m = np.nanmin(a, axis=0)a_s = a  - ma1 = np.asarray([a[f:t] for f,t in fr_to])p_arr = [_arr_poly_(i, SR) for i in a1]  # ** to reverse np.concatenate(a1)frmt = """Polygon ids:   {}From-to pairs:{}Id_from_to array{}"""print(dedent(frmt).format(ids, fr_to, id_fr_to))Polygon ids:   [0 1 1 2 2 3 4]From-to pairs:[[ 0  5] [ 5 16] [16 26] [26 36] [36 48] [48 57] [57 65]]Id_from_to array[[ 0  0  5] [ 1  5 16] [ 1 16 26] [ 2 26 36] [ 2 36 48] [ 3 48 57] [ 4 57 65]]``

Line 15 shows that poly features 1 and 2 have 2 parts.

Lines 16-23 are the from-to pairs of points needed to reconstruct the arrays to make polygons.

Line 24- just combines the two previous lists.

#### Reconstructing the arrays to poly features

The code below does this.  A helper function, then the code block that uses a search cursor to reassemble the array to something that arcpy can use.  Sadly you have to go from a numpy array, then create points, which are placed in an arcpy Array and from there the arcpy.Arrays are assembled to form polygon or polyline features.  And Finally!!! A single-part featureclass is created, then a multipart featureclass, like the original data that went into the whole process.

Complicated? Not really, the real bottleneck is the cursors.  You need them going in (regardless of the shroud placed around them) and going out. (ie __geo_interface__,  _as_narray,  FeatureClassToNumPyArray.

I should point out here that Numpyarraytofeatureclass works with simple geometry to get poly features back, but who works with simple features all the time.

``def _arr_poly_(pnts, SR):    """Single array to polygon, array can be multipart with or without interior    rings.  Outer rings are ordered clockwise, inner rings (holes) are ordered    counterclockwise.  For polylines, there is no concept of order    Splitting is modelled after _nan_split_(arr)    """    subs = []    s = np.isnan(pnts[:, 0])    if np.any(s):        w = np.where(s)[0]        ss = np.split(pnts, w)        subs = [ss[0]]        subs.extend(i[1:] for i in ss[1:])    else:        subs.append(pnts)    aa = []    for sub in subs:        aa.append([arcpy.Point(*pairs) for pairs in sub])    poly = arcpy.Polygon(arcpy.Array(aa), SR)    return poly       def arr_poly_fc(a, p_type='POLYGON', gdb=None, fname=None, sr=None, ids=None):    """Reform poly features from the list of arrays created by poly_pnts    Parameters    ----------    a : array or list of arrays        Some can be object arrays, normally created by ``pnts_arr``    p_type : string        Uppercase geometry type    gdb : text        Geodatabase name    fname : text        Featureclass name    sr : spatial reference        name or object    ids : list/array        Identifies which feature each input belongs to.  This enables one to        account for multipart shapes.    ``_arr_poly_`` is required    """    if ids is None:        ids = np.arange(len(a)).tolist()    polys = []    for i in a:        p = _arr_poly_(i, sr)  # ---- use _arr_poly        polys.append(p)    out = list(zip(polys, ids))    name = gdb + "\\" + fname    wkspace = arcpy.env.workspace = 'in_memory'    arcpy.management.CreateFeatureclass(wkspace, fname, p_type,                                        spatial_reference=sr)    arcpy.management.AddField(fname, 'ID_arr', 'LONG')    with arcpy.da.InsertCursor(fname, ['SHAPE@','ID_arr']) as cur:        for row in out:            cur.insertRow(row)    out_fname = fname + "_mp"    arcpy.management.Dissolve(fname, out_fname, "ID_arr",                              multi_part="MULTI_PART",                              unsplit_lines="DISSOLVE_LINES")    arcpy.management.CopyFeatures(out_fname, name)    del cur    return``

#### Some Results

So just to compare the geometries, I will compare the areas calculated using arcpy and those calculated using the function below.

``polys = fc_shapes(in_fc)areas = [p.area for p in polys]areas1 = poly_area(a, ids, fr_to)areas   #   [100.0, 78.0, 155.0, 52.0, 36.0]areas1  #   array([100.,  78., 155.,  52.,  36.])``

Looks good.  The helper function employs my favorite numpy function einsum, to implement the shoelace formula.  A tad overkill for these teeny polygons, but it works blazingly fast for huge point arrays representing real world polygons.

``def poly_area(a, ids, fr_to=None):    """Calculate of a 2D array sliced into sections using an indices of the    bounds.  ``a`` is created from ``poly_pnts``    """    def _e_area(a):        """mini e_area with a twist, shoelace formula using einsum"""        x0, y1 = (a.T)[:, 1:]        x1, y0 = (a.T)[:, :-1]        e0 = np.einsum('...i,...i->...i', x0, y0)        e1 = np.einsum('...i,...i->...i', x1, y1)        return np.nansum((e0-e1)*0.5)    # ----     subs = [_e_area(a[f:t]) for f,t in fr_to]    totals = np.bincount(ids, weights=subs)    return totals``

I could calculate the area and length properties as I construct the arrays, like cursors do so that the property is readily available.  I find if I wanted those properties I would calculate them, and I don't need them floating around unused.

poly_area is only one of my helper functions for calculating poly properties.  These will be assembled at some stage and posted on GitHub and the Code Sharing site.

So why the whole exercise of converting the polygons to arrays in the first place.

Remember, the output is a 2D array of points.  Shifting, rotating, scaling, thinning, anything-ing is done on the whole array at once, not one point at a time.

Also note, that the Polygon and Polyline classes have properties and methods like intersect, union etc.  Simple functions entailing geometry are not shown/available directly from within those classes.  How do you shift a polygon by a finite amount using arcpy? Your homework.

### Coming soon

The pesky attributes

# Geometry ... Attributes actually... the other bits  # 5

Posted by Dan_Patterson Apr 17, 2019

### Part 5 ... The attributes attached to the geometry

---- As before, the inputs ----

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

1. A square, 5 points, first and last duplicates
2. A lake on an island in a lake...
3. A multipart with a two different shaped donut holes
4. The letter 'C'
5. The letter 'V'

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

---- An Alternate Geometry Reconstructor ----
---- Arrays to Poly* Features ----

Never use <null> in a table.  To many posts on the forum on how to trap them, find them, replace them.

Always put in a value to represent all conditions.  Too many people use None <null> as the catchall category.  In reality all observations need to be classified exactly, there really is no such thing as <null>.  You either made and observation or you didn't.  If you didn't, your classification scheme should provide a key indicating that.

If an observation was made but the phenomenon/parameter/whatever was actually not there, there should be a key for that.  Similar, for 'I forgot', 'Wasn't my job' or whatever other excuses exist.

• For observations recorded as floating point numbers, that truly yielded 'nothing/None/nadda/zilch', I use 'nan' (np.nan) since it is a recognized number.
• For text/string observations, I use None since None the object translates to 'None' the text easily in most tables.
• For time, there is now 'not a time' (NaT), but I don't work with time, preferring to use the string incarnations of those observations
• For integers... sadly there is no 'Nint'.  You have to provide an actual integer value to represent no value observed, although you desperately tried.   You can use the old school -999, or even 2**8, 2**16 etcetera.
• For all of the above, anything that is truly doesn't represent an observation where no value was observed, you will have to provide alternatives

#### Making none/null/real nothingness

You can add to, or remove from, the list below.  These are some that I use.  I will draw your attention to the NumPy incarnations for integers.  Equivalent values exist for floats.  Any value that ensures that you will take a second look if a calculation looks weird is good.  However if your table contains <nulls> even after my lecture above, this will help mitigate your... stupidity is such a harsh word... but you get my drift

``def _make_nulls_(in_fc, int_null=-999):    """Return null values for a list of fields objects, excluding objectid    and geometry related fields.  Throw in whatever else you want.    Parameters    ----------    in_flds : list of field objects        arcpy field objects, use arcpy.ListFields to get a list of featureclass        fields.    int_null : integer        A default to use for integer nulls since there is no ``nan`` equivalent        Other options include    >>> np.iinfo(np.int32).min # -2147483648    >>> np.iinfo(np.int16).min # -32768    >>> np.iinfo(np.int8).min  # -128    [i for i in cur.__iter__()]    [[j if j is not None else -999 for j in i ] for i in cur.__iter__() ]    """    nulls = {'Double': np.nan, 'Single': np.nan, 'Float': np.nan,             'Short': int_null, 'SmallInteger': int_null, 'Long': int_null,             'Integer': int_null, 'String':str(None), 'Text':str(None),             'Date': np.datetime64('NaT')}    #    desc = arcpy.da.Describe(in_fc)    if desc['dataType'] != 'FeatureClass':        print("Only Featureclasses are supported")        return None, None    in_flds = desc['fields']    shp = desc['shapeFieldName']    good = [f for f in in_flds if f.editable and f.name != shp]    fld_dict = {f.name: f.type for f in good}    fld_names = list(fld_dict.keys())    null_dict = {f: nulls[fld_dict[f]] for f in fld_names}    # ---- insert the OBJECTID field    return null_dict, fld_names``

#### My favorite way of getting just the attributes

Such nice functions, FeatureClassToNumPyArray, TableToNumPyArray, and back the other way.

I am sure many of you have explored where it all comes from only to find it all buried in a *.pyd file

``import arcpy.da as apdapd.__file__ # ---- 'C:\\...install path...\\Resources\\ArcPy\\arcpy\\da.py'# ---- which is actually just imports arcgisscripting# ---- so import it directlyimport arcgisscripting as agsags.__file__ # ---- 'C:\\...install path...\\bin\\Python\\envs\\arcgispro-py3\\lib\\# site-packages\\arcgisscripting.pyd'``

No bother, since you can pull out data for the attributes nicely accounting for <null> records.

``def fc_data(in_fc, verbose=False):    """Pull all editable attributes from a featureclass tables.  During the    process, <null> values are changed to an appropriate type.    Parameters    ----------    in_fc : text        Path to the input featureclass    verbose : boolean        Requires ``'prn_rec' in locals().keys()`` in order to set to ``True``.        ``prn_rec`` is imported from arraytools.frmts    """    flds = ['OID@', 'SHAPE@X', 'SHAPE@Y']    null_dict, fld_names = _make_nulls_(in_fc, int_null=-999)    fld_names = flds + fld_names    new_names = ['OID_orig', 'X_centroid', 'Y_centroid']    a = arcpy.da.FeatureClassToNumPyArray(in_fc, fld_names, skip_nulls=False,                                          null_value=null_dict)    a.dtype.names = new_names + fld_names[3:]    if verbose:        try: prn_rec(a)  # ** prn_rec imported from arraytools.frmts        except: print(a)    return np.asarray(a)``

The explorers amongst us, may have discovered a few searchcursor shortcuts

``fld_names = ['OBJECTID', 'Long_1', 'Short_1', 'Float_1', 'Double_1', 'Text_1']cur = arcpy.da.SearchCursor(in_fc, fld_names, explode_to_points=False)z = cur._as_narray()Traceback (most recent call last): File "<ipython-input-220-b4e724f0982b>", line 1, in <module> z = cur._as_narray()``

Sadly, the integer fields with <nulls> bring the whole shortcut down.

The searchcursor actually has enough information in it to create a structured/recarray.

If you have a clean table with no nulls, the actual calls to _dtype and fields show that you can clearly link cursors and NumPy arrays.  Too bad, the whole integer fix isn't incorporated, but _as_narray and FeatureClassToNumPyArray and TableToNumPyArray yield the same results on a 'clean' dataset.

``dir(cur)[...snip... '_as_narray', '_dtype', 'fields', 'next', 'reset']``

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

### Coming soon

• Work with the geometry... and/or … work with the attributes, then put it all back together.

# Geometry ... Deconstructing poly* features  # 3

Posted by Dan_Patterson Apr 10, 2019

Part 3... Alternate geometry deconstructors

---- 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.

### ---- An Alternate Geometry Deconstructor ----

#### ---- Poly* Features to Points ----

This one may seem a bit odd, so I will explain it abit.Returns

•  Line 13, gets the spatial reference for the featureclass
• A null point is created on line 16... none of this None stuff or a
• not_null_point = arcpy.Point()
not_null_point.X   #  0.0
not_null_point.Y   #  0.0
• Line 17, your familiar SearchCursor is used to pull the geometry object out for further work. (I use an enumerator since I use 'i' in a different incarnation of the script...more in a different blog)
• Line 21 pulls our arcpy Point coordinates unless 'None' is found.  If it is, then it is replaced with my null_pnt.
• Line 24 records the number of points in the shape.  This value is used later for poly* reconstruction and is returned is 'as_obj' is False.
• Line 26 produces a numpy array from the coordinates and null_pnts.
• Line 27 does a cumulative sum of the points in each shape.  This is a key... 'c' allows us to split a 2D array of points into subgroups.  A subgroup may contain a null_pnt which further allows one to split the points into parts whether they be interior rings and/or array parts.
• And the final lines.  If 'as_obj' is True a list of arrays is returned.  The arrays may be conventional 2D arrays, 3D arrays or object arrays.  Not for you to worry though.  If 'as_obj' is false, then a 2D list of points is returned as well as the 'c' array should you need to split it later.

#### Why the difference

• With the 2D array of points, it is blooming easy to determine the parameters for the whole array dataset multiple times without having to repeat things.  For example, means, std deviations.
• The 2D array allows you to do things that you can't do easily.  Every tried to rotate polygons about an origin point?  How about a tiny shift in the coordinates.  A quick densification but you only have a Standard license?  Think of something that you can't do easily with the plain arcpy geometry objects... they are easy in array format.
• With the object array option, you can pick and choose which geometry objects to operate on and you don't need to reconstruct them prior to sending them back to arcpy geometries.

====================================================================================

### The Script

``def poly_pnts(in_fc, as_obj=False):    """Poly features to points while maintaining the location of the points in    the input features. null points are replaces with their numpy equivalent.    Parameters    ----------    in_fc : text        Full path to the featureclass.    as_obj : boolean        True returns a list of arrays of the shapes.  The array types may be        mixed.  False, returns a 2D array of points and an array of indices.    """    SR = getSR(in_fc)        out = []    c = []    null_pnt = (np.nan, np.nan)  #  inf_pnt = [np.PINF, np.PINF]    with arcpy.da.SearchCursor(in_fc, 'SHAPE@', spatial_reference=SR) as cur:        for i, row in enumerate(cur):            sub = []            poly = row[0]            pnts = [[pt.X, pt.Y] if pt else null_pnt                     for arr in poly for pt in arr]            sub.extend(pnts)  # removed from next line  sub.append(inf_pnt)            c.append(len(pnts))            out.extend(sub)    out = np.asarray(out)    c = np.cumsum(c)    if as_obj:        return np.split(out, c)[:-1]  # cut off the last empty slice    return out, c``

---- Results ....

``a, c = poly_pnts(in_fc, as_obj=False)a.shape  # (58, 2) ... 58 points of X,Y coordinatesaarray([[ 300010., 5000000.],       [ 300000., 5000000.],       [ 300000., 5000010.],       [ 300010., 5000010.],       [ 300010., 5000000.],  # ---- 5th point, last of the first polygon       [ 300010., 5000020.],       [ 300010., 5000010.],       [ 300000., 5000010.],       [ 300000., 5000020.],       [ 300010., 5000020.],       [     nan,      nan],  # ---- end of the first part of the second shape    .... big snipcarray([ 5, 21, 41, 50, 58], dtype=int32)  # note '5' ``

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

Coming soon

The next blog post will go over reconstructing the geometries after processing their array representation.

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

Posted by Dan_Patterson 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)polysarray([<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)``

#### Filter Blog

By date: By tag: