.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples/05_preprocessing.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_examples_05_preprocessing.py: Preprocessing ============= Raw geospatial vector data is often not ready to use directly in mesh generation: * Polygon data often do not form a valid planar partition: polygons are overlapping, or neighboring polygons have small gaps between them. * Polygon boundaries or linestring segments intersect each other. * Points may be located on polygon boundaries or lines. Since floating point numbers are not exact, points seemingly located on a line are computationally just left or just right of the line and form an extremely thin triangle. * Points may be located extremely close together, thereby generating tiny triangles. Such problems either lead to a generated mesh with extremely small elements, or worse, they lead to a crash of the meshing program. Pandamesh provides a ``Preprocessor`` class to assist with cleaning up some common faults. This example will illustrate some common problems and how to resolve them. .. GENERATED FROM PYTHON SOURCE LINES 24-33 .. code-block:: Python import geopandas as gpd import matplotlib.pyplot as plt import numpy as np import shapely import shapely.geometry as sg import pandamesh as pm .. GENERATED FROM PYTHON SOURCE LINES 37-45 Polygons -------- When generating a mesh, we often have a general area which may be meshed coarsely and an area of interest, which should be meshed more finely. Generally, the fine inner zone is located within the coarse outer zone, but this requires a hole in the outer zone that exactly matches up with the exterior of the inner zone. .. GENERATED FROM PYTHON SOURCE LINES 45-70 .. code-block:: Python outer = sg.Polygon( [ [0.0, 0.0], [10.0, 0.0], [10.0, 10.0], [0.0, 10.0], ] ) inner = sg.Polygon( [ [5.0, 2.0], [8.0, 5.0], [5.0, 8.0], [2.0, 5.0], ] ) gdf = gpd.GeoDataFrame(geometry=[outer, inner]) gdf["cellsize"] = [2.0, 1.0] fig, (ax0, ax1) = plt.subplots(ncols=2, sharex=True, sharey=True) gdf.iloc[[0]].plot(ax=ax0) gdf.iloc[[1]].plot(ax=ax1) .. image-sg:: /examples/images/sphx_glr_05_preprocessing_001.png :alt: 05 preprocessing :srcset: /examples/images/sphx_glr_05_preprocessing_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 71-73 In this case, we have two conflicting specified cell sizes in the inner square. We can resolve this as follows: .. GENERATED FROM PYTHON SOURCE LINES 73-80 .. code-block:: Python resolved = ( pm.Preprocessor(geometry=gdf.geometry, values=gdf.cellsize) .unify_polygons() .to_geodataframe() ).rename(columns={"values": "cellsize"}) .. GENERATED FROM PYTHON SOURCE LINES 81-85 Note that the Preprocessor supports method chaining, allowing you to flexibly execute a set of operations. The resulting geodataframe's geometries are valid planar partition: .. GENERATED FROM PYTHON SOURCE LINES 85-90 .. code-block:: Python fig, (ax0, ax1) = plt.subplots(ncols=2, sharex=True, sharey=True) resolved.iloc[[0]].plot(ax=ax0) resolved.iloc[[1]].plot(ax=ax1) .. image-sg:: /examples/images/sphx_glr_05_preprocessing_002.png :alt: 05 preprocessing :srcset: /examples/images/sphx_glr_05_preprocessing_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 91-92 And we can use it directly to generate a mesh: .. GENERATED FROM PYTHON SOURCE LINES 92-96 .. code-block:: Python vertices, faces = pm.TriangleMesher(resolved).generate() pm.plot(vertices, faces) .. image-sg:: /examples/images/sphx_glr_05_preprocessing_003.png :alt: 05 preprocessing :srcset: /examples/images/sphx_glr_05_preprocessing_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 97-99 Alternatively, multiple polygons with the same cell size specification might be overlapping .. GENERATED FROM PYTHON SOURCE LINES 99-107 .. code-block:: Python inner0 = shapely.affinity.translate(inner, xoff=-1.0) inner1 = shapely.affinity.translate(inner, xoff=1.0) gdf = gpd.GeoDataFrame(geometry=[outer, inner0, inner1]) gdf["cellsize"] = [2.0, 1.0, 1.0] fig, ax = plt.subplots() gdf.plot(ax=ax, facecolor="none") .. image-sg:: /examples/images/sphx_glr_05_preprocessing_004.png :alt: 05 preprocessing :srcset: /examples/images/sphx_glr_05_preprocessing_004.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 108-109 These will also be resolved by ``.unify_polygons``. .. GENERATED FROM PYTHON SOURCE LINES 109-122 .. code-block:: Python resolved = ( pm.Preprocessor(geometry=gdf.geometry, values=gdf.cellsize) .unify_polygons() .to_geodataframe() ).rename(columns={"values": "cellsize"}) vertices, faces = pm.TriangleMesher(resolved).generate() fig, ax = plt.subplots() pm.plot(vertices, faces, ax=ax) resolved.plot(ax=ax, facecolor="none", edgecolor="red") .. image-sg:: /examples/images/sphx_glr_05_preprocessing_005.png :alt: 05 preprocessing :srcset: /examples/images/sphx_glr_05_preprocessing_005.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 123-125 Note, however, that the internal boundaries of the inner polygons are forced into the triangulation. We can rid of these by calling ``.merge_polygons``: .. GENERATED FROM PYTHON SOURCE LINES 125-139 .. code-block:: Python resolved = ( pm.Preprocessor(geometry=gdf.geometry, values=gdf.cellsize) .unify_polygons() .merge_polygons() .to_geodataframe() ).rename(columns={"values": "cellsize"}) vertices, faces = pm.TriangleMesher(resolved).generate() fig, ax = plt.subplots() pm.plot(vertices, faces, ax=ax) resolved.plot(ax=ax, facecolor="none", edgecolor="red") .. image-sg:: /examples/images/sphx_glr_05_preprocessing_006.png :alt: 05 preprocessing :srcset: /examples/images/sphx_glr_05_preprocessing_006.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 140-142 An alternative problem is when polygons are touching, but do not actually share vertices along the boundary. .. GENERATED FROM PYTHON SOURCE LINES 142-166 .. code-block:: Python first = sg.Polygon( [ [0.0, 0.0], [10.0, 0.0], [10.0, 10.0], [0.0, 10.0], ] ) second = sg.Polygon( [ [10.0, 2.0], [18.0, 2.0], [18.0, 8.0], [10.0, 8.0], ] ) gdf = gpd.GeoDataFrame(geometry=[first, second]) gdf["cellsize"] = [4.0, 2.0] vertices, faces = pm.GmshMesher(gdf, intersecting_edges="warn").generate(finalize=True) pm.plot(vertices, faces) .. image-sg:: /examples/images/sphx_glr_05_preprocessing_007.png :alt: 05 preprocessing :srcset: /examples/images/sphx_glr_05_preprocessing_007.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none /home/runner/work/pandamesh/pandamesh/pandamesh/common.py:407: UserWarning: 2 unresolved intersections between polygon boundary or line segments. Run pandamesh.find_edge_intersections(gdf.geometry) to identify the intersection locations. Intersections can be resolved using the pandamesh.Preprocessor. warnings.warn(msg) .. GENERATED FROM PYTHON SOURCE LINES 167-177 At x=10.0, the generated triangles are disconnected. This is caused by the the fact that the polygons do not share an edge: * The polygon on the left has an edge from (10.0, 0.0) to (10.0, 10.0) * The polygon on the right has an edge from (10.0, 2.0) to (10.0, 8.0) In fact, the vertices of the right polygon are intersecting the (edge) of the left polygon. We can identify these intersections with :func:`pandamesh.find_edge_intersections`: .. GENERATED FROM PYTHON SOURCE LINES 177-184 .. code-block:: Python intersections = pm.find_edge_intersections(gdf.geometry) fig, ax = plt.subplots() pm.plot(vertices, faces, ax=ax) intersections.plot(ax=ax) .. image-sg:: /examples/images/sphx_glr_05_preprocessing_008.png :alt: 05 preprocessing :srcset: /examples/images/sphx_glr_05_preprocessing_008.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 185-187 Calling ``.unify_polygons()`` ensures that the vertices of touching polygons are inserted, such that the polygons share an edge. .. GENERATED FROM PYTHON SOURCE LINES 187-201 .. code-block:: Python resolved = ( pm.Preprocessor(geometry=gdf.geometry, values=gdf.cellsize) .unify_polygons() .to_geodataframe() ).rename(columns={"values": "cellsize"}) vertices, faces = pm.TriangleMesher(resolved).generate() polygon0_coords = shapely.get_coordinates(resolved.geometry[0]) fig, ax = plt.subplots() pm.plot(vertices, faces, ax=ax) ax.scatter(*polygon0_coords.T) .. image-sg:: /examples/images/sphx_glr_05_preprocessing_009.png :alt: 05 preprocessing :srcset: /examples/images/sphx_glr_05_preprocessing_009.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 202-206 Lines ----- Lines may only be only partially present, or present in holes: .. GENERATED FROM PYTHON SOURCE LINES 206-240 .. code-block:: Python donut = sg.Polygon( [ [0.0, 0.0], [10.0, 0.0], [10.0, 10.0], [0.0, 10.0], ], holes=[ [ [2.0, 5.0], [5.0, 8.0], [8.0, 5.0], [5.0, 2.0], ] ], ) line0 = shapely.LineString( [ [-2.0, 0.0], [12.0, 10.0], ] ) line1 = shapely.LineString( [ [5.5, 9.0], [9.0, 5.5], ] ) gdf = gpd.GeoDataFrame(geometry=[donut, line0, line1]) gdf["cellsize"] = [2.0, 1.0, 1.0] gdf.plot(edgecolor="k") .. image-sg:: /examples/images/sphx_glr_05_preprocessing_010.png :alt: 05 preprocessing :srcset: /examples/images/sphx_glr_05_preprocessing_010.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 241-243 We can identify these problematic intersections again using :func:`pandamesh.find_edge_intersections`: .. GENERATED FROM PYTHON SOURCE LINES 243-249 .. code-block:: Python intersections = pm.find_edge_intersections(gdf.geometry) fig, ax = plt.subplots() gdf.plot(ax=ax, facecolor="none") intersections.plot(ax=ax) .. image-sg:: /examples/images/sphx_glr_05_preprocessing_011.png :alt: 05 preprocessing :srcset: /examples/images/sphx_glr_05_preprocessing_011.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 250-251 A first step is to remove line segments that do not fall in any polygon: .. GENERATED FROM PYTHON SOURCE LINES 251-259 .. code-block:: Python resolved = ( pm.Preprocessor(geometry=gdf.geometry, values=gdf.cellsize) .clip_lines() .to_geodataframe() ).rename(columns={"values": "cellsize"}) resolved.plot(edgecolor="k") .. image-sg:: /examples/images/sphx_glr_05_preprocessing_012.png :alt: 05 preprocessing :srcset: /examples/images/sphx_glr_05_preprocessing_012.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 260-264 However, this doesn't create suitable input for meshing. The ``GmshMesher`` appears to hang on this input, and Triangle generates a grid with very small triangles. Pandamesh errors on these intersections by default, but way may proceed: .. GENERATED FROM PYTHON SOURCE LINES 264-268 .. code-block:: Python vertices, faces = pm.TriangleMesher(resolved, intersecting_edges="warn").generate() pm.plot(vertices, faces) .. image-sg:: /examples/images/sphx_glr_05_preprocessing_013.png :alt: 05 preprocessing :srcset: /examples/images/sphx_glr_05_preprocessing_013.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none /home/runner/work/pandamesh/pandamesh/pandamesh/common.py:407: UserWarning: 5 unresolved intersections between polygon boundary or line segments. Run pandamesh.find_edge_intersections(gdf.geometry) to identify the intersection locations. Intersections can be resolved using the pandamesh.Preprocessor. warnings.warn(msg) .. GENERATED FROM PYTHON SOURCE LINES 269-283 A better approach here is to ensure all intersections are present in all linework: * First we clip. * Then we call ``unify_lines`` to ensure that the intersection of line0 and line1 at (7.625 6.875) is represented. * Next we call ``unify_polygons``. This ensures the intersections of the lines with the poygon exterior is represented as well. * The result of ``unify_polygons`` is that the line splits the polygon in two parts. These are merged back together with ``merge_polygons``. If we plot the vertices of the resolved polygon, we see that the intersection vertices have been inserted into the polygon boundaries, and that the tiny triangles around the line intersection have disappeared: .. GENERATED FROM PYTHON SOURCE LINES 283-300 .. code-block:: Python resolved = ( pm.Preprocessor(geometry=gdf.geometry, values=gdf.cellsize) .clip_lines() .unify_lines() .unify_polygons() .merge_polygons() .to_geodataframe() ).rename(columns={"values": "cellsize"}) vertices, faces = pm.GmshMesher(resolved).generate(finalize=True) polygon0_coords = shapely.get_coordinates(resolved.geometry[0]) fig, ax = plt.subplots() pm.plot(vertices, faces, ax=ax) ax.scatter(*polygon0_coords.T) .. image-sg:: /examples/images/sphx_glr_05_preprocessing_014.png :alt: 05 preprocessing :srcset: /examples/images/sphx_glr_05_preprocessing_014.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 301-304 In some cases, having line segments terminate exactly on polygon boundaries still causes trouble. We may also ensure that lines are some distance removed from any polygon boundary by providing a distance to ``clip_lines``: .. GENERATED FROM PYTHON SOURCE LINES 304-320 .. code-block:: Python resolved = ( pm.Preprocessor(geometry=gdf.geometry, values=gdf.cellsize) .unify_lines() .clip_lines(distance=0.5) .to_geodataframe() ).rename(columns={"values": "cellsize"}) vertices, faces = pm.GmshMesher(resolved).generate(finalize=True) polygon0_coords = shapely.get_coordinates(resolved.geometry[0]) fig, ax = plt.subplots() pm.plot(vertices, faces, ax=ax) resolved.plot(facecolor="none", edgecolor="red", ax=ax) ax.scatter(*polygon0_coords.T) .. image-sg:: /examples/images/sphx_glr_05_preprocessing_015.png :alt: 05 preprocessing :srcset: /examples/images/sphx_glr_05_preprocessing_015.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 321-323 Another pragmatic approach is to convert any line into interpolated points. Points cannot intersect each other, which sidesteps a large number of problems. .. GENERATED FROM PYTHON SOURCE LINES 323-337 .. code-block:: Python resolved = ( pm.Preprocessor(geometry=gdf.geometry, values=gdf.cellsize) .interpolate_lines_to_points(distance=0.25) .clip_points() .to_geodataframe() ).rename(columns={"values": "cellsize"}) vertices, faces = pm.GmshMesher(resolved).generate(finalize=True) fig, ax = plt.subplots() pm.plot(vertices, faces, ax=ax) resolved.plot(facecolor="none", edgecolor="red", ax=ax) .. image-sg:: /examples/images/sphx_glr_05_preprocessing_016.png :alt: 05 preprocessing :srcset: /examples/images/sphx_glr_05_preprocessing_016.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 338-345 Points ------ Note that the start and end points of the lines are still on, or very near the polygon edges. We can remove those points by providing a distance to ``clip_points``. .. GENERATED FROM PYTHON SOURCE LINES 345-359 .. code-block:: Python resolved = ( pm.Preprocessor(geometry=gdf.geometry, values=gdf.cellsize) .interpolate_lines_to_points(distance=0.25) .clip_points(distance=0.5) .to_geodataframe() ).rename(columns={"values": "cellsize"}) vertices, faces = pm.GmshMesher(resolved).generate(finalize=True) fig, ax = plt.subplots() pm.plot(vertices, faces, ax=ax) resolved.plot(facecolor="none", edgecolor="red", ax=ax) .. image-sg:: /examples/images/sphx_glr_05_preprocessing_017.png :alt: 05 preprocessing :srcset: /examples/images/sphx_glr_05_preprocessing_017.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 360-362 A problem with points is that they may be very close together, thereby generating very small triangles. Let's generate 200 random points to illustrate: .. GENERATED FROM PYTHON SOURCE LINES 362-377 .. code-block:: Python rng = np.random.default_rng() points = gpd.points_from_xy(*rng.random((2, 200)) * 10.0) gdf = gpd.GeoDataFrame(geometry=np.concatenate([[donut], points])) gdf["cellsize"] = 2.0 resolved = ( pm.Preprocessor(geometry=gdf.geometry, values=gdf.cellsize) .clip_points(distance=0.5) .to_geodataframe() ).rename(columns={"values": "cellsize"}) vertices, faces = pm.GmshMesher(resolved).generate(finalize=True) pm.plot(vertices, faces) .. image-sg:: /examples/images/sphx_glr_05_preprocessing_018.png :alt: 05 preprocessing :srcset: /examples/images/sphx_glr_05_preprocessing_018.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 378-380 We can solve this by snapping points together that are located some distance from each other: .. GENERATED FROM PYTHON SOURCE LINES 380-391 .. code-block:: Python resolved = ( pm.Preprocessor(geometry=gdf.geometry, values=gdf.cellsize) .clip_points(distance=0.5) .snap_points(distance=0.5) .to_geodataframe() ).rename(columns={"values": "cellsize"}) vertices, faces = pm.GmshMesher(resolved).generate(finalize=True) pm.plot(vertices, faces) .. image-sg:: /examples/images/sphx_glr_05_preprocessing_019.png :alt: 05 preprocessing :srcset: /examples/images/sphx_glr_05_preprocessing_019.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 392-405 Flexibility and composability ----------------------------- The Preprocessor class in Pandamesh is designed with flexibility and composability in mind through method chaining. By combining various preprocessing steps in any order, you can address a wide range of geometric issues. For instance, you might start by unifying polygons, then clip lines, interpolate them to points, and finally snap those points together. The steps required depend on the nature of geometrical input, and may require experimenting with various methods. The intermediate output can be checked and visualized at any moments, by calling ``to_geodataframe``. For example, to check the intermediate result after clipping but prior to snapping: .. GENERATED FROM PYTHON SOURCE LINES 405-414 .. code-block:: Python check = ( pm.Preprocessor(geometry=gdf.geometry, values=gdf.cellsize) .clip_points(distance=0.5) .to_geodataframe() ) check.plot(facecolor="none") .. image-sg:: /examples/images/sphx_glr_05_preprocessing_020.png :alt: 05 preprocessing :srcset: /examples/images/sphx_glr_05_preprocessing_020.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 415-420 This also makes it easy to apply the preprocessor in steps. Some steps may be relatively costly, such as unifying a large number of detailed polygons. The intermediate result can be stored as e.g. a GeoPackage. Then, in a separate processing step, the intermediate result can be read again, and other processing steps (such as filtering points) can be applied. .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 1.658 seconds) .. _sphx_glr_download_examples_05_preprocessing.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: 05_preprocessing.ipynb <05_preprocessing.ipynb>` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: 05_preprocessing.py <05_preprocessing.py>` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: 05_preprocessing.zip <05_preprocessing.zip>` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_