[gdal-dev] Algorithm to clip the inner loop of a polygon
Chao YUE
chaoyuejoy at gmail.com
Wed Apr 15 06:08:16 PDT 2020
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> 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.
>
> [image: 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
> Mobile: +86 18709187546
>
> ************************************************************************************
>
> _______________________________________________
> gdal-dev mailing listgdal-dev at lists.osgeo.orghttps://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 listgdal-dev at lists.osgeo.orghttps://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.sehttp://www.t-kartor.com
>
>
> _______________________________________________
> gdal-dev mailing listgdal-dev at lists.osgeo.orghttps://lists.osgeo.org/mailman/listinfo/gdal-dev
>
> _______________________________________________
> gdal-dev mailing list
> 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
Mobile: +86 18709187546
************************************************************************************
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20200415/08fdea69/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/08fdea69/attachment-0001.png>
More information about the gdal-dev
mailing list