[gdal-dev] Algorithm to clip the inner loop of a polygon
Andreas Oxenstierna
ao at t-kartor.se
Tue Apr 14 23:03:40 PDT 2020
Hi
Have you tested PostGIS?
ST_ExteriorRing((ST_Dump(ST_LineMerge(ST_Union(ST_Force2D(<your
geometry>))))).geom)
should extract the points you want to keep by splitting at intersections
- but it may fail in very complex cases. Then I would start with a
ST_SnapToGrid with a quite high tolerance to lower the coordinate
resolution.
> 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
>
>
> _______________________________________________
> gdal-dev mailing list
> gdal-dev at lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/gdal-dev
--
Hälsningar
Andreas Oxenstierna
T-Kartor Geospatial AB
Olof Mohlins väg 12 Kristianstad
mobile: +46 733 206831
mailto: ao at t-kartor.se
http://www.t-kartor.com
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20200415/87965824/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/87965824/attachment-0001.png>
More information about the gdal-dev
mailing list