[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