[gdal-dev] Algorithm to clip the inner loop of a polygon
Ari Jolma
ari.jolma at gmail.com
Tue Apr 14 21:22:11 PDT 2020
Hi,
I wrote a method to a Polygon class (in Python) that solves that problem
because the ones in GEOS via Shapely did not seem good enough. My method
was based on no published algorithm but just my own reasoning, it is
probably not optimal speed-wise. I did not perform formal testing on it
but it worked on my data, which was quite large. The basic idea was to
build the polygon by adding points and when adding test for self
intersections / collinearity against existing segments. The method
returns the loops while points are added. The code is not in a public
repo but I'll add the method below. It needs some basic methods and
function which should be simple to add.
Ari
Chao YUE kirjoitti 15.4.2020 klo 6.30:
> Dear all,
>
> Does anyone have some experience or is aware of some algorithm that
> can find and clip the inner loop formed in a polygon ? I attach one
> example here. In this case I would only keep the outer points and drop
> the ones that make an inner loop. I am developing some algorithm to
> simulate wildland fire propagation. The algorithm is based on Richards
> 1990.
>
> In the paper he described an algorithm based on two steps: (1) find
> the points where a concave curvature is made. (2) search for both
> sides of this point to see where any two line segments cross over each
> other.
> But I am wondering whether there is already some existing solutions or
> other better ones.
>
> Thanks a lot for the kind help for any hints on this !
> Kind regards,
> Chao
>
> Gwynfor Richards, 1990. An elliptical growth model of forest fire
> fronts and its numerical solution. International journal for Numerical
> Methods in Engineering, Vol. 30, 1163-1179.
>
> InnerLoop.png
>
> --
> ***********************************************************************************
> Chao YUE（岳超）
> 西北农林科技大学水土保持研究所 研究员
> 黄土高原土壤侵蚀与旱地农业国家重点实验室
> State Key Laboratory of Soil Erosion and Dryland Farming on the Loess
> Plateau
> Institute of Soil and Water Conservation
> Northwest A&F University, Yangling, Shaanxi 712100, P.R. China
> chaoyue at ms.iswc.ac.cn <mailto:chaoyue at ms.iswc.ac.cn>
> Mobile: +86 18709187546
> ************************************************************************************
>
> _______________________________________________
> gdal-dev mailing list
> gdal-dev at lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/gdal-dev
Code copyright Simosol Oy, licence MIT
def add_point(self, p=None):
""" Add a point to the end of points (before root)
The first level ear (hole or in-a-point-intersecting polygon)
is returned
if there is such.
p -- point to add, if None, close the ring
self is a Polygon, which has as attributes a dict of points, where
point = [index of previous point, index of next point, coords]
no duplicates, no null indexes in points
"""
n = self.n_points()
if p is None and n < 3:
raise ValueError("Can't close a polygon with less than
three points.")
if self.root is None:
self.root = 0
self.points[self.root] = [self.root, self.root, p]
return
i0 = self.points[self.root][0]
if p is not None and p[0] == self.points[i0][2][0] and p[1] ==
self.points[i0][2][1]:
# same point, skip
return None
index = i0 + 1
while index in self.points:
index += 1
ear = None
close = False
if n > 2:
# test for self-intersection,
# current - new against (0-1, )1-2, ... prev(current)-current
current_point = self.get_point(i0)
j0 = self.root
if p is None: # close
p = self.get_point(self.root)
j0 = self.index_of_next_point(j0)
close = True
while True:
j1 = self.index_of_next_point(j0)
if j1 == i0:
break
pj0 = self.get_point(j0)
pj1 = self.get_point(j1)
x = get_line_intersection(pj0, pj1, current_point, p)
if type(x) is tuple: # Intersection
#print('intersection',pj0,pj1,current_point,p)
# snip away and return eventually the closed part
ear = self.get_noose(j1, i0, first_point=x)
# set new j1 and delete from j1+1 ... current
self.set_point(j1, x)
j2 = self.index_of_next_point(j1)
self.delete_points(j2, i0)
i0 = j0
index = j1
break
elif abs(x) == 1: # Collinear
if point_on_line(p, pj0, pj1): # new is on j0 - j1
#print('collinear new on
j0-j1',pj0,pj1,current_point,p)
# is j1 on cp - new?
if point_in_box(pj1, current_point, p):
k = j1 # ear is j1 ... cp
else: # j0 is on cp - new
k = j0 # ear is j0 ... cp
ear = self.get_noose(k, i0)
self.delete_points(j1, i0)
i0 = j0
index = j1
break
elif point_on_line(pj1, current_point, p): # j1 is
on cp - new
#print('collinear j1 on
cp-new',pj0,pj1,current_point,p)
# ear is j1 ... cp
ear = self.get_noose(j1, i0)
self.delete_points(j1, i0)
i0 = j0
index = j1
break
j0 = j1
if not close:
self.points[index] = [i0, self.root, p]
self.points[i0][1] = index
self.points[self.root][0] = index
return ear
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20200415/6aab9715/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: InnerLoop.png
Type: image/png
Size: 30127 bytes
Desc: not available
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20200415/6aab9715/attachment-0001.png>
More information about the gdal-dev
mailing list