Latest Contributions by DanPatterson

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

Take two polygons, an area of interest rectangle/square and an overlay polygon with the same extent as the first polygon, but rotated within it's bounds.  Now this could be the basis of clipping or more generally boolean operations on polygons.  From this data we can determine all of the associated overlay cases, like symmetrical difference, intersection, union, etcetera. Points 0, 1, 2, 3 belong to the overlay polygon.  A 4th point is the closing duplicate of point 0. Points 5, 6, 7, 8 belong to the polygon on the bottom.  Again, a 9th point is a duplicate of the first.   As arrays, you can slice the first and last values using [0, -1] where -1 means count backwards to the last.       aoi[[0, -1]], aoi0[[0, -1]] # slice the first and last points (array([[ 0.00, 0.00], [ 0.00, 0.00]]), array([[ 0.00, 4.00], [ 0.00, 4.00]])) When these shapes are intersected you end up with a combined geometry with each part consisting of from-to points as shown below (extra information on clockwise/counterclockwise and part and/or bit IDs are included) prn_tbl(geom2) ... OID_ Fr_pnt To_pnt CW/CCW Part_ID Bit_ID ---------------------------------------------------------------- 000 0 0 4 1 1 0 001 1 4 8 1 1 0 002 2 8 12 1 1 0 003 3 12 16 1 1 0 004 4 16 21 1 1 0 The points are arranged in clockwise order and are numbered as shown in the first figure. pnt shape part X Y -------------------------------- 000 0 0.00 4.00 001 0 0.00 10.00 002 0 4.00 10.00 003 0 ___ 0.00 4.00 004 1 -o 4.00 10.00 005 1 10.00 10.00 006 1 10.00 6.00 007 1 ___ 4.00 10.00 008 2 -o 10.00 6.00 009 2 10.00 0.00 010 2 6.00 0.00 011 2 ___ 10.00 6.00 012 3 -o 6.00 0.00 013 3 0.00 0.00 014 3 0.00 4.00 015 3 ___ 6.00 0.00 016 4 -o 0.00 4.00 017 4 4.00 10.00 018 4 10.00 6.00 019 4 6.00 0.00 020 4 0.00 4.00 The points, and hence their segments, can be classed as either on, out, or in.  Duplicate ids are kept in this example so you can follow the construction of the sub-polygons. on    :    [0, 1, 2, 3] out   :   [5, 6, 7, 8] in_    :   [] So an `out` point would be used to identify polygon parts that are exterior to the common overlay/clipped area. segs_out :   [[0, 5, 1, 1, 0],    [1, 6, 2, 2, 1],    [2, 7, 3, 3, 2],    [3, 8, 0, 0, 3]] segs_in :   [0, 1, 1, 2, 2, 3, 3, 0] In array format # sequence ID values [array([0, 1]), array([0, 5, 1]), array([1, 2]), array([1, 6, 2]), array([2, 3]), array([2, 7, 3]), array([3, 0]), array([3, 8, 0])] # lexicographic sort. Increasing x and decreasing y [array([0, 1]), array([0, 5, 1]), array([0, 3]), array([0, 8, 3]), array([1, 2]), array([1, 6, 2]), array([3, 2]), array([3, 7, 2])] In visual form This forms the basis of "sweep-line" or "plane-sweep" algorithms. In the figure to the right, there are 4 sweep lines which are identified where the two polygons share a common point.  The common points are 0, 1, 2, 3). Note that the original order differs from the lexicographic order since point id 3 is encounted before number 2.         Lexicographic sorting uses "keys" to specify which values to use in considering the sort.  Consider these examples of how the key value (in this case, x or y column) and their value (positive or negative) affects point sort order. kys0 = np.lexsort((aoi[:, -1], aoi[:, 0])) # increasing x, decreasing y kys1 = np.lexsort((-aoi[:, -1], aoi[:, 0])) # increasing x, increasing y kys2 = np.lexsort((aoi[:, 0], aoi[:, -1])) # increasing y, increasing x kys3 = np.lexsort((-aoi[:, 0], aoi[:, -1])) # increasing y, increasing y kys0, kys1, kys2, kys3 (array([0, 4, 1, 3, 2]), # -- note 0 and 4 are duplicate array([1, 0, 4, 2, 3]), array([0, 4, 3, 1, 2]), array([3, 0, 4, 2, 1])) aoi[kys0] array([[ 0.00, 0.00], [ 0.00, 0.00], [ 0.00, 10.00], [ 10.00, 0.00], [ 10.00, 10.00]]) aoi[kys1] array([[ 0.00, 10.00], [ 0.00, 0.00], [ 0.00, 0.00], [ 10.00, 10.00], [ 10.00, 0.00]]) aoi[kys2] array([[ 0.00, 0.00], [ 0.00, 0.00], [ 10.00, 0.00], [ 0.00, 10.00], [ 10.00, 10.00]]) aoi[kys3] array([[ 10.00, 0.00], [ 0.00, 0.00], [ 0.00, 0.00], [ 10.00, 10.00], [ 0.00, 10.00]]) Normally in a lexicograph sort, the last key is the primary sort key and preceding keys are used as the secondary/tertiary/etc sort keys. So don't just use "plain" sorting... consider "plane" sort or its fancy term "lexicographic" sorting. Things can get rather complex when multiple points share common x and/or y values... have a look.                 Or more simply.  Four lines and sorted with various orders. a # -- common start point (0, 0), increasing end points array([[ 0, 0, 10, 0], [ 0, 0, 10, 3], [ 0, 0, 10, 6], [ 0, 0, 10, 10]]) # sort by decreasing y value of the end point wh1_ = np.lexsort((a[:, 2], -a[:, 3], a[:, 1], a[:, 0])) # the keys a[wh1_] array([[ 0, 0, 10, 10], [ 0, 0, 10, 6], [ 0, 0, 10, 3], [ 0, 0, 10, 0]]) ... View more