Common segments
Take a few polygons. Hexs will do.
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.
