.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples/04_triangle-geospatial.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_04_triangle-geospatial.py: Geospatial Triangle Example =========================== In this example we'll illustrate how to generate a mesh from a "real-world" geospatial vector dataset. .. GENERATED FROM PYTHON SOURCE LINES 9-16 .. code-block:: Python import geopandas as gpd import matplotlib.pyplot as plt import pandas as pd import shapely.geometry as sg import pandamesh as pm .. GENERATED FROM PYTHON SOURCE LINES 17-25 Overlap ------- We will get the data of a GeoJSON file describing the provinces of the Netherlands, and select only the name and geometry columns. We'll set the coordinate reference system to the Dutch national standard (EPSG:28992). Finally we set the name column to be used as index, so we can select provinces on name. .. GENERATED FROM PYTHON SOURCE LINES 25-31 .. code-block:: Python provinces = pm.data.provinces_nl().loc[:, ["name", "geometry"]] provinces = provinces.to_crs("epsg:28992") provinces.index = provinces["name"] gdf = provinces.copy() .. GENERATED FROM PYTHON SOURCE LINES 32-35 The mesh generation software cannot deal with overlap of polygons. To get rid of overlap, we can use the spatial functionality that geopandas provides. Let's check the polygons for overlap first. .. GENERATED FROM PYTHON SOURCE LINES 35-43 .. code-block:: Python overlap = gdf.overlay(gdf, how="intersection", keep_geom_type=True) overlap = overlap.loc[overlap["name_1"] != overlap["name_2"]] fig, ax = plt.subplots() gdf.plot(ax=ax) overlap.plot(edgecolor="red", ax=ax) .. image-sg:: /examples/images/sphx_glr_04_triangle-geospatial_001.png :alt: 04 triangle geospatial :srcset: /examples/images/sphx_glr_04_triangle-geospatial_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 44-50 Clean-up -------- There are many small overlaps visible at the province borders. We can generate a consistent polygon using a unary union. .. GENERATED FROM PYTHON SOURCE LINES 50-55 .. code-block:: Python union = sg.Polygon(gdf.unary_union) union_gdf = gpd.GeoDataFrame(geometry=[union]) union_gdf["cellsize"] = 10_000.0 .. rst-class:: sphx-glr-script-out .. code-block:: none /home/runner/work/pandamesh/pandamesh/examples/04_triangle-geospatial.py:51: DeprecationWarning: The 'unary_union' attribute is deprecated, use the 'union_all()' method instead. union = sg.Polygon(gdf.unary_union) .. GENERATED FROM PYTHON SOURCE LINES 56-72 Unfortunately, the province boundaries of this dataset no do align neatly and there are a number of small holes present. Some of these holes are not formed by inconsistencies, but by a small number of Belgian exclaves, `Baarle-Hertog`_. Simplify -------- We'll ignore the subtleties of international law for now and use geopandas to remove all blemishes by: * squeezing out the holes with ``.buffer`` * dissolving the buffered polygons into a single polygon with ``.dissolve`` * simplifying the dissolved polygon to avoid over-refinement with ``.simplify`` This creates a clean, and simpler, geometry. .. GENERATED FROM PYTHON SOURCE LINES 72-82 .. code-block:: Python simplified = gdf.copy() simplified.geometry = simplified.geometry.buffer(500.0) simplified["dissolve_column"] = 0 simplified = simplified.dissolve(by="dissolve_column") simplified.geometry = simplified.geometry.simplify(5_000.0) simplified["cellsize"] = 10_000.0 simplified.plot() .. image-sg:: /examples/images/sphx_glr_04_triangle-geospatial_002.png :alt: 04 triangle geospatial :srcset: /examples/images/sphx_glr_04_triangle-geospatial_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 83-85 Using this clean geometry, we can generate an unstructured grid with a fairly constant cell size. .. GENERATED FROM PYTHON SOURCE LINES 85-90 .. code-block:: Python mesher = pm.TriangleMesher(simplified) vertices, triangles = mesher.generate() pm.plot(vertices, triangles) .. image-sg:: /examples/images/sphx_glr_04_triangle-geospatial_003.png :alt: 04 triangle geospatial :srcset: /examples/images/sphx_glr_04_triangle-geospatial_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 91-107 For real work, buffering and simplifying will likely not suffice. See the preprocessing example to for an overview of common issues and how to apply pandamesh's Preprocessor class to resolve them. Local refinement ---------------- To set a zone of refinement, we can define an additional polygon. We need to ensure that no overlap occurs in the follwing steps: * select the geometry of a single province; * simplify its geometry to an appropriate level of detail; * specify a smaller cell size; * remove this province from the enveloping polygon; * collect the two polygons in a single geodataframe. .. GENERATED FROM PYTHON SOURCE LINES 107-117 .. code-block:: Python utrecht = gdf.loc[["Utrecht"]] utrecht.geometry = utrecht.geometry.simplify(2_500.0) utrecht["cellsize"] = 5000.0 envelope = simplified.overlay(utrecht, how="difference") refined = pd.concat([envelope, utrecht]) refined.index = [0, 1] refined.plot(column="name") .. image-sg:: /examples/images/sphx_glr_04_triangle-geospatial_004.png :alt: 04 triangle geospatial :srcset: /examples/images/sphx_glr_04_triangle-geospatial_004.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 118-119 This results in a mesh with a smaller cell size in the province of Utrecht. .. GENERATED FROM PYTHON SOURCE LINES 119-124 .. code-block:: Python mesher = pm.TriangleMesher(refined) vertices, triangles = mesher.generate() pm.plot(vertices, triangles) .. image-sg:: /examples/images/sphx_glr_04_triangle-geospatial_005.png :alt: 04 triangle geospatial :srcset: /examples/images/sphx_glr_04_triangle-geospatial_005.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 125-137 Conclusion ---------- This example provides a taste of how to convert a geospatial vector dataset into an unstructured grid with a locally refined part. Real-world data generally come with their own idiosyncrasies and inconsistencies. Depending on the nature of the necessary fixes, they can be solved with geopandas functionality, but sometimes manual editing is required. Fortunately, geopandas provides easy input and output for many file formats, which can be opened by e.g. QGIS. .. _Baarle-Hertog: https://en.wikipedia.org/wiki/Baarle-Hertog .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.528 seconds) .. _sphx_glr_download_examples_04_triangle-geospatial.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: 04_triangle-geospatial.ipynb <04_triangle-geospatial.ipynb>` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: 04_triangle-geospatial.py <04_triangle-geospatial.py>` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: 04_triangle-geospatial.zip <04_triangle-geospatial.zip>` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_