remove pygeos & work on mesh verification

Since pygeos is merged into shapely 2 we remove it as dependency. Also took the opportunity to work out the mesh verification process. The commit resolves ec-jrc#23 & ec-jrc#60

pyposeidon/mesh.py

@@ -593,9 +593,10 @@ def validate(self, **kwargs):
     def verify(self, **kwargs):
 
         shp = kwargs.get("coastlines", None)
+        thorough = kwargs.get("thorough", False)
 
         if shp is not None:
-            r = verify(self, shp)
+            r = verify(self, shp, thorough)
             return r
         else:
             logger.warning("No coastlines provided")

pyposeidon/utils/fix.py

@@ -11,7 +11,6 @@
 import numpy as np
 import geopandas as gp
 import shapely
-import pygeos
 import pyresample
 import pandas as pd
 import xarray as xr
@@ -95,10 +94,7 @@ def fix(dem, coastline, **kwargs):
 
     g = block.unary_union.symmetric_difference(grp)  # get the diff
 
-    try:
-        t = gp.GeoDataFrame({"geometry": g})
-    except:
-        t = gp.GeoDataFrame({"geometry": [g]})
+    t = gp.GeoDataFrame({"geometry": [g]}).explode(index_parts=True).droplevel(0)
 
     t["length"] = t["geometry"][:].length  # optional
 
@@ -145,29 +141,31 @@ def fix(dem, coastline, **kwargs):
         df = dem.elevation.to_dataframe().reset_index()
 
     # ---------------------------------------------------------------------
-    logger.debug("invoke pygeos\n")
+    logger.debug("invoke shapely\n")
     # ---------------------------------------------------------------------
 
-    spoints_ = pygeos.points(list(df.loc[:, ["longitude", "latitude"]].values))  # create pygeos objects for the points
+    spoints_ = shapely.points(
+        list(df.loc[:, ["longitude", "latitude"]].values)
+    )  # create shapely objects for the points
 
-    # Add land boundaries to a pygeos object
+    # Add land boundaries to a shapely object
     try:
         lbs = []
         for l in range(len(land.boundary.geoms)):
-            z = pygeos.linearrings(land.boundary.geoms[l].coords[:])
+            z = shapely.linearrings(land.boundary.geoms[l].coords[:])
             lbs.append(z)
     except:
-        lbs = pygeos.linearrings(land.boundary.coords[:])
+        lbs = shapely.linearrings(land.boundary.coords[:])
 
-    bp = pygeos.polygons(lbs)
+    bp = shapely.polygons(lbs)
 
     # ---------------------------------------------------------------------
     logger.debug("find wet and dry masks\n")
     # ---------------------------------------------------------------------
 
     # find the points on land
 
-    tree = pygeos.STRtree(spoints_)
+    tree = shapely.STRtree(spoints_)
 
     try:
         wl = []

pyposeidon/utils/hplot.py

@@ -19,7 +19,6 @@
 import cartopy.feature as cfeature
 import xarray as xr
 import pandas as pd
-import pygeos
 import spatialpandas
 import geopandas as gp
 from hvplot import xarray
@@ -178,7 +177,7 @@ def mesh(self, tiles=False, **kwargs):
 
             quads["coordinates"] = coords
 
-            qpolys = pygeos.polygons(quads.coordinates.to_list())
+            qpolys = shapely.polygons(quads.coordinates.to_list())
 
             df = gp.GeoDataFrame(geometry=qpolys)
 

pyposeidon/utils/postgrid.py

@@ -6,7 +6,6 @@
 import numpy as np
 import pandas as pd
 import geopandas as gp
-import pygeos
 import xarray as xr
 from tqdm import tqdm
 import shapely
@@ -175,12 +174,12 @@ def check(g, shp, bad):
         nodes, elems, bnodes = drop(nodes, elems, bnodes, dq)
 
     # ### Find the invalid nodes (that cross the coasts)
-    cos = pygeos.from_shapely(c.geometry)
-    cos_ = pygeos.set_operations.union_all(cos)
+    cos = c.geometry
+    cos_ = shapely.set_operations.union_all(cos)
 
-    gps = pygeos.points(list(nodes.values))
+    gps = shapely.points(list(nodes.values))
 
-    gtree = pygeos.STRtree(gps)
+    gtree = shapely.STRtree(gps)
 
     invs = gtree.query(cos_, predicate="contains").tolist()
 
@@ -200,26 +199,20 @@ def check(g, shp, bad):
     coords = [[l[i : i + n] for i in range(0, len(l), n)] for l in al]
     elems["coordinates"] = coords
 
-    jig = pygeos.polygons(coords)
+    jig = shapely.polygons(coords)
 
-    jtree = pygeos.STRtree(jig)
+    jtree = shapely.STRtree(jig)
 
-    jig_ = pygeos.set_operations.union_all(jig)
+    jig_ = shapely.set_operations.union_all(jig)
 
-    cross = pygeos.set_operations.intersection(jig_, cos_)
+    cross = shapely.set_operations.intersection(jig_, cos_)
 
     # #### convert to dataframe
 
-    fd = pd.DataFrame({"overlap": pygeos.to_wkt(cross)}, index=[0])
-
-    fd["overlap"] = fd["overlap"].apply(shapely.wkt.loads)
-
-    gover = gp.GeoDataFrame(fd, geometry="overlap")
+    fd = gp.GeoDataFrame({"geometry": cross}, index=[0])
 
     # #### Reject small injuctions
-    ipols = gover.explode(index_parts=True).loc[0]
-
-    ipols.columns = ["geometry"]
+    ipols = fd.explode(index_parts=True).loc[0]
 
     mask = ipols.area.values == 0.0
 
@@ -230,7 +223,7 @@ def check(g, shp, bad):
 
     # ## Get overlapping elements for each contour
 
-    jcos = pygeos.from_shapely(ipols.geometry)
+    jcos = ipols.geometry
 
     maxb = bnodes.id.min()
 
@@ -239,17 +232,17 @@ def check(g, shp, bad):
         for l in tqdm(range(len(jcos))):
             con = jcos[l]
             m += 1
-            jig = pygeos.polygons(elems.coordinates.to_list())
-            jtree = pygeos.STRtree(jig)
+            jig = shapely.polygons(elems.coordinates.to_list())
+            jtree = shapely.STRtree(jig)
             ci = jtree.query(con, predicate="intersects").tolist()
             if not ci:
                 continue
             delems = elems.loc[ci].index.values
             dnodes = np.unique(elems.loc[ci, ["a", "b", "c"]].values.flatten())
             #    print (ci, dnodes, delems)
 
-            inp = pygeos.points(list(nodes.loc[dnodes].values))
-            itree = pygeos.STRtree(inp)
+            inp = shapely.points(list(nodes.loc[dnodes].values))
+            itree = shapely.STRtree(inp)
             ipoints = itree.query(cos_, predicate="contains").tolist()
             ipoints = nodes.loc[dnodes].index[ipoints].values
 
@@ -328,9 +321,9 @@ def check(g, shp, bad):
         nodes, elems, bnodes = drop(nodes, elems, bnodes, dq)
 
     # Get the largest continous area
-    jels = pygeos.polygons(elems.coordinates.values.tolist())
-    wat = pygeos.set_operations.coverage_union_all(jels)
-    w = pd.DataFrame({"overlap": pygeos.to_wkt(wat)}, index=[0])
+    jels = shapely.polygons(elems.coordinates.values.tolist())
+    wat = shapely.set_operations.coverage_union_all(jels)
+    w = pd.DataFrame({"overlap": shapely.to_wkt(wat)}, index=[0])
     w["overlap"] = w["overlap"].apply(shapely.wkt.loads)
     gw = gp.GeoDataFrame(w, geometry="overlap")
 
@@ -342,18 +335,18 @@ def check(g, shp, bad):
     gw = gw.reset_index(drop=True)
 
     # indentify the elements of the large polygon
-    cgw = pygeos.from_shapely(gw.loc[1:].geometry)
-    cgw_ = pygeos.set_operations.union_all(cgw)
-    jtree_ = pygeos.STRtree(jels)
+    cgw = shapely.from_shapely(gw.loc[1:].geometry)
+    cgw_ = shapely.set_operations.union_all(cgw)
+    jtree_ = shapely.STRtree(jels)
     invs = jtree_.query(cgw_, predicate="intersects").tolist()
 
     if len(invs) > 0:
         # Sort the elements (some shouldn't be there)
 
         qnodes = np.unique([elems.loc[x, ["a", "b", "c"]].values.astype(int) for x in invs])
-        nnodes = pygeos.points(list(nodes.loc[qnodes].values))
-        ntree = pygeos.STRtree(nnodes)
-        nels_ = pygeos.from_shapely(gw.loc[0].geometry.buffer(0.00001))
+        nnodes = shapely.points(list(nodes.loc[qnodes].values))
+        ntree = shapely.STRtree(nnodes)
+        nels_ = gw.loc[0].geometry.buffer(0.00001)
         nevs = ntree.query(nels_, predicate="intersects").tolist()
         pns = qnodes[nevs]
         dpoints = [x for x in qnodes if x not in pns]

pyposeidon/utils/spline.py

@@ -2,7 +2,6 @@
 import geopandas as gp
 import numpy as np
 import shapely
-import pygeos
 from scipy.interpolate import interp1d
 
 # based on https://stackoverflow.com/questions/52014197/how-to-interpolate-a-2d-curve-in-python

pyposeidon/utils/verify.py

@@ -6,7 +6,6 @@
 import numpy as np
 import pandas as pd
 import geopandas as gp
-import pygeos
 import xarray as xr
 from tqdm import tqdm
 import shapely
@@ -18,7 +17,7 @@
 
 def verify(g, shp, thorough=False):
     # ---------------------------------------------------------------------
-    logger.info(" Verify grid against coastline\n")
+    logger.info(" Verify mesh against coastline\n")
     # ---------------------------------------------------------------------
 
     lon_min = g.Dataset.SCHISM_hgrid_node_x.values.min()
@@ -43,17 +42,17 @@ def verify(g, shp, thorough=False):
     bnodes = g.Dataset[["node", "id", "type"]].to_dataframe()
 
     # ### Find the invalid nodes (that cross the coasts)
-    cos = pygeos.from_shapely(c.geometry)
-    cos_ = pygeos.set_operations.union_all(cos)
+    cos = c.geometry
+    cos_ = shapely.set_operations.union_all(cos)
 
-    gps = pygeos.points(list(nodes.values))
+    gps = shapely.points(list(nodes.values))
 
-    gtree = pygeos.STRtree(gps)
+    gtree = shapely.STRtree(gps)
 
     invs = gtree.query(cos_, predicate="contains").tolist()
 
     # ---------------------------------------------------------------------
-    logger.info("Number of nodes within the coastlines {}\n".format(len(invs)))
+    logger.info("Number of nodes within/touching the coastlines {}\n".format(len(invs)))
     # ---------------------------------------------------------------------
 
     nps = len(invs)
@@ -78,27 +77,21 @@ def verify(g, shp, thorough=False):
         coords = [[l[i : i + n] for i in range(0, len(l), n)] for l in al]
         elems["coordinates"] = coords
 
-        jig = pygeos.polygons(coords)
+        jig = shapely.polygons(coords)
 
-        jtree = pygeos.STRtree(jig)
+        jtree = shapely.STRtree(jig)
 
-        jig_ = pygeos.set_operations.union_all(jig)
+        jig_ = shapely.set_operations.union_all(jig)
 
-        cross = pygeos.set_operations.intersection(jig_, cos_)
+        cross = shapely.set_operations.intersection(jig_, cos_)
 
         # #### convert to dataframe
 
-        fd = pd.DataFrame({"overlap": pygeos.to_wkt(cross)}, index=[0])
+        fd = gp.GeoDataFrame({"geometry": cross}, index=[0])
 
-        fd["overlap"] = fd["overlap"].apply(shapely.wkt.loads)
-
-        gover = gp.GeoDataFrame(fd, geometry="overlap")
+        ipols = fd.explode(index_parts=True).loc[0]
 
         # #### Reject small injuctions
-        ipols = gover.explode(index_parts=True).loc[0]
-
-        ipols.columns = ["geometry"]
-
         mask = ipols.area.values == 0.0
 
         ipols = ipols[~mask].reset_index(drop=True)
@@ -107,21 +100,29 @@ def verify(g, shp, thorough=False):
         # ---------------------------------------------------------------------
         logger.info("Number of elements intersecting the coastlines {}\n".format(ipols.shape[0]))
         # ---------------------------------------------------------------------
+        logger.info("Maximum intersecting area {}\n".format(ipols.area.max()))
+        # ---------------------------------------------------------------------
 
         nels = ipols.shape[0]
 
+        if ipols.area.max() < 1e-3:
+            # ---------------------------------------------------------------------
+            logger.info("Mesh is verified against the coastline")
+            # ---------------------------------------------------------------------
+            return True
+
     if nps == 0 and nels == 0:
         # ---------------------------------------------------------------------
-        logger.info("Grid is verified against the coastline")
+        logger.info("Mesh is verified against the coastline")
         # ---------------------------------------------------------------------
         return True
     elif nps == 0:
         # ---------------------------------------------------------------------
-        logger.info("Grid is node verified against the coastline")
+        logger.info("Mesh is node verified against the coastline")
         # ---------------------------------------------------------------------
         return True
     else:
         # ---------------------------------------------------------------------
-        logger.warning("Grid is not verified against the coastline")
+        logger.warning("Mesh can't be verified against the coastline")
         # ---------------------------------------------------------------------
         return False