[gdal-dev] Algorithm to clip the inner loop of a polygon
Ari Jolma
ari.jolma at gmail.com
Wed Apr 15 00:04:15 PDT 2020
Andreas Oxenstierna kirjoitti 15.4.2020 klo 9.03:
> 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.
That PostGIS solution looks quite nice. Unfortunately I was not able to
use PostGIS in my case but were restricted to Shapely, which uses GEOS.
But then, PostGIS is mostly based on GEOS, so it should perhaps be
doable that way too. I'm not working or using that solution I gave right
now so won't test it.
Best,
Ari
>
>
>> 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
>
> _______________________________________________
> gdal-dev mailing list
> gdal-dev at lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/gdal-dev
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20200415/dabf2d26/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/dabf2d26/attachment-0001.png>
More information about the gdal-dev
mailing list