Latest Contributions by DanPatterson

Take a square representing a polygon or a closed-loop polyline.  Let us call it aoi, perhaps out area of interest. aoi array([[ 0.00, 0.00], [ 0.00, 10.00], [ 10.00, 10.00], [ 10.00, 0.00], [ 0.00, 0.00]]) aoi.shape (5, 2) Notes: The first and last points are the same, hence it is a closed geometry The shape of the poly* feature is 5 rows (0 through 4 inclusive) and 2 columns ( 0 and 1),  since python numbering is zero-based Conceptually we know from out maths, that we just need to take the sequential difference between the points to get the differences in the x's and y's.  Whip in a quick euclidean calculation courtesy of the right-angle triangle and other cringeful recollections of things numeric. So it comes as no surprise that anything that involves rows and columns would be painfully slow to do one calculation at a time.  Enter einsum notation and its NumPy  representation.  (also useful for you Tensor types out there) diff = aoi[1:] - aoi[:-1] np.sqrt(np.einsum('ij,ij->i', diff, diff)) array([ 10.00, 10.00, 10.00, 10.00]) Which means multiply the two columns, and summarize, just like diff array([[ 0.00, 10.00], [ 10.00, 0.00], [ 0.00, -10.00], [-10.00, 0.00]]) diff * diff array([[ 0.00, 100.00], [ 100.00, 0.00], [ 0.00, 100.00], [ 100.00, 0.00]]) np.sum(diff*diff, axis= 1) # -- axis 1 means by row, since axis 0 is columns array([ 100.00, 100.00, 100.00, 100.00]) So you want the hypotenuse of your square?  Just use einsum summation on the other axis. np.sqrt(np.einsum('ij,ij->j', diff, diff)) # -- note summation is by 'j' this time array([ 14.14, 14.14]) The einsum notation string, just needs to be consistent when referencing the axes to perform actions on np.sqrt(np.einsum('ab,ab->a', diff, diff))  is the same as np.sqrt(np.einsum('ij,ij->j', diff, diff)) I use einsum in a lot of my Geo array calculations.  For example to determine the area of a polygon "bit" (either an inner or outer ring or parts of multipart polygons). def _bit_area_(a): """Mini e_area, used by `areas` and `centroids`. Negative areas are holes. This is intentionally reversed from the `shoelace` formula. """ a = _get_base_(a) x0, y1 = (a.T)[:, 1:] # cross set up as follows x1, y0 = (a.T)[:, :-1] e0 = np.einsum('...i,...i->...i', x0, y0) # 2024-03-28 modified e1 = np.einsum('...i,...i->...i', x1, y1) return np.sum((e0 - e1) * 0.5) This essentially sets up the 2D crossproduct for a modified 'shoelace' formula. How about determining the distance from one point to a bunch of other points? I will use 'aoi' as my other points and define a point 'p'. def _e_2d_(a, p): """Return distance from a point to an array of points""" diff = a - p[None, :] return np.sqrt(np.einsum('ij,ij->i', diff, diff)) p = np.array([3., 4.]) _e_2d_(aoi, p) array([ 5.00, 6.71, 9.22, 8.06, 5.00]) And to end this off, how about a translation and rotation of a polygon?  A little use of einsum to simplify the matrix calculations. def _r_(a, cent, angle, clockwise): """Rotate by `angle` in degrees, about the center.""" angle = np.radians(angle) if clockwise: angle = -angle c, s = np.cos(angle), np.sin(angle) R = np.array(((c, -s), (s, c))) return np.einsum('ij,jk->ik', a - cent, R) + cent Using the function above with these inputs _r_(aoi, cent=[1., 1.], angle=45., clockwise=True) array([[ 1.00, -0.41], [ -6.07, 6.66], [ 1.00, 13.73], [ 8.07, 6.66], [ 1.00, -0.41]]) yields the above coordinates which looks like this einsum_rotate_poly.png Enough for now. ... View more

It began as a simple coding exercise, but it soon morphed into one of those "things". Offset Buffer      I began with the seemingly simple exercise of the offset buffer, which is created by 'offsetting' the segments of a geometry object, polygons in this case, by a finite distance. I will spare you the sheer excitement when I managed to implement an offset buffer for an axis aligned square.  Here is a concave polygon with an offset buffer. offset_simple.png   For the most part the buffer is equally spaced except for those outward pointing corners. On to the rounded corners.  That involved arcs and angles and intersections and the math stuff that I had forgotten or long not used.  The effort was breathtaking.  I left the arc nodes to show I worked really hard. rounded_simple.png On to the next shapes! What could go wrong. Here is another example of a more convoluted concave polygon as an offset and conventional buffer. d1_with_np_wn_algo.png d1_with_np_wn_rounded.png                 I wonder what a really cool geometry would look like with the new buffering code? maple0.png The journey continue smoothly there.  Producing rounded corners when the pointy bits are close together proved to a challenge, but that involves another blog and the code links.  But in the interim, something to look forward to and share with your neighbours. maple_leaf.png   ... View more