Latest Contributions by DanPatterson

A followup to the recent einsum post, but with a field calculator twist   # ---- a numeric field, python expression which uses the shape field to return the longest edge on the perimeter longest_seg(!Shape!) # ---- the code block ----- import json import numpy as np # def longest_seg(a): """Longest edge""" v = json.loads(a.JSON) r = v['rings'] a = np.array(r).squeeze() a = a[0] if a.ndim > 2 else a diff = a[1:] - a[:-1] ret = np.sqrt(np.einsum('ij,ij->i', diff, diff)) val = np.max(ret) return val   lines 12-13 :  really? having to convert to JSON format just to get at the array values? line 14 :  the json format assumes a multipolygon so a nested, nested set of nests of values... numpy squeeze() gets rid of all the extra nesting line 15 :  a check, but used for dimensions > 2, it will take the first part (the outer ring of the first shape) lines 16-17 :  einsum .... read the previous article line 18:  you will get an array of the segment lengths forming the perimeter... this demo is using the max, so modify the code if you want the minimum or the average or the 3rd from the end, etc etc.  Just dont' over-fluff field calculator code to try and do everything line 19  return the value   There really has to be a better way of converting an arcpy shape into another format which can be used without having to use cursors and deal with the featureclass rather than the feature itself.  like make __geo_interface__ read the shape without having to use that in cursors.   ... View more

Start with a polygons boundary.  Where are the unnecessary points? they can be points that are duplicate segments that overlap and/or reverse directions (squish a Z to a _ ... hard to see the middle 2 points ehh?) and other mishaps Can you spot the redundant points?  Original image on the left, cleaned image on the right?   edgy1.png edgy1_clean.png                 Perhaps if I make it more obvious.  A densified version of the inputs and the cleaned version. edgy1_dens_labels.png edgy1_clean.png                 Better? The code def _clean_segments_(a, tol=1e-06): """Remove overlaps and extra points on poly* segments. In `npg_geom_hlp`. Parameters ---------- a : array The input array or a bit from a Geo array. tol : float The tolerance for determining whether a point deviates from a line. Notes ----- - Segments along a straight line can overlap (a construction error). [[0,0], [5, 5], [2, 2], [7, 7]] # points out of order - Extraneous points can exist along a segment. [[0,0], [2, 2], [5, 5], [7, 7]] # extra points not needed for line. """ cr, ba, bc = _bit_crossproduct_(a, extras=True) # -- avoid duplicating the 1st point (0). whr = [i for i in np.nonzero(np.abs(cr) > tol)[0] if i != 0] vals = np.concatenate((a[0][None, :], a[whr], a[-1][None, :]), axis=0) return vals Which calls def _bit_crossproduct_(a, is_closed=True, extras=False): """Cross product. Used by `is_convex`, `_angles_3pnt_` and `_clean_segments_`. .. note:: np.cross for 2D arrays was deprecated in numpy 2.0, use `cross_2d` """ def cross2d(x, y): return x[..., 0] * y[..., 1] - x[..., 1] * y[..., 0] if is_closed: if np.allclose(a[0], a[-1]): # closed loop, remove dupl. a = a[:-1] ba = a - np.concatenate((a[-1][None, :], a[:-1]), axis=0) bc = a - np.concatenate((a[1:], a[0][None, :]), axis=0) # cr = np.cross(ba, bc) + 0.0 # deprecated cr = cross2d(ba, bc) if extras: return cr, ba, bc return cr   Now vals is a NumPy array which can be converted back to an arcpy shape using _arr_poly_ and _poly_to_array_ can be used to get the array(s) from arcpy geometry. import arcpy from arcpy import Array, Exists, Multipoint, Point, Polygon, Polyline def _arr_poly_(arr, SR, as_type): """Slice the array where nan values appear, splitting them off.""" aa = [Point(*pairs) for pairs in arr] if as_type.upper() == 'POLYGON': poly = Polygon(Array(aa), SR) elif as_type.upper() == 'POLYLINE': poly = Polyline(Array(aa), SR) return poly def _poly_to_array_(polys): """Convert polyline or polygon shapes to arrays for use in numpy. Parameters ---------- polys : tuple, list Polyline or polygons in a list/tuple """ def _p2p_(poly): """Convert a single ``poly`` shape to numpy arrays or object.""" sub = [] pt = Point() # arcpy.Point() for arr in poly: pnts = [[p.X, p.Y] for p in arr if pt] sub.append(np.asarray(pnts, dtype='O')) return sub # ---- if not isinstance(polys, (list, tuple)): polys = [polys] out = [] for poly in polys: out.extend(_p2p_(poly)) # or append, extend it is return out   Enough for now.  More on my github page numpy geometry    ... View more