Latest Contributions by DanPatterson

Back to the simple point example of comparing two shapes for common points where there might be floating point issues.  Note that in the b1 array, I introduces a really tiny difference in the second point so that it differs from its previous incarnation.... b a = np.array([[0., 0.], [1., 4.], [4., 3.], [5., 0.],[0., 0.]]) b = np.array([[1., 0.], [1., 4.], [5., 0.], [1., 0.]]) b1 = np.array([[1., 0.], [1., 4.000001], [5., 0.], [1., 0.]]) So if we compare a to b for equality, then we would have an issue because of the small floating point difference.  To account for this, we need to step up the check to account for this.  Here is the code.  Read the header. def a_eq_b(a, b, atol=1.0e-8, rtol=1.0e-5, return_pnts=False): r"""Return indices of, or points where, two point arrays are equal. Parameters ---------- a, b : 2d arrays No error checking so ensure that a and b are 2d, but their shape need not be the same. Notes ----- Modified from `np.isclose`, stripping out the nan and ma checks. Adjust atol and rtol appropriately. >>> np.nonzero(a_eq_b(a, b))[0] One could use:: >>> np.equal(A, B).all(1).nonzero()[0] >>> np.where((a[:, None] == b).all(-1).any(1))[0] If you aren't worried about floating-point issues in equality checks. """ a = np.atleast_2d(a) b = np.atleast_2d(b) if b.size > 2: b = b[:, None] w = np.less_equal(np.abs(a - b), atol + rtol * abs(b)).all(-1) if w.ndim > 1: w = w.any(-1).squeeze() if return_pnts: return b[w].squeeze() return w.squeeze() The magic is to account for tolerances (covered in the NumPy documentation).  This allows one to identify points that are ridiculously close that they should be considered the same.  If the difference is due to errors, then they can be fixed. So compare a to b and a to b1 with and without tolerance checks. # -- a to b comparison, both are correct a_eq_b(a, b, atol=1.0e-8, rtol=1.0e-5, return_pnts=False) array([0, 1, 1, 0]) a_eq_b(a, b, atol=0., rtol=0., return_pnts=False) array([0, 1, 1, 0]) # -- a to b1 comparison, tolerance check is correct a_eq_b(a, b1, atol=1.0e-8, rtol=1.0e-5, return_pnts=False) Out[160]: array([0, 1, 1, 0]) # -- exactitude check only identifies one 'equality' a_eq_b(a, b1, atol=0., rtol=0., return_pnts=False) Out[161]: array([0, 0, 1, 0]) So when you things are all good, make sure that you account for the small differences. If coordinates/values are the result of a calculation, it is really important to round/truncate your results to the expected precision that is associated with the rest of your data. Enough for now.   Addendum see Common points for a fuller discussion Common points ... View more

Start with two simple shapes that have common points, but this time, we are just interested in the points and not the segments. Shape 'a' is the black outline with the red point identifiers and shape 'b' has the blue outline and the green labels. There coordinate values are: a = np.array([[0., 0.], [1., 4.], [4., 3.], [5., 0.],[0., 0.]]) b = np.array([[1., 0.], [1., 4.], [5., 0.], [1., 0.]]) Being the quick people we are, it is obvious that they share points [1., 4.] and [5., 0.] which are pairs [1, 1] and [3, 2] using zero-based enumeration.   So how would you find out what points are duplicates and where they are using python and numpy.  Let's break it down. Python approach ids_ = [] out = [] for i, a_ in enumerate(b): sub = [] for j, b_ in enumerate(a): chk = (a_ == b_) sub.append(chk) if chk.all(): ids_.append((i, j)) out.append(sub) out = np.array(out) ids_ = np.array(ids_)   Nothing like 'enumerate' since you can utilize the cycle to get an index as well as a value.  Line 6 compares the value in 'a' to that in 'b' and saves the result to a subarray.  If they are both equal, then a separate list is appended with the indices from the respective arrays.  We are left with two output arrays.... 'out' and 'ids_' Now, on to NumPy, with a little side note.  Normally if you have two arrays of the same shape (ie. the number of rows and columns), you can just use a direct comparison. Arrays of equal shape # -- take the first 4 points of array a_ and the 4 points of b and compare chk0 = np.equal(a[:-1], b) # which is the same as chk1 = (a[:-1] == b) # -- both yield array([[0, 1], [1, 1], [0, 0], [0, 1]]) # # now 'where' they are all equal, this occurs is simply np.where((a[:-1] == b).all(-1)) # or, more simply np.nonzero((a[:-1] == b).all(-1)) # both yield, just one match (array([1]),)   Arrays of unequal shape # -- since a_ has 5 points and b_ has 4, we can do a direct comparison # b_ has to have another dimension added to it so you can compare 2D elments # to a 2D array compare_ = (a == b[:, None]) # -- b, copy all elements and add a newaxis all_chk_ = compare_.all(-1) whr = np.nonzero(all_chk_)   Both the Python and NumPy approaches yield the same results. # compare array([[[0, 1], [1, 0], [0, 0], [0, 1], [0, 1]], [[0, 0], [1, 1], [0, 0], [0, 0], [0, 0]], [[0, 1], [0, 0], [0, 0], [1, 1], [0, 1]], [[0, 1], [1, 0], [0, 0], [0, 1], [0, 1]]]) # all_chk_ array([[0, 0, 0, 0, 0], [0, 1, 0, 0, 0], [0, 0, 0, 1, 0], [0, 0, 0, 0, 0]]) # whr -- where the values are located in the inputs (array([1, 2]), array([1, 3])) # # and the values obtained by slicing from the original array b[whr[0]] array([[ 1.00, 4.00], [ 5.00, 0.00]]) a[whr[1]] array([[ 1.00, 4.00], [ 5.00, 0.00]])   So the approaches yield the same results but the NumPy approach is substantially faster and can simplify the code notation.  For large arrays, NumPy can be substantially faster.  Run your own tests. Create 1000 random points a1 = np.random.random(size=(1000, 2)) * 10 b1 = np.random.random(size=(1000, 2)) * 10 %%timeit using python approach and NumPy approach for unequal sizes is over 100x faster.  This factor varies with array sizes, but it is useful if speed and storage requirements are at a premium. Add the above to your toolset. ADDENDUM When things aren't quite that perfect, there is always a workaround.... see Really close points - Esri Community   ... View more

Take a few polygons.  Hexs will do. I took the liberty of labelling the segments.  The extent origin is the bottom left.  Polygons are oriented clockwise. All that needs to be done is identify which pairs of points are duplicates.  The procedure, simplified, is as follows: Take each polygon and 'segment' it. Ravel/flatten each segment, aka 'point-pair' to the form (x0, y0, x1, y1).  The first segment (0) has the values [ -3.000, -2.598, -1.500, 0.000] "simply" 😂  identify the segments that share the same coordinates.   This is where numpy comes into play since the identification can be done all at once. A listing of the coordinates and other information follows:       ID : Shape: ID by part R : ring : outer 1, inner 0 P : part : 1 or more ID R P x y ID R P x y 0 1 1 [ -3.00 -2.60] 2 1 1 [ -3.00 2.60] [ -1.50 0.00] [ -1.50 5.20] [ 1.50 0.00] [ 1.50 5.20] [ 3.00 -2.60] [ 3.00 2.60] [ 1.50 -5.20] [ 1.50 0.00] [ -1.50 -5.20] [ -1.50 0.00] [ -3.00 -2.60] [ -3.00 2.60] 1 1 1 [ 1.50 0.00] 3 1 1 [ 1.50 5.20] [ 3.00 2.60] [ 3.00 7.79] [ 6.00 2.60] [ 6.00 7.79] [ 7.50 0.00] [ 7.50 5.20] [ 6.00 -2.60] [ 6.00 2.60] [ 3.00 -2.60] [ 3.00 2.60] [ 1.50 -0.00] [ 1.50 5.20]   The code to identify the initial segments and the common segments is as follows. # -- bts ... is the fr_to = np.concatenate([np.concatenate((b[:-1], b[1:]), axis=1) for b in bts], axis=0) tst = np.concatenate((fr_to[:, 2:4], fr_to[:, :2]), axis=1) idx0, idx1 = np.nonzero((fr_to == tst[:, None]).all(-1)) idx_01 = np.asarray((idx0, idx1)).T common = fr_to[idx0]   Some explanation:   fr_to :  the flattened from-to coordinates as previously indicated for all the segments in all the shapes for b in bts : get all the coordinates from the first to the last join them together (concatenate) as a pair using :         np.concatenate((b[:-1], b[1:]), axis=1)          where  b[:-1]      from the first to the end but not including it                       b[1:]      from the second to the end including the end                       axis=1   concatenate along the rows np.concatenate( ... the above ..., axis=0) take all the segment pairs and join them by columns (axis=0) take the fr_to data and reverse the point pairs to form the 'tst' array, for example   fr_to[0]   =  array([ -3.00, -2.60, -1.50, 0.00])   # -- (x0, y0, x1, y1) tst[0]       =  array([ -1.50, 0.00, -3.00, -2.60])   # -- (x1, y1, x0, y0) Note the reversal of the pairs. Now line 5... compare all the fr_to coordinates to all the tst coordinates item by item and if the are 'all' equal, then note the index id value in each array. This returns the segments uncorrected for directionality common # -- step 1 array([[ -1.50, 0.00, 1.50, 0.00], [ 1.50, 0.00, 3.00, -2.60], [ 1.50, 0.00, 3.00, 2.60], [ 3.00, 2.60, 6.00, 2.60], [ 3.00, -2.60, 1.50, -0.00], [ 1.50, 5.20, 3.00, 2.60], [ 3.00, 2.60, 1.50, 0.00], [ 1.50, 0.00, -1.50, 0.00], [ 6.00, 2.60, 3.00, 2.60], [ 3.00, 2.60, 1.50, 5.20]])   The last step is to just get a set of coordinates common_pair_ids = np.asarray([i for i in idx_01 if i[0] < i[1]]) pair_0 = common_pair_ids[:, 0] common = fr_to[pair_0] Retrieve the information you want common_pair_ids array([[ 1, 16], [ 2, 11], [ 6, 15], [ 7, 22], [14, 23]]) common # step 2 array([[ -1.50, 0.00, 1.50, 0.00], [ 1.50, 0.00, 3.00, -2.60], [ 1.50, 0.00, 3.00, 2.60], [ 3.00, 2.60, 6.00, 2.60], [ 1.50, 5.20, 3.00, 2.60]]) # -- for example, look at the reversal of segments 1 and 16 fr_to[[1, 16]] array([[ -1.50, 0.00, 1.50, 0.00], [ 1.50, 0.00, -1.50, 0.00]])   Although seemingly complex, it really isn't.  The hardest thing to wrap your head around is the 'broadcasting' that is used in the 'np.nonzero' that is used to return the indices.    For the arcpy types, it is the equivalent of taking a segment, comparing it to all the other reversed segments looking for a match...  That is all for now, but add this to your arsenal when you need to find 'common' or 'shared' elements.  The principles go way beyond this simple example.   ... View more