From 895afc232b1ff0c1ad7bb78f056a552c78efe816 Mon Sep 17 00:00:00 2001 From: Ben van Basten Date: Fri, 13 Dec 2024 16:35:05 +0100 Subject: [PATCH 1/4] stub fragment export --- README.rst | 24 ++++++ tests/conftest.py | 8 ++ tests/test_exporter.py | 13 +++ .../fragments/gridadmin_fragments.h5 | 3 + .../gridadmin_fuerthen.h5:Zone.Identifier | 2 + threedigrid/admin/fragments/_init__.py | 0 threedigrid/admin/fragments/exporters.py | 82 +++++++++++++++++++ threedigrid/admin/fragments/models.py | 16 ++++ threedigrid/admin/gridadmin.py | 7 ++ 9 files changed, 155 insertions(+) create mode 100644 tests/test_files/fragments/gridadmin_fragments.h5 create mode 100644 tests/test_files/fragments/gridadmin_fuerthen.h5:Zone.Identifier create mode 100644 threedigrid/admin/fragments/_init__.py create mode 100644 threedigrid/admin/fragments/exporters.py create mode 100644 threedigrid/admin/fragments/models.py diff --git a/README.rst b/README.rst index a98fc43f..fe7af1f8 100644 --- a/README.rst +++ b/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 ------- diff --git a/tests/conftest.py b/tests/conftest.py index d40bb789..d59f0a89 100644 --- a/tests/conftest.py +++ b/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(): """ diff --git a/tests/test_exporter.py b/tests/test_exporter.py index fd4568f4..1f0cee53 100644 --- a/tests/test_exporter.py +++ b/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) diff --git a/tests/test_files/fragments/gridadmin_fragments.h5 b/tests/test_files/fragments/gridadmin_fragments.h5 new file mode 100644 index 00000000..49ad9228 --- /dev/null +++ b/tests/test_files/fragments/gridadmin_fragments.h5 @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:251535f6f199b4e808e0cd61c329f8dd4242dac7cdd273c6c32496c7895fa0f1 +size 1236821 diff --git a/tests/test_files/fragments/gridadmin_fuerthen.h5:Zone.Identifier b/tests/test_files/fragments/gridadmin_fuerthen.h5:Zone.Identifier new file mode 100644 index 00000000..a45e1ac4 --- /dev/null +++ b/tests/test_files/fragments/gridadmin_fuerthen.h5:Zone.Identifier @@ -0,0 +1,2 @@ +[ZoneTransfer] +ZoneId=3 diff --git a/threedigrid/admin/fragments/_init__.py b/threedigrid/admin/fragments/_init__.py new file mode 100644 index 00000000..e69de29b diff --git a/threedigrid/admin/fragments/exporters.py b/threedigrid/admin/fragments/exporters.py new file mode 100644 index 00000000..f9121531 --- /dev/null +++ b/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() diff --git a/threedigrid/admin/fragments/models.py b/threedigrid/admin/fragments/models.py new file mode 100644 index 00000000..845ce883 --- /dev/null +++ b/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), + ] diff --git a/threedigrid/admin/gridadmin.py b/threedigrid/admin/gridadmin.py index 41d031a2..3c229469 100644 --- a/threedigrid/admin/gridadmin.py +++ b/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( From c3daeb381007f26bbbf4c9efb3616f5cf0c535d4 Mon Sep 17 00:00:00 2001 From: Ben <97883955+benvanbasten-ns@users.noreply.github.com> Date: Fri, 13 Dec 2024 16:36:46 +0100 Subject: [PATCH 2/4] Delete tests/test_files/fragments/gridadmin_fuerthen.h5:Zone.Identifier --- .../test_files/fragments/gridadmin_fuerthen.h5:Zone.Identifier | 2 -- 1 file changed, 2 deletions(-) delete mode 100644 tests/test_files/fragments/gridadmin_fuerthen.h5:Zone.Identifier diff --git a/tests/test_files/fragments/gridadmin_fuerthen.h5:Zone.Identifier b/tests/test_files/fragments/gridadmin_fuerthen.h5:Zone.Identifier deleted file mode 100644 index a45e1ac4..00000000 --- a/tests/test_files/fragments/gridadmin_fuerthen.h5:Zone.Identifier +++ /dev/null @@ -1,2 +0,0 @@ -[ZoneTransfer] -ZoneId=3 From 997fc1f42fa7f0bc4cef9cba8dc82889e985d89c Mon Sep 17 00:00:00 2001 From: Ben van Basten Date: Fri, 13 Dec 2024 17:20:53 +0100 Subject: [PATCH 3/4] Geopackage export --- tests/conftest.py | 11 +++++++---- tests/test_exporter.py | 6 ++++++ threedigrid/admin/exporter_constants.py | 1 + threedigrid/admin/exporters/geopackage/default.py | 3 +++ threedigrid/admin/exporters/geopackage/exporter.py | 7 +++++++ threedigrid/admin/fragments/models.py | 6 ++++++ 6 files changed, 30 insertions(+), 4 deletions(-) diff --git a/tests/conftest.py b/tests/conftest.py index d59f0a89..bde58d77 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -36,10 +36,13 @@ def ga(request): @pytest.fixture() -def ga_fragments(): - with GridH5Admin( - os.path.join(test_file_dir, "fragments", "gridadmin_fragments.h5") - ) as ga: +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 diff --git a/tests/test_exporter.py b/tests/test_exporter.py index 1f0cee53..b76f4cbe 100644 --- a/tests/test_exporter.py +++ b/tests/test_exporter.py @@ -4,6 +4,7 @@ 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 @@ -49,6 +50,11 @@ def test_fragments_gpgk_export(ga_fragments, tmp_path): ) # minus dummy +def test_fragments_gpkg_conversion_smoke(ga_fragment_path): + exporter = GeopackageExporter(ga_fragment_path, "fragment_export.gpkg") + exporter.export() + + def test_meta_data_gpgk_export(ga, tmp_path): path = str(tmp_path / ("exporter_meta.gpkg")) exporter = GpkgExporter(ga) diff --git a/threedigrid/admin/exporter_constants.py b/threedigrid/admin/exporter_constants.py index 35918750..7e968259 100644 --- a/threedigrid/admin/exporter_constants.py +++ b/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" diff --git a/threedigrid/admin/exporters/geopackage/default.py b/threedigrid/admin/exporters/geopackage/default.py index a81b5599..79f4bb5e 100644 --- a/threedigrid/admin/exporters/geopackage/default.py +++ b/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' diff --git a/threedigrid/admin/exporters/geopackage/exporter.py b/threedigrid/admin/exporters/geopackage/exporter.py index 685450af..1e5bfa3e 100644 --- a/threedigrid/admin/exporters/geopackage/exporter.py +++ b/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) diff --git a/threedigrid/admin/fragments/models.py b/threedigrid/admin/fragments/models.py index 845ce883..448c5598 100644 --- a/threedigrid/admin/fragments/models.py +++ b/threedigrid/admin/fragments/models.py @@ -8,6 +8,12 @@ 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) From 8b5b62af70ba190b07bccaab095b6121aee5fbbf Mon Sep 17 00:00:00 2001 From: Ben van Basten Date: Tue, 17 Dec 2024 11:47:11 +0100 Subject: [PATCH 4/4] Added tests --- tests/test_exporter.py | 51 ++++++++++++++++++++++++- tests/test_gridadmin.py | 26 +++++++++++++ threedigrid/admin/exporter_constants.py | 3 ++ threedigrid/admin/prepare.py | 14 +++++++ threedigrid/admin/serializers.py | 22 ++++++++++- 5 files changed, 113 insertions(+), 3 deletions(-) diff --git a/tests/test_exporter.py b/tests/test_exporter.py index b76f4cbe..644ef06f 100644 --- a/tests/test_exporter.py +++ b/tests/test_exporter.py @@ -1,5 +1,6 @@ import json import os +from pathlib import Path import pytest from osgeo import ogr @@ -50,9 +51,54 @@ def test_fragments_gpgk_export(ga_fragments, tmp_path): ) # minus dummy -def test_fragments_gpkg_conversion_smoke(ga_fragment_path): - exporter = GeopackageExporter(ga_fragment_path, "fragment_export.gpkg") +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): @@ -157,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): diff --git a/tests/test_gridadmin.py b/tests/test_gridadmin.py index d77eaa73..f8661b8c 100644 --- a/tests/test_gridadmin.py +++ b/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) diff --git a/threedigrid/admin/exporter_constants.py b/threedigrid/admin/exporter_constants.py index 7e968259..6335ee09 100644 --- a/threedigrid/admin/exporter_constants.py +++ b/threedigrid/admin/exporter_constants.py @@ -179,6 +179,8 @@ "z_coordinate", ] +FRAGMENTS_EXPORT_FIELDS = ["id", "node_id"] + DEFAULT_EXPORT_FIELDS = { "Lines": "ALL", "Pipes": PIPES_EXPORT_FIELDS, @@ -194,4 +196,5 @@ "Breaches": BREACHES_EXPORT_FIELDS, "Levees": LEVEES_EXPORT_FIELDS, "Pumps": PUMPS_EXPORT_FIELDS, + "Fragments": FRAGMENTS_EXPORT_FIELDS, } diff --git a/threedigrid/admin/prepare.py b/threedigrid/admin/prepare.py index bf0b0a0d..e915b559 100644 --- a/threedigrid/admin/prepare.py +++ b/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 diff --git a/threedigrid/admin/serializers.py b/threedigrid/admin/serializers.py index afbe9b3f..9ee7ac7e 100644 --- a/threedigrid/admin/serializers.py +++ b/threedigrid/admin/serializers.py @@ -5,7 +5,7 @@ import numpy as np from threedigrid.admin import constants -from threedigrid.geo_utils import raise_import_exception, transform_bbox +from threedigrid.geo_utils import raise_import_exception, transform_bbox, transform_xys from threedigrid.orm.base.encoder import NumpyEncoder from threedigrid.orm.base.models import Model @@ -127,6 +127,26 @@ def geos_iter(self): ) properties = fill_properties(self.fields, data, i, model_type) yield geojson.Feature(geometry=polygon, properties=properties) + elif content_type == "fragments": + for i in range(data["id"].shape[-1]): + coords = np.round( + data["coords"][i].reshape(2, -1).astype("float64"), + constants.LONLAT_DIGITS, + ) + if ( + self._model.reproject_to_epsg is not None + and self._model.reproject_to_epsg != self._model.epsg_code + ): + # Pick reproject_to_epsg or original model epsg_code + coords = transform_xys( + np.array(coords[0]), + np.array(coords[1]), + self._model.epsg_code, + self._model.reproject_to_epsg, + ) + polygon = geojson.Polygon(coords.T.tolist()) + properties = fill_properties(self.fields, data, i, model_type) + yield geojson.Feature(geometry=polygon, properties=properties) elif content_type == "levees": for i in range(data["id"].shape[-1]): coords = np.round(