<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
</head>
<body>
<p>Hi,</p>
<p>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.<br>
</p>
<p>Ari<br>
</p>
<div class="moz-cite-prefix">Chao YUE kirjoitti 15.4.2020 klo 6.30:<br>
</div>
<blockquote type="cite"
cite="mid:CAAN-aRFQsyoVTBceU0QhGcoPr-obmeVhUujAccYroaJwBBPUCw@mail.gmail.com">
<meta http-equiv="content-type" content="text/html; charset=UTF-8">
<div dir="ltr">
<div class="gmail_default" style="font-size:small">
<div>Dear all,</div>
<div><br>
</div>
<div>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.</div>
<div><br>
</div>
<div>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. </div>
<div>But I am wondering whether there is already some existing
solutions or other better ones. </div>
<div><br>
</div>
<div>Thanks a lot for the kind help for any hints on this !</div>
<div>Kind regards,</div>
<div>Chao</div>
<div><br>
</div>
<div>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.</div>
<div><br>
</div>
<div><img src="cid:part1.1F4D8229.94F736F6@gmail.com"
alt="InnerLoop.png" class="" width="562" height="415"><br>
</div>
</div>
<div><br>
</div>
-- <br>
<div dir="ltr" class="gmail_signature"
data-smartmail="gmail_signature">
<div dir="ltr">
<div>
<div dir="ltr">
<div>
<div dir="ltr">
<div dir="ltr">
<div dir="ltr">
<div><span style="font-size:12.8px">***********************************************************************************</span><br>
</div>
<div><font color="#000000">Chao YUE（岳超）</font></div>
<font color="#000000">西北农林科技大学水土保持研究所 研究员</font></div>
<div dir="ltr"><font color="#000000">黄土高原土壤侵蚀与旱地农业国家重点实验室</font></div>
<div dir="ltr"><span
style="font-family:Helvetica;font-size:12px"><font
color="#000000">State Key Laboratory of Soil
Erosion and Dryland Farming on the Loess
Plateau</font></span></div>
<div dir="ltr"><span
style="font-family:Helvetica;font-size:12px"><font
color="#000000">Institute of Soil and Water
Conservation</font></span></div>
<div dir="ltr"><font color="#000000"><span
style="font-family:Helvetica;font-size:12px">Northwest
A&F University, Yangling, Shaanxi
712100, P.R. China</span><br>
</font></div>
<div dir="ltr">
<div><font color="#000000"><a
href="mailto:chaoyue@ms.iswc.ac.cn"
style="font-family:Helvetica;font-size:12px"
target="_blank" moz-do-not-send="true">chaoyue@ms.iswc.ac.cn</a><br>
</font></div>
<div><font color="#000000">Mobile: +86
18709187546</font></div>
<div>************************************************************************************<br>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
<br>
<fieldset class="mimeAttachmentHeader"></fieldset>
<pre class="moz-quote-pre" wrap="">_______________________________________________
gdal-dev mailing list
<a class="moz-txt-link-abbreviated" href="mailto:gdal-dev@lists.osgeo.org">gdal-dev@lists.osgeo.org</a>
<a class="moz-txt-link-freetext" href="https://lists.osgeo.org/mailman/listinfo/gdal-dev">https://lists.osgeo.org/mailman/listinfo/gdal-dev</a></pre>
</blockquote>
<p><br>
</p>
<p>Code copyright Simosol Oy, licence MIT<br>
</p>
<p>def add_point(self, p=None):<br>
""" Add a point to the end of points (before root)<br>
The first level ear (hole or in-a-point-intersecting
polygon) is returned<br>
if there is such.<br>
p -- point to add, if None, close the ring</p>
<p> self is a Polygon, which has as attributes a dict of points,
where<br>
point = [index of previous point, index of next point, coords]<br>
no duplicates, no null indexes in points<br>
</p>
<p> """<br>
n = self.n_points()<br>
if p is None and n < 3:<br>
raise ValueError("Can't close a polygon with less than
three points.")<br>
if self.root is None:<br>
self.root = 0<br>
self.points[self.root] = [self.root, self.root, p]<br>
return<br>
i0 = self.points[self.root][0]<br>
if p is not None and p[0] == self.points[i0][2][0] and
p[1] == self.points[i0][2][1]:<br>
# same point, skip<br>
return None<br>
index = i0 + 1<br>
while index in self.points:<br>
index += 1<br>
ear = None<br>
close = False<br>
if n > 2:<br>
# test for self-intersection,<br>
# current - new against (0-1, )1-2, ...
prev(current)-current<br>
current_point = self.get_point(i0)<br>
j0 = self.root<br>
if p is None: # close<br>
p = self.get_point(self.root)<br>
j0 = self.index_of_next_point(j0)<br>
close = True<br>
while True:<br>
j1 = self.index_of_next_point(j0)<br>
if j1 == i0:<br>
break<br>
pj0 = self.get_point(j0)<br>
pj1 = self.get_point(j1)<br>
x = get_line_intersection(pj0, pj1, current_point,
p)<br>
<br>
if type(x) is tuple: # Intersection<br>
#print('intersection',pj0,pj1,current_point,p)<br>
# snip away and return eventually the closed
part<br>
ear = self.get_noose(j1, i0, first_point=x)<br>
# set new j1 and delete from j1+1 ... current<br>
self.set_point(j1, x)<br>
j2 = self.index_of_next_point(j1)<br>
self.delete_points(j2, i0)<br>
i0 = j0<br>
index = j1<br>
break<br>
<br>
elif abs(x) == 1: # Collinear<br>
if point_on_line(p, pj0, pj1): # new is on j0
- j1<br>
#print('collinear new on
j0-j1',pj0,pj1,current_point,p)<br>
# is j1 on cp - new?<br>
if point_in_box(pj1, current_point, p):<br>
k = j1 # ear is j1 ... cp<br>
else: # j0 is on cp - new<br>
k = j0 # ear is j0 ... cp<br>
ear = self.get_noose(k, i0)<br>
self.delete_points(j1, i0)<br>
i0 = j0<br>
index = j1<br>
break<br>
elif point_on_line(pj1, current_point, p): #
j1 is on cp - new<br>
#print('collinear j1 on
cp-new',pj0,pj1,current_point,p)<br>
# ear is j1 ... cp<br>
ear = self.get_noose(j1, i0)<br>
self.delete_points(j1, i0)<br>
i0 = j0<br>
index = j1<br>
break<br>
j0 = j1<br>
<br>
if not close:<br>
self.points[index] = [i0, self.root, p]<br>
self.points[i0][1] = index<br>
self.points[self.root][0] = index<br>
return ear<br>
</p>
</body>
</html>