Deltares / GEOLib

GEOLib: Python wrappers around the input and output files of the Deltares D-Serie models
https://deltares.github.io/GEOLib/
MIT License
22 stars 17 forks source link

DStability - Shapely accuracy leads to invalid new layers #172

Open breinbaas opened 10 months ago

breinbaas commented 10 months ago

In the current version of geolib adding new layers to a DStability geometry can lead to errors in the geometry. This seems to be related to the accuracy of the coordinates. Here is an example of a layer that has been added programmatically;

01

As you can see the rightmost point of the new layer is not connected to the geometry which leads to an error if you want to calculate this using DStability.

breinbaas commented 10 months ago

The solution is easy and implemented in https://github.com/Deltares/GEOLib/pull/171

Since it was a very small code change I have added it to the other PR but if you want to keep this seperate the code to be changed is;

def connect_layers(self, layer1: PersistableLayer, layer2: PersistableLayer):
        """Connects two polygons by adding a the missing points on the polygon edges. Returns the two new polygons."""
        linestring1 = self.to_shapely_linestring(layer1.Points)
        linestring2 = self.to_shapely_linestring(layer2.Points)

        # Create a union of the two polygons and polygonize it creating two connected polygons
        union = linestring1.union(linestring2, grid_size=1e-3) # ADDED GRIDSIZE PARAMETER
        result = [geom for geom in polygonize(union)]

        # If the result has two polygons, we return them, otherwise we return the original polygons
        if len(result) == 2:
            return result[0].exterior, result[1].exterior
        else:
            return linestring1, linestring2

Now everything looks fine;

02