[gdal-dev] Algorithm to clip the inner loop of a polygon
Andreas Oxenstierna
ao at t-kartor.se
Wed Apr 15 06:31:45 PDT 2020
Hi
I see that I forgot to create polygons from the splitted lines
Sp the syntax should be something like:
ST_ExteriorRing((ST_Dump(ST_Polygonize(ST_LineMerge(ST_Union(ST_Force2D(<your
geometry>)))))).geom)
> Hi Ari and Andreas,
>
> Thanks for the information. I will check it out. Today I advanced a
> little by almost implementing what's been proposed in my message there.
>
> Regards,
> Chao
>
> On Wed, Apr 15, 2020 at 3:04 PM Ari Jolma <ari.jolma at gmail.com
> <mailto:ari.jolma at gmail.com>> wrote:
>
> 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 <mailto: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 <mailto: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 <mailto:ao at t-kartor.se>
>> http://www.t-kartor.com
>>
>> _______________________________________________
>> gdal-dev mailing list
>> gdal-dev at lists.osgeo.org <mailto:gdal-dev at lists.osgeo.org>
>> https://lists.osgeo.org/mailman/listinfo/gdal-dev
> _______________________________________________
> gdal-dev mailing list
> gdal-dev at lists.osgeo.org <mailto:gdal-dev at lists.osgeo.org>
> https://lists.osgeo.org/mailman/listinfo/gdal-dev
>
>
>
> --
> ***********************************************************************************
> 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
> ************************************************************************************
--
Best Regards
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/f823edf6/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/f823edf6/attachment-0001.png>
More information about the gdal-dev
mailing list