stub fragment export

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,14 @@ def ga(request):
         yield ga
 
 
+@pytest.fixture()
+def ga_fragments():
+    with GridH5Admin(
+        os.path.join(test_file_dir, "fragments", "gridadmin_fragments.h5")
+    ) as ga:
+        yield ga
+
+
 @pytest.fixture
 def gr_bergen_with_boundaries():
     """

tests/test_exporter.py

@@ -5,6 +5,7 @@
 from osgeo import ogr
 
 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 +37,18 @@ 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_meta_data_gpgk_export(ga, tmp_path):
     path = str(tmp_path / ("exporter_meta.gpkg"))
     exporter = GpkgExporter(ga)

tests/test_files/fragments/gridadmin_fuerthen.h5:Zone.Identifier

@@ -0,0 +1,2 @@
+[ZoneTransfer]
+ZoneId=3

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,16 @@
+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()
+
+    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(