Merge pull request #82 from noaa-ocs-modeling/depend/shapely2

Add Support for Shapely 2

.github/workflows/documentation.yml

@@ -16,7 +16,7 @@ jobs:
 
     steps:
     - name: Checkout repository
-      uses: actions/checkout@v2
+      uses: actions/checkout@v3
     - name: Setup Mamba environment for Python
       uses: mamba-org/provision-with-micromamba@main
       with:

.github/workflows/functional_test.yml

@@ -15,7 +15,7 @@ jobs:
         python-version: [ '3.8', '3.9', '3.10' ]
 
     steps:
-    - uses: actions/checkout@v2
+    - uses: actions/checkout@v3
 
     - name: Setup Mamba environment for Python
       uses: mamba-org/provision-with-micromamba@main

.github/workflows/functional_test_2.yml

@@ -16,15 +16,15 @@ jobs:
 
     steps:
     - name: Checkout 
-      uses: actions/checkout@v2
+      uses: actions/checkout@v3
 
     - name: OS binaries
       shell: bash -l {0}
       run: |
         sudo apt install gdal-bin
 
     - name: Install Python
-      uses: actions/setup-python@v2
+      uses: actions/setup-python@v4
       with:
         python-version: ${{ matrix.python-version }}
 

.github/workflows/pylint.yml

@@ -8,7 +8,7 @@ jobs:
     runs-on: ubuntu-latest
 
     steps:
-    - uses: actions/checkout@v2
+    - uses: actions/checkout@v3
     - name: Setup Mamba environment for Python
       uses: mamba-org/provision-with-micromamba@main
       with:

.github/workflows/release.yml

@@ -15,9 +15,9 @@ jobs:
         python-version: [ '3.10' ]
     steps:
       - name: checkout repository
-        uses: actions/checkout@v2
+        uses: actions/checkout@v3
       - name: install Python
-        uses: actions/setup-python@v2
+        uses: actions/setup-python@v4
         with:
           python-version: ${{ matrix.python-version }}
       - name: load cached Python installation

environment.yml

@@ -2,17 +2,16 @@ name: ocsmesh
 channels:
   - conda-forge
 dependencies:
-  - python>=3.8, <3.11 # 3.11 not supported by numba
+  - python>=3.8, <3.11 # for numba
   - gdal
   - geos
   - proj
   - netcdf4
   - udunits2
   - pyproj
-  - shapely>=1.8, <2
+  - shapely
   - rasterio
   - fiona
-  - pygeos
   - geopandas
   - utm
   - scipy

ocsmesh/geom/collector.py

@@ -213,7 +213,8 @@ def __init__(
             elif isinstance(in_item, BaseMesh):
                 geom = MeshGeom(in_item)
 
-            elif isinstance(in_item, str):
+            elif isinstance(in_item, (Path, str)):
+                in_item = str(in_item)
                 if in_item.endswith('.tif'):
                     raster = Raster(in_item)
                     if self._base_shape:
@@ -306,8 +307,10 @@ def get_multipolygon(self, **kwargs: Any) -> MultiPolygon:
                 # TODO: Make sure all calcs are in EPSG:4326
                 base_multipoly = self._base_mesh.hull.multipolygon()
 
-            feather_files.append(self._extract_global_boundary(
-                temp_path, base_multipoly))
+            raster_geom_path = self._extract_global_boundary(
+                temp_path, base_multipoly)
+            if raster_geom_path.is_file():
+                feather_files.append(raster_geom_path)
             feather_files.extend(self._extract_nonraster_boundary(
                 temp_path, base_multipoly))
             feather_files.extend(self._extract_features(
@@ -447,11 +450,12 @@ def _type_chk(self, input_list) -> None:
             If the any of the inputs are not supported
         """
 
-        valid_types = (str, Raster, BaseGeom, BaseMesh)
+        valid_types = (str, Path, Raster, BaseGeom, BaseMesh)
         if not all(isinstance(item, valid_types) for item in input_list):
             raise TypeError(
-                f'Input list items must be of type {", ".join(valid_types)}'
-                f', or a derived type.')
+                f'Input list items must be of type'
+                f' {", ".join(str(i) for i in valid_types)},'
+                f' or a derived type.')
 
     def _get_raster_sources(self) -> List[Union[RasterGeom, Raster]]:
         """Get the list rasters from inputs to the object constructor
@@ -629,7 +633,7 @@ def _extract_nonraster_boundary(
             multipoly = self._get_valid_multipolygon(
                     geom.get_multipolygon())
             gdf_non_raster = gpd.GeoDataFrame(
-                    {'geometry': multipoly}, crs=crs)
+                    {'geometry': multipoly.geoms}, crs=crs)
             if crs != CRS.from_user_input("EPSG:4326"):
                 gdf_non_raster = gdf_non_raster.to_crs("EPSG:4326")
 
@@ -723,7 +727,7 @@ def _apply_patch(
             if ptch_defn:
                 patch_mp, crs = ptch_defn.get_multipolygon()
                 gdf_patch = gpd.GeoDataFrame(
-                        {'geometry': patch_mp}, crs=crs)
+                        {'geometry': patch_mp.geoms}, crs=crs)
                 if crs != CRS.from_user_input("EPSG:4326"):
                     gdf_patch = gdf_patch.to_crs("EPSG:4326")
                 combine_poly = MultiPolygon(list(gdf_patch.geometry))

ocsmesh/hfun/collector.py

@@ -1285,7 +1285,7 @@ def add_patch(
     def add_feature(
             self,
             shape: Union[MultiLineString, LineString, None] = None,
-            line_defn: Optional[LineString] = None,
+            line_defn: Optional[LineFeature] = None,
             shapefile: Union[None, str, Path] = None,
             expansion_rate: float = 0.01,
             target_size: Optional[float] = None,

ocsmesh/hfun/raster.py

@@ -836,9 +836,9 @@ def add_patch(
 
         # For expansion_rate
         if expansion_rate is not None:
-            exteriors = [ply.exterior for ply in multipolygon]
+            exteriors = [ply.exterior for ply in multipolygon.geoms]
             interiors = [
-                inter for ply in multipolygon for inter in ply.interiors]
+                inter for ply in multipolygon.geoms for inter in ply.interiors]
 
             features = MultiLineString([*exteriors, *interiors])
             # pylint: disable=E1123, E1125
@@ -862,7 +862,7 @@ def add_patch(
                 start = time()
                 try:
                     mask, _, _ = rasterio.mask.raster_geometry_mask(
-                        self.src, multipolygon,
+                        self.src, multipolygon.geoms,
                         all_touched=True, invert=True)
                     mask = mask[rasterio.windows.window_index(window)]
 

ocsmesh/mesh/mesh.py

@@ -745,6 +745,7 @@ def __call__(self) -> gpd.GeoDataFrame:
         data = []
         bnd_id = 0
         for poly in polys.geoms:
+            # TODO: Enforce CCW vs CW here
             data.append({
                     "geometry": poly.exterior,
                     "bnd_id": bnd_id,
@@ -1034,7 +1035,7 @@ def implode(self) -> gpd.GeoDataFrame:
 
         return gpd.GeoDataFrame(
             {"geometry": MultiPolygon([polygon.geometry for polygon
-                                       in self().itertuples()])},
+                                       in self().itertuples()]).geoms},
             crs=self.mesh.crs)
 
     def multipolygon(self) -> MultiPolygon:
@@ -1596,6 +1597,7 @@ def nodeidxlist(self):
 
         ext_bdry_nodes = []
         for ring in self.mesh.hull.rings.exterior().itertuples():
+            # TODO: Enforce external rings to be ccw as https://github.com/noaa-ocs-modeling/OCSMesh/issues/65
             ext_ring_coo = ring.geometry.coords
             ext_ring = np.array([
                     (coo_to_idx[ext_ring_coo[e]],
@@ -1607,6 +1609,7 @@ def nodeidxlist(self):
 
         int_bdry_nodes = []
         for ring in self.mesh.hull.rings.interior().itertuples():
+            # TODO: Internal rings should be cw?
             if not LinearRing(ring.geometry).is_ccw:
                 ring.geometry = ring.geometry.reverse()
             int_ring_coo = ring.geometry.coords

ocsmesh/ops/combine_geom.py

@@ -229,7 +229,7 @@ def run(self):
                     [
                         rasters_gdf,
                         gpd.GeoDataFrame(
-                            {'geometry': self._read_multipolygon(feather_f)},
+                            {'geometry': self._read_multipolygon(feather_f).geoms},
                             crs=self._calc_crs
                         ),
                     ],
@@ -289,7 +289,7 @@ def _multipolygon_to_disk(
             # a multipolygon instead of polygon for dataframe creation
             multipolygon = MultiPolygon([multipolygon])
 
-        gpd.GeoDataFrame({'geometry': multipolygon}).to_feather(path)
+        gpd.GeoDataFrame({'geometry': multipolygon.geoms}).to_feather(path)
 
 
     def _read_multipolygon(
@@ -523,7 +523,7 @@ def _write_to_file(
         # TODO: Check for correct extension on out_file
         if out_format == "shapefile":
             gdf = gpd.GeoDataFrame(
-                    {'geometry': multi_polygon},
+                    {'geometry': multi_polygon.geoms},
                     crs=self._calc_crs
                     )
             if not crs.equals(self._calc_crs):
@@ -535,7 +535,7 @@ def _write_to_file(
 
         elif out_format == "feather":
             gdf = gpd.GeoDataFrame(
-                    {'geometry': multi_polygon},
+                    {'geometry': multi_polygon.geoms},
                     crs=self._calc_crs
                     )
             if not crs.equals(self._calc_crs):

ocsmesh/raster.py

@@ -578,9 +578,14 @@ def get_multipolygon(
 
         for win in iter_windows:
             x, y, z = self.get_window_data(win, band=band)
-            new_mask = np.full(z.mask.shape, 0)
-            new_mask[np.where(z.mask)] = -1
-            new_mask[np.where(~z.mask)] = 1
+            if z.mask.ndim == 2:
+                new_mask = np.full(z.mask.shape, 0)
+                new_mask[np.where(z.mask)] = -1
+                new_mask[np.where(~z.mask)] = 1
+            else:
+                # If not mask available
+                # NOTE: We want dtype to be int64, not float64
+                new_mask = np.full((len(y), len(x)), 1)
 
             if zmin is not None:
                 new_mask[np.where(z < zmin)] = -1
@@ -593,9 +598,10 @@ def get_multipolygon(
 
             fig, ax = plt.subplots()
             ax.contourf(x, y, new_mask, levels=[0, 1])
+            mpoly = utils.get_multipolygon_from_pathplot(ax)
+            if mpoly is not None:
+                polygon_collection.extend(mpoly.geoms)
             plt.close(fig)
-            polygon_collection.extend(
-                    utils.get_multipolygon_from_pathplot(ax).geoms)
 
         union_result = ops.unary_union(polygon_collection)
         if not isinstance(union_result, MultiPolygon):

ocsmesh/utils.py

@@ -184,8 +184,9 @@ def get_mesh_polygons(mesh):
 
             res_gdf = polys_gdf[polys_gdf.intersects(pnt)]
             if len(res_gdf) == 0:
-                # How is this possible?!
-                pnts = MultiPoint([*(pnts.geoms[:idx]), *(pnts.geoms[idx + 1:])])
+                # How is this possible for the remaining points not to
+                # intersect with the remaining polys?!
+                pnts = MultiPoint([pt for i, pt in enumerate(pnts.geoms) if i != idx])
                 if pnts.is_empty:
                     break
 
@@ -625,6 +626,7 @@ def inner_ring_collection(mesh):
 def get_multipolygon_from_pathplot(ax):
     # extract linear_rings from plot
     linear_ring_collection = []
+    multipolygon = None
     for path_collection in ax.collections:
         for path in path_collection.get_paths():
             polygons = path.to_polygons(closed_only=True)
@@ -651,7 +653,7 @@ def get_multipolygon_from_pathplot(ax):
             polygon_collection.append(Polygon(outer_ring, inner_rings))
 
         multipolygon = MultiPolygon(polygon_collection)
-    else:
+    elif len(linear_ring_collection) != 0:
         multipolygon = MultiPolygon(
             [Polygon(linear_ring_collection.pop())])
     return multipolygon

pyproject.toml

@@ -25,7 +25,7 @@ dependencies = [
     "jigsawpy", "matplotlib", "netCDF4", "numba",
     "numpy>=1.21", # introduce npt.NDArray
     "pyarrow", "rtree", "pyproj>=3.0", "rasterio", "requests", "scipy",
-    "shapely>=1.8, <2", "tqdm", "typing_extensions", "utm",
+    "shapely", "tqdm", "typing_extensions", "utm",
     ]
 dynamic = ["version"]
 

tests/api/common.py

@@ -3,6 +3,7 @@
 from pyproj import CRS
 from jigsawpy import jigsaw_msh_t
 
+import ocsmesh
 
 def create_rectangle_mesh(nx, ny, holes, x_extent=None, y_extent=None):
     """
@@ -121,3 +122,37 @@ def msht_from_numpy(
     return mesh
 
 
+def topo_2rast_1mesh(r1_path, r2_path, m1_path):
+
+        rast_xy_1 = np.mgrid[-1:0.1:0.1, -0.7:0.1:0.1]
+        rast_xy_2 = np.mgrid[0:1.1:0.1, -0.7:0.1:0.1]
+        nx_z_1, ny_z_1 = rast_xy_1[0].shape
+        rast_z_1 = np.ones_like(rast_xy_1[0]) * 100
+        rast_z_1[:, ny_z_1*1//3:ny_z_1*2//3] = 10
+        rast_z_1[:, ny_z_1*2//3:] = 0
+        rast_z_1[nx_z_1*7//16:nx_z_1*11//16, :] = -10
+        nx_z_2, ny_z_2 = rast_xy_2[0].shape
+        rast_z_2 = np.ones_like(rast_xy_2[0]) * -100
+        rast_z_2[:, :ny_z_2*1//3] = 0
+        rast_z_2[:, ny_z_2*1//3:ny_z_2*2//3] = -10
+        rast_z_2[nx_z_2*7//16:nx_z_2*11//16, :ny_z_2*2//3] = -10
+
+        raster_from_numpy(
+            r1_path, rast_z_1, rast_xy_1, 4326
+        )
+        raster_from_numpy(
+            r2_path, rast_z_2, rast_xy_2, 4326
+        )
+
+        msh_t = create_rectangle_mesh(
+            nx=17, ny=7, holes=[40, 41], x_extent=(-1, 1), y_extent=(0, 1))
+        msh_t.crs = CRS.from_epsg(4326)
+        msh_t.value[:] = 10
+
+        mesh = ocsmesh.Mesh(msh_t)
+        # NOTE: Assuming the interpolation works fine!
+        mesh.interpolate(
+            [ocsmesh.Raster(i) for i in [r1_path, r2_path]],
+            method='nearest'
+        )
+        mesh.write(str(m1_path), format='grd', overwrite=False)