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
chao...@ms.iswc.ac.cn <mailto:chao...@ms.iswc.ac.cn>
Mobile: +86 18709187546
************************************************************************************
_______________________________________________
gdal-dev mailing list
gdal-dev@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
_______________________________________________
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev
