Fragment model

README.rst

@@ -150,6 +150,30 @@ Subscription usage::
     async for item in subscription.enumerate():
           # do something with item
 
+Local development
+-----------------
+
+In order to set up a virtual environment, perform the following steps:
+
+Clone the repo and fetch the LFS objects::
+
+    git lfs fetch origin refs/remotes/origin/master
+    git lfs checkout
+
+Install platform dependencies::
+
+    sudo apt-get update && sudo apt-get install --yes --no-install-recommends libgdal-dev
+
+Create and activate a virtual environment::
+
+    python -m venv ./venv
+    source ./venv/bin/activate
+
+Install the dependencies. For your distribution, check the dependency matrix in .github/workflows/test.yml. For example, for Python 3.10::
+
+    pip install --disable-pip-version-check --upgrade pip setuptools wheel
+    pip install -e .[geo,results] pygdal==$(gdal-config --version).* ipython pytest flake8 sphinx==1.8.5 docutils==0.17.* sphinx_rtd_theme>=0.4.3 numpy==1.23.* h5py==3.7.* shapely==1.8.* pyproj==3.4.* geojson==2.5.* mercantile==1.2.1 cftime==1.6.2
+
 
 Credits
 -------

tests/conftest.py

@@ -35,6 +35,17 @@ def ga(request):
         yield ga
 
 
+@pytest.fixture()
+def ga_fragment_path():
+    return os.path.join(test_file_dir, "fragments", "gridadmin_fragments.h5")
+
+
+@pytest.fixture()
+def ga_fragments(ga_fragment_path):
+    with GridH5Admin(ga_fragment_path) as ga:
+        yield ga
+
+
 @pytest.fixture
 def gr_bergen_with_boundaries():
     """

tests/test_exporter.py

@@ -1,10 +1,13 @@
 import json
 import os
+from pathlib import Path
 
 import pytest
 from osgeo import ogr
 
+from threedigrid.admin.exporters.geopackage import GeopackageExporter
 from threedigrid.admin.exporters.geopackage.exporter import GpkgExporter
+from threedigrid.admin.fragments.exporters import FragmentsOgrExporter
 from threedigrid.admin.lines.exporters import LinesOgrExporter
 
 test_file_dir = os.path.join(os.getcwd(), "tests/test_files")
@@ -36,6 +39,68 @@ def test_nodes_gpgk_export(ga, tmp_path):
     assert layer.GetFeatureCount() == nodes_2d_open_water.id.size
 
 
+def test_fragments_gpgk_export(ga_fragments, tmp_path):
+    path = str(tmp_path / ("exporter_test_fragments.gpkg"))
+    exporter = FragmentsOgrExporter(ga_fragments.fragments)
+    exporter.save(path, "fragments", "4326")
+    assert os.path.exists(path)
+    s = ogr.Open(path)
+    layer = s.GetLayer()
+    assert layer.GetFeatureCount() == (
+        ga_fragments.fragments.id.size - 1
+    )  # minus dummy
+
+
+def test_fragments_to_gpgk_export(ga_fragments, tmp_path):
+    path = str(tmp_path / "only_fragment_export.gpkg")
+    ga_fragments.fragments.to_gpkg(path, "fragments")
+    s = ogr.Open(path)
+    layer = s.GetLayer("fragments")
+    assert layer.GetFeatureCount() == (
+        ga_fragments.fragments.id.size - 1
+    )  # minus dummy
+
+
+def test_fragments_to_geojson_export(ga_fragments, tmp_path):
+    path = str(tmp_path / "only_fragment_export.json")
+    ga_fragments.fragments.to_geojson(path, use_ogr=True)
+    s = ogr.Open(path)
+    layer = s.GetLayer(Path(path).name)
+    assert layer.GetFeatureCount() == (
+        ga_fragments.fragments.id.size - 1
+    )  # minus dummy
+    with open(path) as file:
+        data = json.load(file)
+        assert len(data["features"][0]["properties"]) == 2
+        assert "node_id" in data["features"][0]["properties"]
+
+    path = str(tmp_path / "only_fragment_export_2.json")
+    ga_fragments.fragments.to_geojson(path, use_ogr=False)
+    s = ogr.Open(path)
+    layer = s.GetLayer(Path(path).stem)
+    # We export more properties and use different layer name
+    assert layer.GetFeatureCount() == (
+        ga_fragments.fragments.id.size - 1
+    )  # minus dummy
+    with open(path) as file:
+        data = json.load(file)
+        assert len(data["features"][0]["properties"]) == 3
+        assert "node_id" in data["features"][0]["properties"]
+        assert "model_type" in data["features"][0]["properties"]
+
+
+def test_fragments_gpkg_conversion(ga_fragments, ga_fragment_path, tmp_path):
+    exporter = GeopackageExporter(
+        ga_fragment_path, str(tmp_path / "fragment_export.gpkg")
+    )
+    exporter.export()
+    s = ogr.Open(str(tmp_path / "fragment_export.gpkg"))
+    layer = s.GetLayer("fragment")
+    assert layer.GetFeatureCount() == (
+        ga_fragments.fragments.id.size - 1
+    )  # minus dummy
+
+
 def test_meta_data_gpgk_export(ga, tmp_path):
     path = str(tmp_path / ("exporter_meta.gpkg"))
     exporter = GpkgExporter(ga)
@@ -138,6 +203,7 @@ def test_export_null(ga, tmp_path):
         ("grid", "grid_2D_open_water"),
         ("grid", "grid_2D_groundwater"),
         ("flowlines", "flowlines"),
+        ("fragments", "fragments"),
     ],
 )
 def test_export_method(ga_export, export_method, expected_filename):

tests/test_gridadmin.py

@@ -572,3 +572,29 @@ def test_lines_filter_id_in(self):
         )
         expected = self.parser.lines.data["line_coords"][:, trues]
         np.testing.assert_equal(filtered, expected)
+
+
+class FragmentFilterTests(unittest.TestCase):
+    def setUp(self):
+        grid_admin_h5_file = os.path.join(
+            test_file_dir, "fragments", "gridadmin_fragments.h5"
+        )
+        self.parser = GridH5Admin(grid_admin_h5_file)
+
+    def test_fragments_filter_id_eq(self):
+
+        filtered = self.parser.fragments.filter(id=3).data["node_id"]
+        (trues,) = np.where(self.parser.fragments.data["id"] == 3)
+        expected = self.parser.fragments.data["node_id"][trues]
+        np.testing.assert_equal(filtered, expected)
+
+    def test_fragments_filter_id_in(self):
+        filtered = self.parser.fragments.filter(id__in=list(range(3, 7))).data[
+            "node_id"
+        ]
+        (trues,) = np.where(
+            (self.parser.fragments.data["id"] >= 3)
+            & (self.parser.fragments.data["id"] < 7)
+        )
+        expected = self.parser.fragments.data["node_id"][trues]
+        np.testing.assert_equal(filtered, expected)

threedigrid/admin/exporter_constants.py

@@ -26,6 +26,7 @@
         "line": ogr.wkbLineString,
         "multiline": ogr.wkbLineString,
         "bbox": ogr.wkbPolygon,
+        "polygon": ogr.wkbPolygon,
     }
 
 SHP_DRIVER_NAME = "ESRI Shapefile"
@@ -178,6 +179,8 @@
     "z_coordinate",
 ]
 
+FRAGMENTS_EXPORT_FIELDS = ["id", "node_id"]
+
 DEFAULT_EXPORT_FIELDS = {
     "Lines": "ALL",
     "Pipes": PIPES_EXPORT_FIELDS,
@@ -193,4 +196,5 @@
     "Breaches": BREACHES_EXPORT_FIELDS,
     "Levees": LEVEES_EXPORT_FIELDS,
     "Pumps": PUMPS_EXPORT_FIELDS,
+    "Fragments": FRAGMENTS_EXPORT_FIELDS,
 }

threedigrid/admin/exporters/geopackage/default.py

@@ -61,6 +61,7 @@ def progress_func(count, total):
             lines = ga.lines.filter(id__gte=1)
             nodes = ga.nodes.filter(id__gte=1)
             obstacles = ga.lines.filter(kcu=101)
+            fragments = ga.fragments.filter(id__gte=1)
 
             # Linestring geometry for pumps
             pumps_linestring = pumps.filter(node1_id__ne=-9999, node2_id__ne=-9999)
@@ -76,6 +77,7 @@ def progress_func(count, total):
                     + nodes.count
                     + obstacles_count
                     + pumps_linestring.count
+                    + fragments.count
                 )
                 self.start = 0
 
@@ -104,6 +106,7 @@ def internal_progress_func(item_count, item_total):
             )
             nodes.to_gpkg(self.gpkg_filename, progress_func=internal_func)
             cells.to_gpkg(self.gpkg_filename, progress_func=internal_func)
+            fragments.to_gpkg(self.gpkg_filename, progress_func=internal_func)
 
             if obstacles_count > 0:
                 # override the value for 'cross_pixel_coords'

threedigrid/admin/exporters/geopackage/exporter.py

@@ -13,6 +13,7 @@
 
 from threedigrid.admin import exporter_constants as const
 from threedigrid.geo_utils import get_spatial_reference
+from threedigrid.numpy_utils import reshape_flat_array
 from threedigrid.orm.base.exporters import BaseOgrExporter
 
 logger = logging.getLogger(__name__)
@@ -39,6 +40,12 @@ def get_geometry(field_name, field_type, data, index, **kwargs):
 
         # Create polygon from ring
         geom.AddGeometry(ring)
+    elif field_type == "polygon":
+        ring = ogr.Geometry(ogr.wkbLinearRing)
+        polygon_points = reshape_flat_array(data[field_name][index]).T
+        for x in polygon_points:
+            ring.AddPoint(x[0], x[1])
+        geom.AddGeometry(ring)
     elif field_type == "line":
         view = data[field_name].astype("float64")
         # print(data[field_name][:, index], data[field_name][:, index].dtype)

threedigrid/admin/fragments/exporters.py

@@ -0,0 +1,82 @@
+# (c) Nelen & Schuurmans.  GPL licensed, see LICENSE.rst.
+
+import os
+from collections import OrderedDict
+
+try:
+    from osgeo import ogr
+except ImportError:
+    ogr = None
+
+from threedigrid.admin import exporter_constants as const
+from threedigrid.geo_utils import get_spatial_reference
+from threedigrid.numpy_utils import reshape_flat_array
+from threedigrid.orm.base.exporters import BaseOgrExporter
+
+
+class FragmentsOgrExporter(BaseOgrExporter):
+    """
+    Exports to ogr formats. You need to set the driver explicitly
+    before calling save()
+    """
+
+    def __init__(self, fragments):
+        """
+        :param fragments: fragments.models.Fragments instance
+        """
+        self._fragments = fragments
+        self.supported_drivers = {
+            const.GEO_PACKAGE_DRIVER_NAME,
+            const.SHP_DRIVER_NAME,
+            const.GEOJSON_DRIVER_NAME,
+        }
+        self.driver = None
+
+    def save(self, file_name, layer_name, target_epsg_code, **kwargs):
+        """
+        save to file format specified by the driver, e.g. shapefile
+
+        :param file_name: name of the outputfile
+        :param fragment_data: dict of fragment data
+        """
+        if self.driver is None:
+            self.set_driver(extension=os.path.splitext(file_name)[1])
+
+        assert self.driver is not None
+
+        sr = get_spatial_reference(target_epsg_code)
+
+        self.del_datasource(file_name)
+        data_source = self.driver.CreateDataSource(file_name)
+        layer = data_source.CreateLayer(
+            str(os.path.basename(file_name)), sr, ogr.wkbPolygon
+        )
+
+        fields = OrderedDict(
+            [
+                ("id", "int"),
+                ("node_id", "int"),
+            ]
+        )
+
+        for field_name, field_type in fields.items():
+            layer.CreateField(
+                ogr.FieldDefn(field_name, const.OGR_FIELD_TYPE_MAP[field_type])
+            )
+        _definition = layer.GetLayerDefn()
+
+        for i in range(len(self._fragments.coords)):
+            if self._fragments.id[i] == 0:
+                continue  # skip the dummy element
+            feature = ogr.Feature(_definition)
+            ring = ogr.Geometry(ogr.wkbLinearRing)
+            polygon_points = reshape_flat_array(self._fragments.coords[i]).T
+            for x in polygon_points:
+                ring.AddPoint(x[0], x[1])
+            polygon = ogr.Geometry(ogr.wkbPolygon)
+            polygon.AddGeometry(ring)
+            feature.SetGeometry(polygon)
+            self.set_field(feature, "id", "int", self._fragments.id[i])
+            self.set_field(feature, "node_id", "int", self._fragments.node_id[i])
+            layer.CreateFeature(feature)
+            feature.Destroy()

threedigrid/admin/fragments/models.py

@@ -0,0 +1,22 @@
+from threedigrid.admin.fragments import exporters
+from threedigrid.orm.fields import ArrayField, PolygonArrayField
+from threedigrid.orm.models import Model
+
+
+class Fragments(Model):
+
+    node_id = ArrayField(type=int)
+    coords = PolygonArrayField()
+
+    GPKG_DEFAULT_FIELD_MAP = {
+        "id": "id",
+        "node_id": "node_id",
+        "coords": "the_geom",
+    }
+
+    def __init__(self, *args, **kwargs):
+        super().__init__(*args, **kwargs)
+
+        self._exporters = [
+            exporters.FragmentsOgrExporter(self),
+        ]

threedigrid/admin/gridadmin.py

@@ -10,6 +10,7 @@
 
 from threedigrid.admin.breaches.models import Breaches
 from threedigrid.admin.crosssections.models import CrossSections
+from threedigrid.admin.fragments.models import Fragments
 from threedigrid.admin.h5py_datasource import H5pyGroup
 from threedigrid.admin.levees.models import Levees
 from threedigrid.admin.lines.models import Lines
@@ -115,6 +116,12 @@ def lines(self):
             self.datasource_class(self.h5py_file, "lines"), **self._grid_kwargs
         )
 
+    @property
+    def fragments(self):
+        return Fragments(
+            self.datasource_class(self.h5py_file, "fragments"), **self._grid_kwargs
+        )
+
     @property
     def cross_sections(self):
         return CrossSections(

threedigrid/admin/prepare.py

@@ -508,6 +508,20 @@ def export_pumps(self):
             dest, indent=self._indent
         )
 
+    def export_fragments(self):
+        if not self.ga.fragments.id.size == 0:
+            logger.info(
+                "[*] Model {} does not have fragments, "
+                "skipping export fragments...".format(self.ga.model_name)
+            )
+            return
+
+        dest_ow = os.path.join(self._dest, "fragments" + self._extension)
+        getattr(
+            self.ga.fragments.reproject_to(self._epsg),
+            self._export_method,
+        )(dest_ow, indent=self._indent)
+
     def export_levees(self):
         """
         writes shapefile of levees