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 <mailto: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 <mailto:gdal-dev@lists.osgeo.org>
https://lists.osgeo.org/mailman/listinfo/gdal-dev
