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.jo...@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 > chao...@ms.iswc.ac.cn > Mobile: +86 18709187546 > > ************************************************************************************ > > _______________________________________________ > gdal-dev mailing > listgdal-dev@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@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@t-kartor.sehttp://www.t-kartor.com > > > _______________________________________________ > gdal-dev mailing > listgdal-dev@lists.osgeo.orghttps://lists.osgeo.org/mailman/listinfo/gdal-dev > > _______________________________________________ > gdal-dev mailing list > gdal-dev@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 chao...@ms.iswc.ac.cn Mobile: +86 18709187546 ************************************************************************************
_______________________________________________ gdal-dev mailing list gdal-dev@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/gdal-dev
