Merge pull request #51 from fema-ffrd/bug/mesh_polygons

bug/mesh_cell_polygons

pyproject.toml

@@ -12,7 +12,7 @@ classifiers = [
     "Programming Language :: Python :: 3.11",
     "Programming Language :: Python :: 3.12",
 ]
-version = "0.2.1"
+version = "0.2.2"
 dependencies = ["h5py", "geopandas", "pyarrow"]
 
 [project.optional-dependencies]

src/rashdf/geom.py

@@ -18,7 +18,7 @@
     LineString,
     MultiLineString,
     MultiPolygon,
-    polygonize,
+    polygonize_full,
 )
 
 from typing import Dict, List, Optional
@@ -139,13 +139,19 @@ def mesh_cell_polygons(self) -> GeoDataFrame:
             cell_dict["cell_id"] += cell_ids
             cell_dict["geometry"] += list(
                 np.vectorize(
-                    lambda face_id_list: polygonize(
-                        np.ravel(
-                            mesh_faces[
-                                np.array(face_id_list.strip("[]").split()).astype(int)
-                            ]
+                    lambda face_id_list: (
+                        lambda geom_col: Polygon((geom_col[0] or geom_col[3]).geoms[0])
+                    )(
+                        polygonize_full(
+                            np.ravel(
+                                mesh_faces[
+                                    np.array(face_id_list.strip("[]").split()).astype(
+                                        int
+                                    )
+                                ]
+                            )
                         )
-                    ).geoms[0]
+                    )
                 )(face_id_lists)
             )
         return GeoDataFrame(cell_dict, geometry="geometry", crs=self.projection())

src/rashdf/utils.py

@@ -7,6 +7,7 @@
 from datetime import datetime, timedelta
 import re
 from typing import Any, List, Tuple, Union, Optional
+from shapely import LineString, Polygon, polygonize_full
 
 
 def parse_ras_datetime(datetime_str: str) -> datetime:

tests/data/json/coal-bad-mesh-cell-polygons.json

@@ -0,0 +1 @@
+{"type": "FeatureCollection", "features": [{"id": "8561", "type": "Feature", "properties": {"mesh_name": "Perimeter 1", "cell_id": 8561}, "geometry": {"type": "Polygon", "coordinates": [[[-81.93838041650048, 38.18701291114405], [-81.93351085972894, 38.18747285993424], [-81.93039764009161, 38.19022466753743], [-81.93142853728568, 38.19420939282696], [-81.93232042778128, 38.19513770535859], [-81.93243907136474, 38.19501410050651], [-81.93256378275053, 38.19486462780729], [-81.93267562214464, 38.19470913288595], [-81.93277409700462, 38.19454832876619], [-81.93285875991612, 38.194382878687314], [-81.93292925736498, 38.19421351170956], [-81.9329852757799, 38.194040934003525], [-81.9330266058494, 38.193865863766945], [-81.93305305747408, 38.19368910389221], [-81.93306452560488, 38.19351138459886], [-81.93306138417432, 38.193340722380135], [-81.93303400386331, 38.192876950415375], [-81.9330335544586, 38.19286969504248], [-81.93303252059805, 38.19285404144383], [-81.93300656475326, 38.192482814659584], [-81.93305844587206, 38.19239900605469], [-81.93311257114541, 38.19230723565101], [-81.93316649068662, 38.192205653525775], [-81.93339427107881, 38.19175032466181], [-81.93339433691649, 38.19175016729454], [-81.93345209957121, 38.19173742278804], [-81.93366861672425, 38.19167761870666], [-81.9338796246049, 38.19160695060574], [-81.93390735139756, 38.19159670579108], [-81.93413452496513, 38.19151181880671], [-81.93431141415032, 38.191440869306206], [-81.9345086816885, 38.19134950608091], [-81.93469782292014, 38.19124832417171], [-81.93487799927583, 38.191137804228454], [-81.93504848163037, 38.19101841203771], [-81.93520846619597, 38.19089063226441], [-81.93535734733395, 38.19075502742718], [-81.9354944447435, 38.19061217892142], [-81.93561126721173, 38.19047279197159], [-81.93567580284684, 38.190390388031624], [-81.93568371282369, 38.19038030312414], [-81.93579553880446, 38.190224805296864], [-81.9358939660921, 38.19006399459883], [-81.93597861682495, 38.18989854230042], [-81.93604910279267, 38.18972917346145], [-81.93610511047302, 38.189556594260694], [-81.9361171994328, 38.189511716540366], [-81.93618882069218, 38.18923335433667], [-81.9363174447124, 38.18911861646767], [-81.93641436568782, 38.18902831931168], [-81.93651872114977, 38.18892136619254], [-81.93680590274592, 38.18861151423493], [-81.93683860439668, 38.18857561355812], [-81.93696333585468, 38.18842611341469], [-81.93705621426075, 38.1882986998], [-81.93725476928356, 38.18800839708326], [-81.93727367110456, 38.18798033453258], [-81.93737212690529, 38.18781952662532], [-81.93745677173061, 38.18765407327319], [-81.93752725214433, 38.187484703550176], [-81.93757904878805, 38.18732511046833], [-81.93763294391364, 38.18729627407441], [-81.93781310617955, 38.187185749783346], [-81.93792607873903, 38.18710835661449], [-81.93809919940769, 38.18698459220709], [-81.9381593991755, 38.186993105191426], [-81.93838041650048, 38.18701291114405]]]}}, {"id": "11791", "type": "Feature", "properties": {"mesh_name": "Perimeter 1", "cell_id": 11791}, "geometry": {"type": "Polygon", "coordinates": [[[-81.91297751133456, 38.23134252029514], [-81.91188679329382, 38.23048260485189], [-81.91045769203454, 38.23392751638167], [-81.91046919604935, 38.23389991579975], [-81.91051388163784, 38.233766777848984], [-81.9105234610186, 38.23373483540925], [-81.9105733817855, 38.233741922785526], [-81.91080188008674, 38.23376241066241], [-81.91103161703967, 38.23377132826821], [-81.9111503572862, 38.23377138132072], [-81.9114397520803, 38.23376769636856], [-81.91168412358901, 38.23376653552678], [-81.91183155068104, 38.23376344063458], [-81.91206084240088, 38.23374907150056], [-81.91228842431191, 38.233723121460486], [-81.91251334561395, 38.23368572815638], [-81.91259687745733, 38.23366876239905], [-81.91297978037437, 38.23358709469575], [-81.91297751133456, 38.23134252029514]]]}}, {"id": "17529", "type": "Feature", "properties": {"mesh_name": "Perimeter 1", "cell_id": 17529}, "geometry": {"type": "Polygon", "coordinates": [[[-81.73490117509893, 37.81697512194823], [-81.73287781741101, 37.81759694733074], [-81.72764230259725, 37.815160872140794], [-81.72714478327165, 37.81114739314949], [-81.7271201031621, 37.811310097062474], [-81.72710515815339, 37.811380757602805], [-81.72705835547083, 37.811584861210534], [-81.72704530679083, 37.8116417725805], [-81.72707213544749, 37.8116499648186], [-81.72720204526726, 37.811692517825634], [-81.72724569704283, 37.81170709155368], [-81.72745221137566, 37.81178363802989], [-81.72765179914522, 37.81187053554256], [-81.72784364871796, 37.81196741351272], [-81.72802699346478, 37.8120738516558], [-81.7282009735213, 37.81218936370556], [-81.728364862016, 37.81231350656279], [-81.72842497377707, 37.81236541625037], [-81.72853838860757, 37.812416449949055], [-81.72866555956448, 37.81247364127586], [-81.72867287011745, 37.81247694984793], [-81.72867768167255, 37.81247911157303], [-81.72868882450975, 37.81248416403009], [-81.72875973803754, 37.812516524074574], [-81.72876022505966, 37.81251674643825], [-81.72912177431314, 37.81268173323116], [-81.72928444582216, 37.812755975230054], [-81.72930350604629, 37.812737948048614], [-81.72941815762739, 37.81263501747653], [-81.72957771302501, 37.81250750277632], [-81.72974761891639, 37.812388385097044], [-81.72992721254754, 37.812278163976636], [-81.73011566361012, 37.81217729172387], [-81.73031217114288, 37.81208625160828], [-81.73051592578337, 37.81200538838854], [-81.73072603860624, 37.811935092469625], [-81.73094164684169, 37.811875619815495], [-81.73099635020131, 37.81186244806439], [-81.73119910054578, 37.81181510382175], [-81.73136449440578, 37.81177989249515], [-81.7315881886892, 37.81174282419043], [-81.73181452937332, 37.81171723386156], [-81.73204254685436, 37.81170319989176], [-81.73206663905889, 37.8117023774302], [-81.73264672381931, 37.811684645051514], [-81.73266243632429, 37.81168110343669], [-81.73324548180697, 37.81153994623256], [-81.73328035634157, 37.81148057011723], [-81.73335465229329, 37.81136188541789], [-81.73346632508355, 37.811206590266664], [-81.73359079631211, 37.81105728452808], [-81.73372757919985, 37.81091465350389], [-81.73379379717831, 37.81085206239745], [-81.73391008562447, 37.81074538813885], [-81.7339923526406, 37.81067261121904], [-81.73415189591516, 37.81054509049176], [-81.73420314327197, 37.810507499325055], [-81.73439744165873, 37.81036775830114], [-81.73451608787656, 37.81028622518639], [-81.73455154508562, 37.810264456388], [-81.73545583108775, 37.81112954077671], [-81.73490117509893, 37.81697512194823]]]}}]}

tests/test_geom.py

@@ -9,6 +9,7 @@
 
 TEST_DATA = Path("./tests/data")
 MUNCIE_G05 = TEST_DATA / "ras/Muncie.g05.hdf"
+COAL_G01 = TEST_DATA / "ras/Coal.g01.hdf"
 TEST_JSON = TEST_DATA / "json"
 
 TEST_ATTRS = {"test_attribute1": "test_str1", "test_attribute2": 500}
@@ -59,6 +60,24 @@ def test_mesh_cell_polygons():
         assert _gdf_matches_json(ghdf.mesh_cell_polygons(), mesh_cell_polygons_json)
 
 
+def test_mesh_cell_polygons_coal():
+    """Test with the mesh from the Coal River model.
+
+    The Jan 2024 Coal River model from the Kanawha FFRD pilot project
+    contains some topologically incorrect polygons in the 2D mesh;
+    some of the mesh cell faces overlap.
+
+    See: https://github.com/fema-ffrd/rashdf/issues/31
+    """
+    coal_bad_polygons_json = TEST_JSON / "coal-bad-mesh-cell-polygons.json"
+    with RasGeomHdf(COAL_G01) as geom_hdf:
+        gdf = geom_hdf.mesh_cell_polygons()
+        gdf_bad_polygons = gdf[gdf["cell_id"].isin([8561, 11791, 17529])].to_crs(
+            "EPSG:4326"
+        )  # reproject because the Kanawha FFRD pilot project uses a custom CRS
+        assert _gdf_matches_json(gdf_bad_polygons, coal_bad_polygons_json)
+
+
 def test_bc_lines():
     bc_lines_json = TEST_JSON / "bc_lines.json"
     with RasGeomHdf(MUNCIE_G05) as ghdf: