From 9ae6a6d76f72696e01a9e72f5267c0e93ba6166b Mon Sep 17 00:00:00 2001 From: "Soroosh.Mani" Date: Wed, 13 Sep 2023 17:53:14 -0400 Subject: [PATCH 01/26] Use triangle to improve subsetting transition --- ocsmesh/cli/subset_n_combine.py | 168 +++++++++++++++++--------------- ocsmesh/geom/base.py | 42 +------- ocsmesh/utils.py | 132 +++++++++++++++++++++++++ pyproject.toml | 2 +- 4 files changed, 226 insertions(+), 118 deletions(-) diff --git a/ocsmesh/cli/subset_n_combine.py b/ocsmesh/cli/subset_n_combine.py index be863e56..69f29faa 100755 --- a/ocsmesh/cli/subset_n_combine.py +++ b/ocsmesh/cli/subset_n_combine.py @@ -11,11 +11,12 @@ from jigsawpy.jig_t import jigsaw_jig_t import numpy as np from pyproj import CRS, Transformer +from scipy.interpolate import griddata from scipy.spatial import cKDTree from shapely.geometry import MultiPolygon, Polygon, GeometryCollection from shapely.ops import polygonize, unary_union -from ocsmesh import Raster, Geom, Mesh, Hfun, utils +from ocsmesh import Raster, Geom, Mesh, Hfun, utils, JigsawDriver logging.basicConfig( stream=sys.stdout, @@ -242,96 +243,101 @@ def _add_overlap_to_polygon(self, mesh, polygon): return new_polygon, clipped_mesh - def _calculate_mesh_size_function(self, hires_mesh_clip, lores_mesh_clip, crs): + def _calculate_mesh_size_function( + self, + buffer_domain, + hires_mesh_clip, + lores_mesh_clip, + buffer_crs + ): - # calculate mesh size for clipped bits - hfun_hires = Hfun(Mesh(deepcopy(hires_mesh_clip))) - hfun_hires.size_from_mesh() - hfun_lowres = Hfun(Mesh(deepcopy(lores_mesh_clip))) - hfun_lowres.size_from_mesh() + assert(buffer_crs == hires_mesh_clip.crs == lores_mesh_clip.crs) -# _logger.info("Writing hfun clips...") -# start = time() -# hfun_hires.mesh.write(str(out_dir / "hfun_fine.2dm"), format="2dm", overwrite=True) -# hfun_lowres.mesh.write(str(out_dir / "hfun_coarse.2dm"), format="2dm", overwrite=True) -# _logger.info(f"Done in {time() - start} sec") + # HARDCODED FOR NOW + approx_elem_per_width = 3 - jig_hfun = utils.merge_msh_t( - hfun_hires.msh_t(), hfun_lowres.msh_t(), - out_crs=crs)#jig_geom) ### TODO: CRS MUST BE == GEOM_MSHT CRS + utm = utils.estimate_bounds_utm( + buffer_domain.bounds, buffer_crs + ) - return jig_hfun + # calculate mesh size for clipped bits + msht_hi = deepcopy(hires_mesh_clip) + utils.reproject(msht_hi, utm) + hfun_hi = Hfun(Mesh(msht_hi)) + hfun_hi.size_from_mesh() + + msht_lo = deepcopy(lores_mesh_clip) + utils.reproject(msht_lo, utm) + hfun_lo = Hfun(Mesh(msht_lo)) + hfun_lo.size_from_mesh() + + msht_cdt = utils.triangulate_polygon( + buffer_domain, None, opts='p' + ) + msht_cdt.crs = buffer_crs + utils.reproject(msht_cdt, utm) + hfun_cdt = Hfun(Mesh(msht_cdt)) + hfun_cdt.size_from_mesh() - def _generate_mesh_for_buffer_region( - self, buffer_polygon, jig_hfun, crs): + hfun_cdt_sz = deepcopy(hfun_cdt.msh_t().value) / approx_elem_per_width + msht_cdt.value[:] = hfun_cdt_sz + hfun_rep = Hfun(Mesh(msht_cdt)) - seam = Geom(buffer_polygon, crs=crs) + mesh_domain_rep = JigsawDriver( + geom=Geom(buffer_domain, crs=buffer_crs), + hfun=hfun_rep, + initial_mesh=False + ).run() - jig_geom = seam.msh_t() + msht_domain_rep = deepcopy(mesh_domain_rep.msh_t) + utils.reproject(msht_domain_rep, utm) - # IMPORTANT: Setting these to -1 results in NON CONFORMAL boundary -# jig_geom.vert2['IDtag'][:] = -1 -# jig_geom.edge2['IDtag'][:] = -1 + pts_2mesh = np.vstack( + (hfun_hi.msh_t().vert2['coord'], hfun_lo.msh_t().vert2['coord']) + ) + val_2mesh = np.vstack( + (hfun_hi.msh_t().value, hfun_lo.msh_t().value) + ) + domain_sz_1 = griddata( + pts_2mesh, val_2mesh, msht_domain_rep.vert2['coord'], method='linear' + ) + domain_sz_2 = griddata( + pts_2mesh, val_2mesh, msht_domain_rep.vert2['coord'], method='nearest' + ) + domain_sz = domain_sz_1.copy() + domain_sz[np.isnan(domain_sz_1)] = domain_sz_2[np.isnan(domain_sz_1)] - jig_init = deepcopy(seam.msh_t()) - jig_init.vert2['IDtag'][:] = -1 - jig_init.edge2['IDtag'][:] = -1 + msht_domain_rep.value[:] = domain_sz + hfun_interp = Hfun(Mesh(msht_domain_rep)) - # Calculate length of all edges on geom - geom_edges = jig_geom.edge2['index'] - geom_coords = jig_geom.vert2['coord'][geom_edges, :] - geom_edg_lens = np.sqrt(np.sum( - np.power(np.abs(np.diff(geom_coords, axis=1)), 2), - axis=2)).squeeze() + return hfun_interp - # TODO: Use marche to calculate mesh size in the area between - # the two regions? - _logger.info("Meshing...") - start = time() - opts = jigsaw_jig_t() - opts.hfun_scal = "absolute" - opts.hfun_hmin = np.min(geom_edg_lens) - opts.hfun_hmax = np.max(geom_edg_lens) -# opts.hfun_hmin = np.min(jig_hfun.value.ravel()) -# opts.hfun_hmax = np.max(jig_hfun.value.ravel()) - opts.optm_zip_ = False - opts.optm_div_ = False - opts.mesh_dims = +2 - opts.mesh_rad2 = 1.05 -# opts.mesh_rad2 = 2.0 - - jig_mesh = jigsaw_msh_t() - jig_mesh.mshID = 'euclidean-mesh' - jig_mesh.ndims = 2 - jig_mesh.crs = jig_geom.crs - - jigsawpy.lib.jigsaw( - opts, jig_geom, jig_mesh, -# hfun=jig_hfun, - init=jig_init - ) - - jig_mesh.value = np.zeros((jig_mesh.vert2.shape[0], 1)) - self._transform_mesh(jig_mesh, crs) - - # NOTE: Remove out of domain elements (some corner cases!) - elems = jig_mesh.tria3['index'] - coord = jig_mesh.vert2['coord'] - - gdf_elems = gpd.GeoDataFrame( - geometry=[Polygon(tri) for tri in coord[elems]], - crs=jig_mesh.crs - ) - idx = gdf_elems.representative_point().sindex.query( - seam.get_multipolygon(), predicate='intersects' + def _generate_mesh_for_buffer_region( + self, buffer_polygon, hfun_buffer, buffer_crs): + + mesh_domain_interp = JigsawDriver( + geom=Geom(buffer_polygon, crs=buffer_crs), + hfun=hfun_buffer, + initial_mesh=False + ).run() + + msht_domain_interp = deepcopy(mesh_domain_interp.msh_t) + in_verts_mask = np.ones_like(msht_domain_interp.vert2, dtype=bool) + in_verts_mask[ + np.unique(utils.get_boundary_edges(msht_domain_interp)) + ] = False + + utils.reproject(msht_domain_interp, buffer_crs) + msht_buffer = utils.triangulate_polygon( + shape=buffer_polygon, + aux_pts=msht_domain_interp.vert2[in_verts_mask]['coord'], + opts='p' ) - flag = np.zeros(len(gdf_elems), dtype=bool) - flag[idx] = True - jig_mesh.tria3 = jig_mesh.tria3[flag] + msht_buffer.crs = buffer_crs - return jig_mesh + return msht_buffer def _transform_mesh(self, mesh, out_crs): @@ -623,9 +629,13 @@ def _main( poly_seam = poly_seam_8 - # TODO: Get CRS correctly from geom utm -# jig_hfun = self._calculate_mesh_size_function(jig_clip_hires, jig_clip_lowres, crs) - jig_buffer_mesh = self._generate_mesh_for_buffer_region(poly_seam, None, crs) + hfun_buffer = self._calculate_mesh_size_function( + poly_seam, jig_clip_hires, jig_clip_lowres, crs + ) + jig_buffer_mesh = self._generate_mesh_for_buffer_region( + poly_seam, hfun_buffer, crs + ) + # TODO: UTM TO CRS? _logger.info("Combining meshes...") start = time() diff --git a/ocsmesh/geom/base.py b/ocsmesh/geom/base.py index c38f1625..05091cea 100644 --- a/ocsmesh/geom/base.py +++ b/ocsmesh/geom/base.py @@ -144,42 +144,8 @@ def multipolygon_to_jigsaw_msh_t( transformer = Transformer.from_crs(crs, utm_crs, always_xy=True) multipolygon = ops.transform(transformer.transform, multipolygon) - vert2: List[Tuple[Tuple[float, float], int]] = [] - for polygon in multipolygon.geoms: - if np.all( - np.asarray( - polygon.exterior.coords).flatten() == float('inf')): - raise NotImplementedError("ellispoidal-mesh") - for x, y in polygon.exterior.coords[:-1]: - vert2.append(((x, y), 0)) - for interior in polygon.interiors: - for x, y in interior.coords[:-1]: - vert2.append(((x, y), 0)) - - # edge2 - edge2: List[Tuple[int, int]] = [] - for polygon in multipolygon.geoms: - polygon = [polygon.exterior, *polygon.interiors] - for linear_ring in polygon: - _edge2 = [] - for i in range(len(linear_ring.coords)-2): - _edge2.append((i, i+1)) - _edge2.append((_edge2[-1][1], _edge2[0][0])) - edge2.extend( - [(e0+len(edge2), e1+len(edge2)) - for e0, e1 in _edge2]) - # geom - geom = jigsaw_msh_t() - geom.ndims = +2 - geom.mshID = 'euclidean-mesh' - # TODO: Consider ellipsoidal case. - # geom.mshID = 'euclidean-mesh' if self._ellipsoid is None \ - # else 'ellipsoidal-mesh' - geom.vert2 = np.asarray(vert2, dtype=jigsaw_msh_t.VERT2_t) - geom.edge2 = np.asarray( - [((e0, e1), 0) for e0, e1 in edge2], - dtype=jigsaw_msh_t.EDGE2_t) - geom.crs = crs + msht = utils.shape_to_msh_t(multipolygon) + msht.crs = crs if utm_crs is not None: - geom.crs = utm_crs - return geom + msht.crs = utm_crs + return msht diff --git a/ocsmesh/utils.py b/ocsmesh/utils.py index 685fe14a..80983e1f 100644 --- a/ocsmesh/utils.py +++ b/ocsmesh/utils.py @@ -13,10 +13,12 @@ import numpy as np # type: ignore[import] import numpy.typing as npt import rasterio as rio +import triangle as tr from pyproj import CRS, Transformer # type: ignore[import] from scipy.interpolate import ( # type: ignore[import] RectBivariateSpline, griddata) from scipy import sparse, constants +from shapely import force_2d from shapely.geometry import ( # type: ignore[import] Polygon, MultiPolygon, box, GeometryCollection, Point, MultiPoint, @@ -2190,3 +2192,133 @@ def msht_from_numpy( ) return mesh + + +def shape_to_msh_t(shape: Union[Polygon, MultiPolygon]) -> jigsaw_msh_t: + """Calculate vertex-edge representation of polygon shape + + Calculate `jigsawpy` vertex-edge representation of the input + `shapely` shape. + + Parameters + ---------- + shape : Polygon or MultiPolygon + Input polygon for which the vertex-edge representation is to + be calculated + + Returns + ------- + jigsaw_msh_t + Vertex-edge representation of the input shape + + Raises + ------ + NotImplementedError + """ + + vert2: List[Tuple[Tuple[float, float], int]] = [] + if isinstance(shape, Polygon): + shape = MultiPolygon([shape]) + for polygon in shape.geoms: + if np.all( + np.asarray( + polygon.exterior.coords).flatten() == float('inf')): + raise NotImplementedError("ellispoidal-mesh") + for x, y in polygon.exterior.coords[:-1]: + vert2.append(((x, y), 0)) + for interior in polygon.interiors: + for x, y in interior.coords[:-1]: + vert2.append(((x, y), 0)) + + # edge2 + edge2: List[Tuple[int, int]] = [] + for polygon in shape.geoms: + polygon = [polygon.exterior, *polygon.interiors] + for linear_ring in polygon: + edg = [] + for i in range(len(linear_ring.coords) - 2): + edg.append((i, i+1)) + edg.append((edg[-1][1], edg[0][0])) + edge2.extend( + [(e0+len(edge2), e1+len(edge2)) + for e0, e1 in edg]) + + msht = jigsaw_msh_t() + msht.ndims = +2 + msht.mshID = 'euclidean-mesh' + msht.vert2 = np.asarray(vert2, dtype=jigsaw_msh_t.VERT2_t) + msht.edge2 = np.asarray( + [((e0, e1), 0) for e0, e1 in edge2], + dtype=jigsaw_msh_t.EDGE2_t) + return msht + + +def triangulate_polygon( + shape: Union[Polygon, MultiPolygon, gpd.GeoDataFrame, gpd.GeoSeries], + aux_pts: Union[np.array, Point, MultiPoint, gpd.GeoDataFrame, gpd.GeoSeries] = None, + opts='p', +) -> None: + ''' + + ''' + + # NOTE: Triangle can handle separate closed polygons, + # so no need to have for loop + + if isinstance(shape, (gpd.GeoDataFrame, gpd.GeoSeries)): + shape = shape.unary_union + + if isinstance(shape, Polygon): + shape = MultiPolygon([shape]) + + if not isinstance(shape, (MultiPolygon, Polygon)): + raise ValueError("Input shape must be convertible to polygon!") + + msht_shape = shape_to_msh_t(shape) + coords = msht_shape.vert2['coord'] + edges = msht_shape.edge2['index'] + + holes = [] + for poly in shape.geoms: + for ring in poly.interiors: + holes.append( + np.array(Polygon(ring).representative_point().xy).ravel() + ) + holes = np.array(holes) + + + if aux_pts is not None: + if isinstance(aux_pts, (gpd.GeoDataFrame, gpd.GeoSeries)): + aux_pts = np.array(aux_pts.unary_union.xy).T + + elif isinstance(aux_pts, (Point, MultiPolygon)): + aux_pts = np.array(aux_pts.xy).T + + coords = np.vstack((coords, aux_pts)) + + input_dict = {'vertices': coords, 'segments': edges} + if len(holes): + input_dict['holes'] = holes + cdt = tr.triangulate(input_dict, opts=opts) + + cells = cdt['triangles'] + points = cdt['vertices'] + + msht_tri = jigsaw_msh_t() + msht_tri.mshID = 'euclidean-mesh' + msht_tri.ndims = 2 + + msht_tri.vert2 = np.array( + [(crd, -1) for crd in points], + dtype=jigsaw_msh_t.VERT2_t + ) + msht_tri.tria3 = np.array( + [(cell, -1) for cell in cells], + dtype=jigsaw_msh_t.TRIA3_t + ) + msht_tri.value = np.array( + np.zeros((len(points) ,1)), + dtype=jigsaw_msh_t.REALS_t + ) + + return msht_tri diff --git a/pyproject.toml b/pyproject.toml index bcaafe11..6e60e156 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -25,7 +25,7 @@ dependencies = [ "jigsawpy", "matplotlib", "netCDF4", "numba", "numpy>=1.21", # introduce npt.NDArray "pyarrow", "rtree", "pyproj>=3.0", "rasterio", "scipy", - "shapely", "typing_extensions", "utm", + "shapely", "triangle", "typing_extensions", "utm", ] dynamic = ["version"] From 19506ab13b0ea1471558ebe3781a062013a62e0d Mon Sep 17 00:00:00 2001 From: "Soroosh.Mani" Date: Wed, 13 Sep 2023 20:36:05 -0400 Subject: [PATCH 02/26] Use shapely2 funcs & add new shp2msht func --- ocsmesh/cli/subset_n_combine.py | 22 +++++--- ocsmesh/utils.py | 94 ++++++++++++++++++++++++++++++++- pyproject.toml | 2 +- 3 files changed, 109 insertions(+), 9 deletions(-) diff --git a/ocsmesh/cli/subset_n_combine.py b/ocsmesh/cli/subset_n_combine.py index 69f29faa..294c1e79 100755 --- a/ocsmesh/cli/subset_n_combine.py +++ b/ocsmesh/cli/subset_n_combine.py @@ -14,7 +14,7 @@ from scipy.interpolate import griddata from scipy.spatial import cKDTree from shapely.geometry import MultiPolygon, Polygon, GeometryCollection -from shapely.ops import polygonize, unary_union +from shapely.ops import polygonize, unary_union, transform from ocsmesh import Raster, Geom, Mesh, Hfun, utils, JigsawDriver @@ -259,6 +259,8 @@ def _calculate_mesh_size_function( utm = utils.estimate_bounds_utm( buffer_domain.bounds, buffer_crs ) + transformer = Transformer.from_crs(buffer_crs, utm, always_xy=True) + buffer_domain = transform(transformer.transform, buffer_domain) # calculate mesh size for clipped bits msht_hi = deepcopy(hires_mesh_clip) @@ -274,9 +276,9 @@ def _calculate_mesh_size_function( msht_cdt = utils.triangulate_polygon( buffer_domain, None, opts='p' ) - msht_cdt.crs = buffer_crs + msht_cdt.crs = utm - utils.reproject(msht_cdt, utm) +# utils.reproject(msht_cdt, utm) hfun_cdt = Hfun(Mesh(msht_cdt)) hfun_cdt.size_from_mesh() @@ -285,7 +287,7 @@ def _calculate_mesh_size_function( hfun_rep = Hfun(Mesh(msht_cdt)) mesh_domain_rep = JigsawDriver( - geom=Geom(buffer_domain, crs=buffer_crs), + geom=Geom(buffer_domain, crs=utm), hfun=hfun_rep, initial_mesh=False ).run() @@ -317,8 +319,14 @@ def _calculate_mesh_size_function( def _generate_mesh_for_buffer_region( self, buffer_polygon, hfun_buffer, buffer_crs): + utm = utils.estimate_bounds_utm( + buffer_polygon.bounds, buffer_crs + ) + transformer = Transformer.from_crs(buffer_crs, utm, always_xy=True) + buffer_polygon = transform(transformer.transform, buffer_polygon) + mesh_domain_interp = JigsawDriver( - geom=Geom(buffer_polygon, crs=buffer_crs), + geom=Geom(buffer_polygon, crs=utm), hfun=hfun_buffer, initial_mesh=False ).run() @@ -329,13 +337,13 @@ def _generate_mesh_for_buffer_region( np.unique(utils.get_boundary_edges(msht_domain_interp)) ] = False - utils.reproject(msht_domain_interp, buffer_crs) +# utils.reproject(msht_domain_interp, buffer_crs) msht_buffer = utils.triangulate_polygon( shape=buffer_polygon, aux_pts=msht_domain_interp.vert2[in_verts_mask]['coord'], opts='p' ) - msht_buffer.crs = buffer_crs + msht_buffer.crs = utm return msht_buffer diff --git a/ocsmesh/utils.py b/ocsmesh/utils.py index 80983e1f..2dac519a 100644 --- a/ocsmesh/utils.py +++ b/ocsmesh/utils.py @@ -25,6 +25,7 @@ LineString, LinearRing) from shapely.ops import polygonize, linemerge, unary_union import geopandas as gpd +import pandas as pd import utm @@ -2253,6 +2254,96 @@ def shape_to_msh_t(shape: Union[Polygon, MultiPolygon]) -> jigsaw_msh_t: return msht +def shape_to_msh_t_2(shape: Union[Polygon, MultiPolygon]) -> jigsaw_msh_t: + gdf_shape = shape + if not isinstance(shape, gpd.GeoDataFrame): + gdf_shape = gpd.GeoDataFrame(geometry=[shape]) + + df_pts = ( + gdf_shape + .geometry + .transform(force_2d) + .boundary + .explode(ignore_index=True) + .map(lambda i: i.coords) + .explode() + ) + + + df_pts_2 = ( + pd.DataFrame( + df_pts.tolist(), + index=df_pts.index, + columns=['lon', 'lat'] + ) + .reset_index() + .drop_duplicates() # Polygon index is used for duplicate removal + ) + df_seg = ( + df_pts_2.join( + df_pts_2.groupby('index').transform(np.roll, 1, axis=0), + lsuffix='_1', + rsuffix='_2' + ).dropna() + .set_index('index') + ) + + df_nodes = ( + df_pts_2.drop(columns='index') + .drop_duplicates() + .reset_index(drop=True) # renumber nodes + .reset_index() # add node idx as df data column + .set_index(['lon','lat']) + ) + + # CRD Table + df_coo = df_nodes.reset_index().drop(columns='index') + + # CNN Table + df_edg = ( + df_nodes.loc[ + pd.MultiIndex.from_frame(df_seg[['lon_1', 'lat_1']]) + ].reset_index(drop=True) + .join( + df_nodes.loc[ + pd.MultiIndex.from_frame(df_seg[['lon_2', 'lat_2']]) + ] + .reset_index(drop=True), + lsuffix='_1', + rsuffix='_2' + ) + ) + + + ar_edg = np.sort(df_edg.values, axis=1) + df_cnn = ( + pd.DataFrame.from_records( + ar_edg, + columns=['index_1', 'index_2'], + ) + .drop_duplicates() # Remove duplicate edges + .reset_index(drop=True) + ) + + msht = jigsaw_msh_t() + msht.ndims = +2 + msht.mshID = 'euclidean-mesh' + msht.vert2 = np.asarray( + [(crd, 0) for crd in df_coo[['lon', 'lat']].values], + dtype=jigsaw_msh_t.VERT2_t + ) + msht.edge2 = np.asarray( + [((e0, e1), 0) for e0, e1 in df_cnn.values], + dtype=jigsaw_msh_t.EDGE2_t + ) + msht.value = np.array( + np.zeros((len(msht.vert2) ,1)), + dtype=jigsaw_msh_t.REALS_t + ) + return msht + + + def triangulate_polygon( shape: Union[Polygon, MultiPolygon, gpd.GeoDataFrame, gpd.GeoSeries], aux_pts: Union[np.array, Point, MultiPoint, gpd.GeoDataFrame, gpd.GeoSeries] = None, @@ -2274,7 +2365,8 @@ def triangulate_polygon( if not isinstance(shape, (MultiPolygon, Polygon)): raise ValueError("Input shape must be convertible to polygon!") - msht_shape = shape_to_msh_t(shape) +# msht_shape = shape_to_msh_t(shape) + msht_shape = shape_to_msh_t_2(shape) coords = msht_shape.vert2['coord'] edges = msht_shape.edge2['index'] diff --git a/pyproject.toml b/pyproject.toml index 6e60e156..889aced1 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -25,7 +25,7 @@ dependencies = [ "jigsawpy", "matplotlib", "netCDF4", "numba", "numpy>=1.21", # introduce npt.NDArray "pyarrow", "rtree", "pyproj>=3.0", "rasterio", "scipy", - "shapely", "triangle", "typing_extensions", "utm", + "shapely>=2", "triangle", "typing_extensions", "utm", ] dynamic = ["version"] From b0771b47172406cda33bc4a595f5a45a159886d2 Mon Sep 17 00:00:00 2001 From: "Soroosh.Mani" Date: Wed, 13 Sep 2023 20:48:07 -0400 Subject: [PATCH 03/26] Fix CRS issue in new merge --- ocsmesh/cli/subset_n_combine.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/ocsmesh/cli/subset_n_combine.py b/ocsmesh/cli/subset_n_combine.py index 294c1e79..be51d3f6 100755 --- a/ocsmesh/cli/subset_n_combine.py +++ b/ocsmesh/cli/subset_n_combine.py @@ -345,6 +345,8 @@ def _generate_mesh_for_buffer_region( ) msht_buffer.crs = utm + utils.reproject(msht_buffer, buffer_crs) + return msht_buffer From 91fb3c945be404d06d0313f7bbf8c266dc9df349 Mon Sep 17 00:00:00 2001 From: "Soroosh.Mani" Date: Thu, 14 Sep 2023 09:18:57 -0400 Subject: [PATCH 04/26] Fix issue related to separate multipolygon subset buffer regions --- ocsmesh/cli/subset_n_combine.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/ocsmesh/cli/subset_n_combine.py b/ocsmesh/cli/subset_n_combine.py index be51d3f6..94d7a220 100755 --- a/ocsmesh/cli/subset_n_combine.py +++ b/ocsmesh/cli/subset_n_combine.py @@ -290,7 +290,7 @@ def _calculate_mesh_size_function( geom=Geom(buffer_domain, crs=utm), hfun=hfun_rep, initial_mesh=False - ).run() + ).run(sieve=0) msht_domain_rep = deepcopy(mesh_domain_rep.msh_t) utils.reproject(msht_domain_rep, utm) @@ -329,7 +329,7 @@ def _generate_mesh_for_buffer_region( geom=Geom(buffer_polygon, crs=utm), hfun=hfun_buffer, initial_mesh=False - ).run() + ).run(sieve=0) msht_domain_interp = deepcopy(mesh_domain_interp.msh_t) in_verts_mask = np.ones_like(msht_domain_interp.vert2, dtype=bool) From 138f6c11345bbe4e010785f37020c2aa5a0225eb Mon Sep 17 00:00:00 2001 From: "Soroosh.Mani" Date: Thu, 14 Sep 2023 10:06:10 -0400 Subject: [PATCH 05/26] Cleanup and remove shapely2 func --- ocsmesh/utils.py | 25 ++++++++----------------- pyproject.toml | 2 +- 2 files changed, 9 insertions(+), 18 deletions(-) diff --git a/ocsmesh/utils.py b/ocsmesh/utils.py index 2dac519a..e34090f5 100644 --- a/ocsmesh/utils.py +++ b/ocsmesh/utils.py @@ -18,7 +18,6 @@ from scipy.interpolate import ( # type: ignore[import] RectBivariateSpline, griddata) from scipy import sparse, constants -from shapely import force_2d from shapely.geometry import ( # type: ignore[import] Polygon, MultiPolygon, box, GeometryCollection, Point, MultiPoint, @@ -2259,29 +2258,21 @@ def shape_to_msh_t_2(shape: Union[Polygon, MultiPolygon]) -> jigsaw_msh_t: if not isinstance(shape, gpd.GeoDataFrame): gdf_shape = gpd.GeoDataFrame(geometry=[shape]) - df_pts = ( + df_lonlat = ( gdf_shape .geometry - .transform(force_2d) .boundary .explode(ignore_index=True) .map(lambda i: i.coords) .explode() + .apply(lambda s: pd.Series({'lon': s[0], 'lat': s[1]})) + .reset_index() # put polygon index in the dataframe + .drop_duplicates() # drop duplicates within polygons ) - - df_pts_2 = ( - pd.DataFrame( - df_pts.tolist(), - index=df_pts.index, - columns=['lon', 'lat'] - ) - .reset_index() - .drop_duplicates() # Polygon index is used for duplicate removal - ) df_seg = ( - df_pts_2.join( - df_pts_2.groupby('index').transform(np.roll, 1, axis=0), + df_lonlat.join( + df_lonlat.groupby('index').transform(np.roll, 1, axis=0), lsuffix='_1', rsuffix='_2' ).dropna() @@ -2289,8 +2280,8 @@ def shape_to_msh_t_2(shape: Union[Polygon, MultiPolygon]) -> jigsaw_msh_t: ) df_nodes = ( - df_pts_2.drop(columns='index') - .drop_duplicates() + df_lonlat.drop(columns='index') # drop to allow cross poly dupl + .drop_duplicates() # drop duplicates across multiple polygons .reset_index(drop=True) # renumber nodes .reset_index() # add node idx as df data column .set_index(['lon','lat']) diff --git a/pyproject.toml b/pyproject.toml index 889aced1..6e60e156 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -25,7 +25,7 @@ dependencies = [ "jigsawpy", "matplotlib", "netCDF4", "numba", "numpy>=1.21", # introduce npt.NDArray "pyarrow", "rtree", "pyproj>=3.0", "rasterio", "scipy", - "shapely>=2", "triangle", "typing_extensions", "utm", + "shapely", "triangle", "typing_extensions", "utm", ] dynamic = ["version"] From 4e674e740d098a2a3b8ccd2bf80160da2ca72f95 Mon Sep 17 00:00:00 2001 From: "Soroosh.Mani" Date: Thu, 14 Sep 2023 12:04:26 -0400 Subject: [PATCH 06/26] Code reorg & cleanup --- ocsmesh/utils.py | 40 ++++++++++++++++++++++++---------------- 1 file changed, 24 insertions(+), 16 deletions(-) diff --git a/ocsmesh/utils.py b/ocsmesh/utils.py index e34090f5..1175d580 100644 --- a/ocsmesh/utils.py +++ b/ocsmesh/utils.py @@ -2168,6 +2168,7 @@ def msht_from_numpy( coordinates, triangles, quadrilaterals=None, + values=None, crs=CRS.from_epsg(4326) ): mesh = jigsaw_msh_t() @@ -2190,6 +2191,18 @@ def msht_from_numpy( [(index, 0) for index in quadrilaterals], dtype=jigsaw_msh_t.QUAD4_t ) + if values is None: + values = np.array( + np.zeros((len(mesh.vert2), 1)), + dtype=jigsaw_msh_t.REALS_t + ) + elif values.shape != (len(mesh.vert2), 1): + raise ValueError( + "Input for mesh values must either be None or a" + " 2D-array with row count equal to vertex count!" + ) + + mesh.value = values return mesh @@ -2253,7 +2266,7 @@ def shape_to_msh_t(shape: Union[Polygon, MultiPolygon]) -> jigsaw_msh_t: return msht -def shape_to_msh_t_2(shape: Union[Polygon, MultiPolygon]) -> jigsaw_msh_t: +def shape_to_msh_t_2(shape: Union[Polygon, MultiPolygon], crs=None) -> jigsaw_msh_t: gdf_shape = shape if not isinstance(shape, gpd.GeoDataFrame): gdf_shape = gpd.GeoDataFrame(geometry=[shape]) @@ -2316,21 +2329,15 @@ def shape_to_msh_t_2(shape: Union[Polygon, MultiPolygon]) -> jigsaw_msh_t: .reset_index(drop=True) ) - msht = jigsaw_msh_t() - msht.ndims = +2 - msht.mshID = 'euclidean-mesh' - msht.vert2 = np.asarray( - [(crd, 0) for crd in df_coo[['lon', 'lat']].values], - dtype=jigsaw_msh_t.VERT2_t - ) - msht.edge2 = np.asarray( - [((e0, e1), 0) for e0, e1 in df_cnn.values], - dtype=jigsaw_msh_t.EDGE2_t - ) - msht.value = np.array( - np.zeros((len(msht.vert2) ,1)), - dtype=jigsaw_msh_t.REALS_t + + msht = msht_from_numpy( + coordinates=df_coo[['lon', 'lat']].values, + triangles=df_cnn.values, + quadrilaterals=None, + values=None, + crs=crs, ) + return msht @@ -2339,6 +2346,7 @@ def triangulate_polygon( shape: Union[Polygon, MultiPolygon, gpd.GeoDataFrame, gpd.GeoSeries], aux_pts: Union[np.array, Point, MultiPoint, gpd.GeoDataFrame, gpd.GeoSeries] = None, opts='p', + crs=None, ) -> None: ''' @@ -2357,7 +2365,7 @@ def triangulate_polygon( raise ValueError("Input shape must be convertible to polygon!") # msht_shape = shape_to_msh_t(shape) - msht_shape = shape_to_msh_t_2(shape) + msht_shape = shape_to_msh_t_2(shape, crs) coords = msht_shape.vert2['coord'] edges = msht_shape.edge2['index'] From 69b855677c440a26d3c7d2f3d4aa74a59cc72106 Mon Sep 17 00:00:00 2001 From: "Soroosh.Mani" Date: Thu, 14 Sep 2023 14:54:31 -0400 Subject: [PATCH 07/26] Bugfix for cleanup! --- ocsmesh/utils.py | 68 +++++++++++++++++++++--------------------------- 1 file changed, 30 insertions(+), 38 deletions(-) diff --git a/ocsmesh/utils.py b/ocsmesh/utils.py index 1175d580..6ad2f2d9 100644 --- a/ocsmesh/utils.py +++ b/ocsmesh/utils.py @@ -2166,7 +2166,9 @@ def raster_from_numpy( def msht_from_numpy( coordinates, - triangles, + *, # Get everything else as keyword args + edges=None, + triangles=None, quadrilaterals=None, values=None, crs=CRS.from_epsg(4326) @@ -2182,25 +2184,32 @@ def msht_from_numpy( [(coord, 0) for coord in coordinates], dtype=jigsaw_msh_t.VERT2_t ) - mesh.tria3 = np.array( - [(index, 0) for index in triangles], - dtype=jigsaw_msh_t.TRIA3_t + + if edges is not None: + mesh.edge2 = np.array( + [(index, 0) for index in edges], + dtype=jigsaw_msh_t.EDGE2_t + ) + if triangles is not None: + mesh.tria3 = np.array( + [(index, 0) for index in triangles], + dtype=jigsaw_msh_t.TRIA3_t ) if quadrilaterals is not None: mesh.quad4 = np.array( [(index, 0) for index in quadrilaterals], dtype=jigsaw_msh_t.QUAD4_t - ) + ) if values is None: values = np.array( np.zeros((len(mesh.vert2), 1)), dtype=jigsaw_msh_t.REALS_t ) elif values.shape != (len(mesh.vert2), 1): - raise ValueError( - "Input for mesh values must either be None or a" - " 2D-array with row count equal to vertex count!" - ) + raise ValueError( + "Input for mesh values must either be None or a" + " 2D-array with row count equal to vertex count!" + ) mesh.value = values @@ -2266,7 +2275,7 @@ def shape_to_msh_t(shape: Union[Polygon, MultiPolygon]) -> jigsaw_msh_t: return msht -def shape_to_msh_t_2(shape: Union[Polygon, MultiPolygon], crs=None) -> jigsaw_msh_t: +def shape_to_msh_t_2(shape: Union[Polygon, MultiPolygon]) -> jigsaw_msh_t: gdf_shape = shape if not isinstance(shape, gpd.GeoDataFrame): gdf_shape = gpd.GeoDataFrame(geometry=[shape]) @@ -2328,16 +2337,14 @@ def shape_to_msh_t_2(shape: Union[Polygon, MultiPolygon], crs=None) -> jigsaw_ms .drop_duplicates() # Remove duplicate edges .reset_index(drop=True) ) - - + msht = msht_from_numpy( coordinates=df_coo[['lon', 'lat']].values, - triangles=df_cnn.values, - quadrilaterals=None, - values=None, - crs=crs, + edges=df_cnn.values, + values=np.zeros((len(df_coo) ,1)), + crs=None ) - + return msht @@ -2346,7 +2353,6 @@ def triangulate_polygon( shape: Union[Polygon, MultiPolygon, gpd.GeoDataFrame, gpd.GeoSeries], aux_pts: Union[np.array, Point, MultiPoint, gpd.GeoDataFrame, gpd.GeoSeries] = None, opts='p', - crs=None, ) -> None: ''' @@ -2364,8 +2370,7 @@ def triangulate_polygon( if not isinstance(shape, (MultiPolygon, Polygon)): raise ValueError("Input shape must be convertible to polygon!") -# msht_shape = shape_to_msh_t(shape) - msht_shape = shape_to_msh_t_2(shape, crs) + msht_shape = shape_to_msh_t_2(shape) coords = msht_shape.vert2['coord'] edges = msht_shape.edge2['index'] @@ -2392,24 +2397,11 @@ def triangulate_polygon( input_dict['holes'] = holes cdt = tr.triangulate(input_dict, opts=opts) - cells = cdt['triangles'] - points = cdt['vertices'] - - msht_tri = jigsaw_msh_t() - msht_tri.mshID = 'euclidean-mesh' - msht_tri.ndims = 2 - - msht_tri.vert2 = np.array( - [(crd, -1) for crd in points], - dtype=jigsaw_msh_t.VERT2_t - ) - msht_tri.tria3 = np.array( - [(cell, -1) for cell in cells], - dtype=jigsaw_msh_t.TRIA3_t - ) - msht_tri.value = np.array( - np.zeros((len(points) ,1)), - dtype=jigsaw_msh_t.REALS_t + msht_tri = msht_from_numpy( + coordinates=cdt['vertices'], + triangles=cdt['triangles'], + values=np.zeros((len(cdt['vertices']) ,1)), + crs=None, ) return msht_tri From 0d1eb8097bb15152a204d631d2711ea3fefd453c Mon Sep 17 00:00:00 2001 From: "Soroosh.Mani" Date: Thu, 14 Sep 2023 16:10:51 -0400 Subject: [PATCH 08/26] Pin python version to below 3.11 for triangle support --- .github/workflows/documentation.yml | 2 +- .github/workflows/functional_test.yml | 2 +- .github/workflows/functional_test_2.yml | 2 +- .github/workflows/pylint.yml | 2 +- pyproject.toml | 2 +- 5 files changed, 5 insertions(+), 5 deletions(-) diff --git a/.github/workflows/documentation.yml b/.github/workflows/documentation.yml index 98b5d66c..bead88e8 100644 --- a/.github/workflows/documentation.yml +++ b/.github/workflows/documentation.yml @@ -23,7 +23,7 @@ jobs: environment-file: ./environment.yml environment-name: ocsmesh-env extra-specs: | - python=3.* + python=3.10 - name: Install dependencies shell: bash -l {0} run: | diff --git a/.github/workflows/functional_test.yml b/.github/workflows/functional_test.yml index 6a4f753f..728c445e 100644 --- a/.github/workflows/functional_test.yml +++ b/.github/workflows/functional_test.yml @@ -28,7 +28,7 @@ jobs: fail-fast: false matrix: os: [ ubuntu-latest ] - python-version: [ '3.9', '3.10', '3.11' ] + python-version: [ '3.9', '3.10' ] steps: - uses: actions/checkout@v3 diff --git a/.github/workflows/functional_test_2.yml b/.github/workflows/functional_test_2.yml index 92153e6d..fe6a7182 100644 --- a/.github/workflows/functional_test_2.yml +++ b/.github/workflows/functional_test_2.yml @@ -28,7 +28,7 @@ jobs: fail-fast: false matrix: os: [ ubuntu-latest ] - python-version: [ '3.9', '3.10', '3.11' ] + python-version: [ '3.9', '3.10' ] steps: - name: Checkout diff --git a/.github/workflows/pylint.yml b/.github/workflows/pylint.yml index 3a8b0e20..12ee8bb5 100644 --- a/.github/workflows/pylint.yml +++ b/.github/workflows/pylint.yml @@ -31,7 +31,7 @@ jobs: environment-file: ./environment.yml environment-name: ocsmesh-env extra-specs: | - python=3.* + python=3.10 - name: Install dependencies shell: bash -l {0} run: | diff --git a/pyproject.toml b/pyproject.toml index 6e60e156..d0886b8c 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -19,7 +19,7 @@ maintainers = [ description = "Package to generate computational unstructured meshes from planetary modeling." license = {file = "LICENSE"} readme = "README.md" -requires-python = '>=3.9' # 3.8 not supported by scipy +requires-python = '>=3.9, <3.11' # 3.8 -> scipy, 3.11 -> triangle dependencies = [ "colored-traceback", "fiona", "geopandas", "jigsawpy", "matplotlib", "netCDF4", "numba", From 44fbabfc95ad64a2dcd45d2664dbfdc20017e553 Mon Sep 17 00:00:00 2001 From: "Soroosh.Mani" Date: Thu, 14 Sep 2023 17:23:06 -0400 Subject: [PATCH 09/26] Fix tests for the now keyword only args to util func --- tests/api/hfun.py | 4 +++- tests/api/utils.py | 4 +++- 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/tests/api/hfun.py b/tests/api/hfun.py index 3570545b..1f43f9f7 100755 --- a/tests/api/hfun.py +++ b/tests/api/hfun.py @@ -678,7 +678,9 @@ def setUp(self): [0, 2, 6], [5, 4, 7], ]) - msh_t = ocsmesh.utils.msht_from_numpy(crd, tria, crs=4326) + msh_t = ocsmesh.utils.msht_from_numpy( + crd, triangles=tria, crs=4326 + ) mesh = ocsmesh.Mesh(msh_t) mesh.write(str(self.mesh1), format='grd', overwrite=False) diff --git a/tests/api/utils.py b/tests/api/utils.py index d555b066..c343e3e2 100644 --- a/tests/api/utils.py +++ b/tests/api/utils.py @@ -67,7 +67,9 @@ def test_cleanup_mesh_and_generate_valid_mesh(self): quads = msh_t1.quad4['index'] quads = np.vstack((quads, msh_t2.quad4['index'] + len(msh_t1.vert2))) - msh_t = utils.msht_from_numpy(verts, trias, quads) + msh_t = utils.msht_from_numpy( + verts, triangles=trias, quadrilaterals=quads + ) utils.finalize_mesh(msh_t) From 446c0ad4de5da544e7b514a70b5d7b06f30cf551 Mon Sep 17 00:00:00 2001 From: "Soroosh.Mani" Date: Thu, 14 Sep 2023 17:38:16 -0400 Subject: [PATCH 10/26] Add empty placeholder test - RED --- tests/api/utils.py | 53 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 53 insertions(+) diff --git a/tests/api/utils.py b/tests/api/utils.py index c343e3e2..dc10efd2 100644 --- a/tests/api/utils.py +++ b/tests/api/utils.py @@ -590,6 +590,16 @@ def test_crs(self): self.assertEqual(out_msht_2.crs, CRS.from_user_input('esri:102008')) + def test_values(self): + # Test it has values + # Test value input of wrong size + # Test value None input + self.assertTrue(False) + + + def test_kwonly_args(self): + self.assertTrue(False) + class CreateRasterFromNumpy(unittest.TestCase): @@ -655,5 +665,48 @@ def test_data_masking(self): self.assertEqual(rast.src.nodata, fill_value) +class ShapeToMeshT(unittest.TestCase): + + def test_input_output_types(self): + self.assertTrue(False) + + + def test_old_implementation(self): + self.assertTrue(False) + + + def test_new_implementation(self): + self.assertTrue(False) + + + + +class TestTriangulatePolygon(unittest.TestCase): + + + def test_different_input_types(self): + self.assertTrue(False) + + + def test_wrong_input_of_right_type(self): + self.assertTrue(False) + + + def test_output_types(self): + self.assertTrue(False) + + + def test_fixed_boundary(self): + self.assertTrue(False) + + + def test_aux_points(self) + self.assertTrue(False) + + + def test_polygon_holes(self) + self.assertTrue(False) + + if __name__ == '__main__': unittest.main() From 7904b61008f5108782ca2543f9be35cd405e4541 Mon Sep 17 00:00:00 2001 From: "Soroosh.Mani" Date: Fri, 15 Sep 2023 15:28:50 -0400 Subject: [PATCH 11/26] Add tests and more sanity checks to code --- ocsmesh/utils.py | 24 +++++- tests/api/utils.py | 180 ++++++++++++++++++++++++++++++++++++++++++--- 2 files changed, 191 insertions(+), 13 deletions(-) diff --git a/ocsmesh/utils.py b/ocsmesh/utils.py index 6ad2f2d9..63ecaa6d 100644 --- a/ocsmesh/utils.py +++ b/ocsmesh/utils.py @@ -2241,6 +2241,13 @@ def shape_to_msh_t(shape: Union[Polygon, MultiPolygon]) -> jigsaw_msh_t: vert2: List[Tuple[Tuple[float, float], int]] = [] if isinstance(shape, Polygon): shape = MultiPolygon([shape]) + + if not isinstance(shape, MultiPolygon): + raise ValueError(f"Invalid input shape type: {type(shape)}!") + + if not shape.is_valid: + raise ValueError("Input contains invalid (multi)polygons!") + for polygon in shape.geoms: if np.all( np.asarray( @@ -2275,14 +2282,23 @@ def shape_to_msh_t(shape: Union[Polygon, MultiPolygon]) -> jigsaw_msh_t: return msht -def shape_to_msh_t_2(shape: Union[Polygon, MultiPolygon]) -> jigsaw_msh_t: +def shape_to_msh_t_2( + shape: Union[Polygon, MultiPolygon, gpd.GeoDataFrame, gpd.GeoSeries] +) -> jigsaw_msh_t: gdf_shape = shape - if not isinstance(shape, gpd.GeoDataFrame): + if isinstance(shape, gpd.GeoSeries): + gdf_shape = gpd.GeoDataFrame(geometry=shape) + elif not isinstance(shape, gpd.GeoDataFrame): gdf_shape = gpd.GeoDataFrame(geometry=[shape]) + if np.sum(gdf_shape.area) == 0: + raise ValueError("The shape must have an area, such as Polygon!") + + if not np.all(gdf_shape.is_valid): + raise ValueError("Input contains invalid (multi)polygons!") + df_lonlat = ( gdf_shape - .geometry .boundary .explode(ignore_index=True) .map(lambda i: i.coords) @@ -2342,7 +2358,7 @@ def shape_to_msh_t_2(shape: Union[Polygon, MultiPolygon]) -> jigsaw_msh_t: coordinates=df_coo[['lon', 'lat']].values, edges=df_cnn.values, values=np.zeros((len(df_coo) ,1)), - crs=None + crs=gdf_shape.crs ) return msht diff --git a/tests/api/utils.py b/tests/api/utils.py index dc10efd2..d4901e79 100644 --- a/tests/api/utils.py +++ b/tests/api/utils.py @@ -1,12 +1,13 @@ #! python -import unittest +import re import tempfile +import unittest from copy import deepcopy import numpy as np +import geopandas as gpd from jigsawpy import jigsaw_msh_t from pyproj import CRS -import re from shapely.geometry import ( Point, LineString, @@ -15,6 +16,7 @@ MultiPolygon, GeometryCollection, ) +from shapely.ops import polygonize from ocsmesh import Raster, utils @@ -667,22 +669,182 @@ def test_data_masking(self): class ShapeToMeshT(unittest.TestCase): - def test_input_output_types(self): - self.assertTrue(False) + def setUp(self): + self.valid_input_1 = box(0, 0, 1, 1) + self.valid_input_2 = gpd.GeoDataFrame( + geometry=[self.valid_input_1] + ) + self.valid_input_3 = gpd.GeoSeries(self.valid_input_1) + # NOTE: Hole touching boundary is still valid shape for shapely + self.valid_input_4 = Polygon( + [ + [0, 0], + [2, 2], + [4, 0], + [2, -2], + [0, 0], + ], + [ + [ + [0, 0], + [1, -0.5], + [2, 0], + [1, 0.5], + [0, 0] + ] + ] + ) + self.valid_input_5 = MultiPolygon( + [box(0, 0, 1, 1), box(2, 2, 3, 3)] + ) + + self.invalid_input_1 = Point(0, 0) + self.invalid_input_2 = LineString([[0, 0], [1, 1]]) + # NOTE: Unlike hole touching boundary, this is invalid shape!! + self.invalid_input_3 = Polygon( + [ + [0, 0], + [2, 2], + [4, 0], + [2, -2], + [0, 0], + [1, -1], + [2, 0], + [1, 1], + [0, 0] + ] + ) + + + def test_old_input_output_validity(self): + msht = utils.shape_to_msh_t(self.valid_input_1) + with self.assertRaises(ValueError) as exc_1: + utils.shape_to_msh_t(self.invalid_input_1) + + with self.assertRaises(ValueError) as exc_2: + utils.shape_to_msh_t(self.invalid_input_2) + + with self.assertRaises(ValueError) as exc_3: + utils.shape_to_msh_t(self.invalid_input_3) + + self.assertIsInstance(msht, jigsaw_msh_t) + + self.assertIsNotNone( + re.search('invalid.*type', str(exc_1.exception).lower()), + ) + + self.assertIsNotNone( + re.search('invalid.*type', str(exc_2.exception).lower()), + ) + + self.assertIsNotNone( + re.search('invalid.*polygon', str(exc_3.exception).lower()), + ) + + + def test_new_input_output_validity(self): + msht_1 = utils.shape_to_msh_t_2(self.valid_input_1) + msht_2 = utils.shape_to_msh_t_2(self.valid_input_2) + msht_3 = utils.shape_to_msh_t_2(self.valid_input_3) + + with self.assertRaises(ValueError) as exc_1: + utils.shape_to_msh_t_2(self.invalid_input_1) + + with self.assertRaises(ValueError) as exc_2: + utils.shape_to_msh_t_2(self.invalid_input_2) + + with self.assertRaises(ValueError) as exc_3: + utils.shape_to_msh_t_2(self.invalid_input_3) + + self.assertIsInstance(msht_1, jigsaw_msh_t) + self.assertIsInstance(msht_2, jigsaw_msh_t) + self.assertIsInstance(msht_3, jigsaw_msh_t) + + self.assertTrue( + np.all(msht_1.vert2 == msht_2.vert2) + & np.all(msht_2.vert2 == msht_3.vert2) + ) + self.assertTrue( + np.all(msht_1.edge2 == msht_2.edge2) + & np.all(msht_2.edge2 == msht_3.edge2) + ) + self.assertTrue( + np.all(msht_1.value == msht_2.value) + & np.all(msht_2.value == msht_3.value) + ) + + self.assertIsNotNone( + re.search('have.*area', str(exc_1.exception).lower()), + ) + + self.assertIsNotNone( + re.search('have.*area', str(exc_2.exception).lower()), + ) + + self.assertIsNotNone( + re.search('invalid.*polygon', str(exc_3.exception).lower()), + ) def test_old_implementation(self): - self.assertTrue(False) + msht_1 = utils.shape_to_msh_t(self.valid_input_1) + msht_2 = utils.shape_to_msh_t(self.valid_input_4) + msht_3 = utils.shape_to_msh_t(self.valid_input_5) + + self.assertTrue( + list( + polygonize(msht_1.vert2['coord'][msht_1.edge2['index']]) + )[0].equals(self.valid_input_1) + ) + self.assertTrue( + list( + polygonize(msht_2.vert2['coord'][msht_2.edge2['index']]) + )[0].equals(self.valid_input_4) + ) + self.assertTrue( + MultiPolygon( + polygonize(msht_3.vert2['coord'][msht_3.edge2['index']]) + ).equals(self.valid_input_5) + ) + + # Old approach result in duplicate nodes! + self.assertEqual( + len(msht_2.vert2['coord']) - 1, + len(np.unique(msht_2.vert2['coord'], axis=0)) + ) def test_new_implementation(self): - self.assertTrue(False) + msht_1 = utils.shape_to_msh_t_2(self.valid_input_1) + msht_2 = utils.shape_to_msh_t_2(self.valid_input_4) + msht_3 = utils.shape_to_msh_t_2(self.valid_input_5) + self.assertTrue( + list( + polygonize(msht_1.vert2['coord'][msht_1.edge2['index']]) + )[0].equals(self.valid_input_1) + ) + self.assertTrue( + list( + polygonize(msht_2.vert2['coord'][msht_2.edge2['index']]) + )[0].equals(self.valid_input_4) + ) + self.assertTrue( + MultiPolygon( + polygonize(msht_3.vert2['coord'][msht_3.edge2['index']]) + ).equals(self.valid_input_5) + ) + + # New approach removes duplicate nodes! + self.assertEqual( + len(msht_2.vert2['coord']), + len(np.unique(msht_2.vert2['coord'], axis=0)) + ) -class TestTriangulatePolygon(unittest.TestCase): +class TestTriangulatePolygon(unittest.TestCase): def test_different_input_types(self): self.assertTrue(False) @@ -700,11 +862,11 @@ def test_fixed_boundary(self): self.assertTrue(False) - def test_aux_points(self) + def test_aux_points(self): self.assertTrue(False) - def test_polygon_holes(self) + def test_polygon_holes(self): self.assertTrue(False) From 186e5e2609892ae8c7935ab990bb9d94a5d98cba Mon Sep 17 00:00:00 2001 From: "Soroosh.Mani" Date: Fri, 15 Sep 2023 16:56:57 -0400 Subject: [PATCH 12/26] Add tests for triangulation --- ocsmesh/utils.py | 3 + tests/api/utils.py | 183 ++++++++++++++++++++++++++++++++++++++++++--- 2 files changed, 174 insertions(+), 12 deletions(-) diff --git a/ocsmesh/utils.py b/ocsmesh/utils.py index 63ecaa6d..75a399d7 100644 --- a/ocsmesh/utils.py +++ b/ocsmesh/utils.py @@ -2419,5 +2419,8 @@ def triangulate_polygon( values=np.zeros((len(cdt['vertices']) ,1)), crs=None, ) + if aux_pts is not None: + # To make sure unused points are discarded + cleanup_isolates(msht_tri) return msht_tri diff --git a/tests/api/utils.py b/tests/api/utils.py index d4901e79..c73d4bad 100644 --- a/tests/api/utils.py +++ b/tests/api/utils.py @@ -716,7 +716,7 @@ def setUp(self): ) - def test_old_input_output_validity(self): + def test_old_io_validity(self): msht = utils.shape_to_msh_t(self.valid_input_1) with self.assertRaises(ValueError) as exc_1: utils.shape_to_msh_t(self.invalid_input_1) @@ -742,7 +742,7 @@ def test_old_input_output_validity(self): ) - def test_new_input_output_validity(self): + def test_new_io_validity(self): msht_1 = utils.shape_to_msh_t_2(self.valid_input_1) msht_2 = utils.shape_to_msh_t_2(self.valid_input_2) msht_3 = utils.shape_to_msh_t_2(self.valid_input_3) @@ -844,31 +844,190 @@ def test_new_implementation(self): -class TestTriangulatePolygon(unittest.TestCase): +class TriangulatePolygon(unittest.TestCase): - def test_different_input_types(self): - self.assertTrue(False) + def setUp(self): + self.valid_input_1 = box(0, 0, 1, 1) + self.valid_input_2 = gpd.GeoDataFrame( + geometry=[self.valid_input_1] + ) + self.valid_input_3 = gpd.GeoSeries(self.valid_input_1) + # NOTE: Hole touching boundary is still valid shape for shapely + self.valid_input_4 = Polygon( + [ + [0, 0], + [2, 2], + [4, 0], + [2, -2], + [0, 0], + ], + [ + [ + [0, 0], + [1, -0.5], + [2, 0], + [1, 0.5], + [0, 0] + ] + ] + ) + self.valid_input_5 = MultiPolygon( + [box(0, 0, 1, 1), box(2, 2, 3, 3)] + ) - def test_wrong_input_of_right_type(self): - self.assertTrue(False) + self.invalid_input_1 = Point(0, 0) + self.invalid_input_2 = LineString([[0, 0], [1, 1]]) + # NOTE: Unlike hole touching boundary, this is invalid shape!! + self.invalid_input_3 = Polygon( + [ + [0, 0], + [2, 2], + [4, 0], + [2, -2], + [0, 0], + [1, -1], + [2, 0], + [1, 1], + [0, 0] + ] + ) - def test_output_types(self): - self.assertTrue(False) + def test_io_validity(self): + msht_1 = utils.triangulate_polygon(self.valid_input_1) + msht_2 = utils.triangulate_polygon(self.valid_input_2) + msht_3 = utils.triangulate_polygon(self.valid_input_3) + with self.assertRaises(ValueError) as exc_1: + utils.triangulate_polygon(self.invalid_input_1) - def test_fixed_boundary(self): - self.assertTrue(False) + with self.assertRaises(ValueError) as exc_2: + utils.triangulate_polygon(self.invalid_input_2) + + with self.assertRaises(ValueError) as exc_3: + utils.triangulate_polygon(self.invalid_input_3) + + self.assertIsInstance(msht_1, jigsaw_msh_t) + self.assertIsInstance(msht_2, jigsaw_msh_t) + self.assertIsInstance(msht_3, jigsaw_msh_t) + + self.assertTrue( + np.all(msht_1.vert2 == msht_2.vert2) + & np.all(msht_2.vert2 == msht_3.vert2) + ) + self.assertTrue( + np.all(msht_1.tria3 == msht_2.tria3) + & np.all(msht_2.tria3 == msht_3.tria3) + ) + self.assertTrue( + np.all(msht_1.value == msht_2.value) + & np.all(msht_2.value == msht_3.value) + ) + + self.assertEqual(len(msht_1.edge2), 0) + self.assertEqual(len(msht_1.quad4), 0) + self.assertEqual(len(msht_1.vert2), len(msht_1.value)) + self.assertTrue((msht_1.value == 0).all()) + + self.assertIsNotNone( + re.search( + 'must be convertible to polygon', + str(exc_1.exception).lower() + ), + ) + + self.assertIsNotNone( + re.search( + 'must be convertible to polygon', + str(exc_2.exception).lower() + ), + ) + + self.assertIsNotNone( + re.search('invalid.*polygon', str(exc_3.exception).lower()), + ) + + + def test_preserves_boundaries(self): + bx = Polygon( + np.array([ + [0, 0], [0, 1], [1, 1], [2, 1], [3, 1], + [4, 1], [4, 0], [3, 0], [2, 0], [1, 0] + ]) + + np.random.random((10, 2)) * 0.49 # to avoid exactness! + ) + msht = utils.triangulate_polygon(bx) + bdry_lines = utils.get_boundary_segments(msht) + + self.assertTrue( + list(polygonize(bdry_lines))[0].equals(bx) + ) def test_aux_points(self): - self.assertTrue(False) + bx = Polygon( + np.array([ + [0, 0], [0, 1], [1, 1], [2, 1], [3, 1], + [4, 1], [4, 0], [3, 0], [2, 0], [1, 0] + ]) + + np.random.random((10, 2)) * 0.49 # to avoid exactness! + ) + aux_1 = [[1, 0.5], [2, 0.5], [3, 0.5]] + aux_2 = [*aux_1, [10, 0.5]] # Out of domain points + + msht_1 = utils.triangulate_polygon(bx, aux_1, opts='p') + msht_2 = utils.triangulate_polygon(bx, aux_2, opts='p') + + self.assertTrue( + np.all([ + np.any([pt == v.tolist() for v in msht_1.vert2['coord']]) + for pt in aux_1 + ]) + ) + # Out of domain points are discarded + self.assertFalse( + np.all([ + np.any([pt == v.tolist() for v in msht_2.vert2['coord']]) + for pt in aux_2 + ]) + ) def test_polygon_holes(self): + poly = Polygon( + [[0, 0], [4, 0], [4, 4], [0, 4]], + [[[1, 1], [2, 0], [2, 2], [1, 2]]] + ) + msht = utils.triangulate_polygon(poly, opts='p') + mesh_poly = utils.get_mesh_polygons(msht) + + self.assertTrue(poly.equals(mesh_poly)) + + + def test_multipolygon(self): + mpoly = MultiPolygon([box(0, 0, 1, 1), box(2, 2, 3, 3)]) + + msht = utils.triangulate_polygon(mpoly, opts='p') + mesh_poly = utils.get_mesh_polygons(msht) + + self.assertTrue(mpoly.equals(mesh_poly)) + + + + +class MeshTFromNumpy(unittest.TestCase): + + def test_io_validity(self): self.assertTrue(False) + def test_values(self): + self.assertTrue(False) + + def test_crs(self): + self.assertTrue(False) + + if __name__ == '__main__': unittest.main() From 132e1141be6f9f702faaf0ec2227530908fc62f1 Mon Sep 17 00:00:00 2001 From: "Soroosh.Mani" Date: Wed, 20 Sep 2023 21:21:12 -0400 Subject: [PATCH 13/26] Fix triangulation and mesh poly calc --- ocsmesh/utils.py | 103 +++++++++------------------------------------ tests/api/utils.py | 25 +++++++---- 2 files changed, 36 insertions(+), 92 deletions(-) diff --git a/ocsmesh/utils.py b/ocsmesh/utils.py index 75a399d7..389f27b4 100644 --- a/ocsmesh/utils.py +++ b/ocsmesh/utils.py @@ -193,86 +193,14 @@ def get_boundary_segments(mesh): def get_mesh_polygons(mesh): + elm_polys = [] + for elm_type in ELEM_2D_TYPES: + elems = getattr(mesh, elm_type)['index'] + elm_polys.extend( + [Polygon(mesh.vert2['coord'][cell]) for cell in elems] + ) - - # TODO: Copy mesh? - target_mesh = mesh - result_polys = [] - - # 2-pass find, first find using polygons that intersect non-boundary - # vertices, then from the rest of the mesh find polygons that - # intersect any vertex - for find_pass in range(2): - - - coords = target_mesh.vert2['coord'] - - if len(coords) == 0: - continue - - boundary_edges = get_boundary_edges(target_mesh) - - lines = get_boundary_segments(target_mesh) - - poly_gen = polygonize(lines) - polys = list(poly_gen) - polys = sorted(polys, key=lambda p: p.area, reverse=True) - - - bndry_verts = np.unique(boundary_edges) - - if find_pass == 0: - non_bndry_verts = np.setdiff1d( - np.arange(len(coords)), bndry_verts) - pnts = MultiPoint(coords[non_bndry_verts]) - else: - pnts = MultiPoint(coords[bndry_verts]) - - - # NOTE: This logic requires polygons to be sorted by area - pass_valid_polys = [] - while len(pnts.geoms) > 0: - - - idx = np.random.randint(len(pnts.geoms)) - pnt = pnts.geoms[idx] - - polys_gdf = gpd.GeoDataFrame( - {'geometry': polys, 'list_index': range(len(polys))}) - - - res_gdf = polys_gdf[polys_gdf.intersects(pnt)] - if len(res_gdf) == 0: - # 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 - - continue - - poly = res_gdf.geometry.iloc[0] - polys.pop(res_gdf.iloc[0].list_index) - - - - pass_valid_polys.append(poly) - pnts = pnts.difference(poly) - if pnts.is_empty: - break - if isinstance(pnts, Point): - pnts = MultiPoint([pnts]) - - - result_polys.extend(pass_valid_polys) - target_mesh = clip_mesh_by_shape( - target_mesh, - shape=MultiPolygon(pass_valid_polys), - inverse=True, fit_inside=True) - - - - return MultiPolygon(result_polys) + return unary_union(elm_polys) def repartition_features(linestring, max_verts): @@ -2390,12 +2318,19 @@ def triangulate_polygon( coords = msht_shape.vert2['coord'] edges = msht_shape.edge2['index'] + # To make sure islands formed by two polygons touching at multiple + # points are also considered as holes, instead of getting interiors, + # we get the negative of the domain and calculate points in all + # negative polygons! + + # buffer by 1/100 of shape length scale = sqrt(area)/100 + world = box(*shape.buffer(np.sqrt(shape.area) / 100).bounds) + neg_shape = world.difference(shape) + if isinstance(neg_shape, Polygon): + neg_shape = MultiPolygon([neg_shape]) holes = [] - for poly in shape.geoms: - for ring in poly.interiors: - holes.append( - np.array(Polygon(ring).representative_point().xy).ravel() - ) + for poly in neg_shape.geoms: + holes.append(np.array(poly.representative_point().xy).ravel()) holes = np.array(holes) diff --git a/tests/api/utils.py b/tests/api/utils.py index c73d4bad..cb96a507 100644 --- a/tests/api/utils.py +++ b/tests/api/utils.py @@ -597,10 +597,12 @@ def test_values(self): # Test value input of wrong size # Test value None input self.assertTrue(False) + #TODO: def test_kwonly_args(self): self.assertTrue(False) + #TODO: class CreateRasterFromNumpy(unittest.TestCase): @@ -960,6 +962,8 @@ def test_preserves_boundaries(self): msht = utils.triangulate_polygon(bx) bdry_lines = utils.get_boundary_segments(msht) + # TODO: Make sure equal means equal in all vertices, not just + # combined shape (i.e. no edge split) self.assertTrue( list(polygonize(bdry_lines))[0].equals(bx) ) @@ -1014,17 +1018,22 @@ def test_multipolygon(self): self.assertTrue(mpoly.equals(mesh_poly)) + def test_polygons_touching_two_points_no_hole(self): + poly1 = Polygon( + [[0, 0], [0, 4], [6, 4], [6, 0], [4, 2], [2, 2], [0, 0]], + ) + poly2 = Polygon( + [[0, 0], [0, -4], [6, -4], [6, 0], [4, -2], [2, -2], [0, 0]], + ) + multpoly = MultiPolygon([poly1, poly2]) + msht = utils.triangulate_polygon(multpoly, opts='p') + mesh_poly = utils.get_mesh_polygons(msht) + self.assertTrue(multpoly.equals(mesh_poly)) -class MeshTFromNumpy(unittest.TestCase): - - def test_io_validity(self): - self.assertTrue(False) - - def test_values(self): - self.assertTrue(False) - def test_crs(self): +class GetMeshPolygon(unittest.TestCase): + def test_something(self): self.assertTrue(False) From 6dbd530033a420f539d39faf858d3183151da7cc Mon Sep 17 00:00:00 2001 From: "Soroosh.Mani" Date: Wed, 20 Sep 2023 21:37:23 -0400 Subject: [PATCH 14/26] Bugfix for new mesh polygon calc --- ocsmesh/utils.py | 6 +++++- tests/api/utils.py | 2 +- 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/ocsmesh/utils.py b/ocsmesh/utils.py index 389f27b4..d08d319b 100644 --- a/ocsmesh/utils.py +++ b/ocsmesh/utils.py @@ -200,7 +200,11 @@ def get_mesh_polygons(mesh): [Polygon(mesh.vert2['coord'][cell]) for cell in elems] ) - return unary_union(elm_polys) + poly = unary_union(elm_polys) + if isinstance(poly, Polygon): + poly = MultiPolygon([poly]) + + return poly def repartition_features(linestring, max_verts): diff --git a/tests/api/utils.py b/tests/api/utils.py index cb96a507..eb737ca7 100644 --- a/tests/api/utils.py +++ b/tests/api/utils.py @@ -1033,7 +1033,7 @@ def test_polygons_touching_two_points_no_hole(self): class GetMeshPolygon(unittest.TestCase): - def test_something(self): + def test_always_returns_multipolygon(self): self.assertTrue(False) From c80221e55de9d17a6c42c54353c328fff36e362b Mon Sep 17 00:00:00 2001 From: Parallel Works app-run user Date: Thu, 21 Sep 2023 20:33:56 +0000 Subject: [PATCH 15/26] Fix subsetting logic for bufferzone meshing --- ocsmesh/cli/subset_n_combine.py | 31 +++++++++++++++++++------------ 1 file changed, 19 insertions(+), 12 deletions(-) diff --git a/ocsmesh/cli/subset_n_combine.py b/ocsmesh/cli/subset_n_combine.py index 94d7a220..b0955bb4 100755 --- a/ocsmesh/cli/subset_n_combine.py +++ b/ocsmesh/cli/subset_n_combine.py @@ -13,7 +13,7 @@ from pyproj import CRS, Transformer from scipy.interpolate import griddata from scipy.spatial import cKDTree -from shapely.geometry import MultiPolygon, Polygon, GeometryCollection +from shapely.geometry import MultiPolygon, Polygon, GeometryCollection, MultiPoint from shapely.ops import polygonize, unary_union, transform from ocsmesh import Raster, Geom, Mesh, Hfun, utils, JigsawDriver @@ -325,23 +325,30 @@ def _generate_mesh_for_buffer_region( transformer = Transformer.from_crs(buffer_crs, utm, always_xy=True) buffer_polygon = transform(transformer.transform, buffer_polygon) - mesh_domain_interp = JigsawDriver( + mesh_buf_apprx = JigsawDriver( geom=Geom(buffer_polygon, crs=utm), hfun=hfun_buffer, initial_mesh=False ).run(sieve=0) - msht_domain_interp = deepcopy(mesh_domain_interp.msh_t) - in_verts_mask = np.ones_like(msht_domain_interp.vert2, dtype=bool) - in_verts_mask[ - np.unique(utils.get_boundary_edges(msht_domain_interp)) - ] = False - -# utils.reproject(msht_domain_interp, buffer_crs) + msht_buf_apprx = deepcopy(mesh_buf_apprx.msh_t) + # If vertices are too close to buffer geom boundary, + # it's going to cause issues (thin elements) + if msht_buf_apprx.crs != hfun_buffer.crs: + ocsmesh.utils.reproject(msht_buf_apprx, hfun_buffer.crs) + gdf_pts = gpd.GeoDataFrame( + geometry=[MultiPoint(msht_buf_apprx.vert2['coord'])], + crs=msht_buf_apprx.crs + ).explode() + gdf_aux_pts = gdf_pts[ + (~gdf_pts.intersects( + buffer_polygon.boundary.buffer(hfun_buffer.hmin)) + ) & (gdf_pts.within(buffer_polygon)) + ] + +# utils.reproject(msht_buf_apprx, buffer_crs) msht_buffer = utils.triangulate_polygon( - shape=buffer_polygon, - aux_pts=msht_domain_interp.vert2[in_verts_mask]['coord'], - opts='p' + shape=buffer_polygon, aux_pts=gdf_aux_pts, opts='p' ) msht_buffer.crs = utm From cc0b82a245754dc57997172675636b5ee93c42f1 Mon Sep 17 00:00:00 2001 From: "Soroosh.Mani" Date: Thu, 21 Sep 2023 20:12:48 -0400 Subject: [PATCH 16/26] Fix getting aux point XY for triangulation --- ocsmesh/utils.py | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/ocsmesh/utils.py b/ocsmesh/utils.py index d08d319b..40e7b9f1 100644 --- a/ocsmesh/utils.py +++ b/ocsmesh/utils.py @@ -2339,11 +2339,14 @@ def triangulate_polygon( if aux_pts is not None: - if isinstance(aux_pts, (gpd.GeoDataFrame, gpd.GeoSeries)): - aux_pts = np.array(aux_pts.unary_union.xy).T - - elif isinstance(aux_pts, (Point, MultiPolygon)): - aux_pts = np.array(aux_pts.xy).T + if isinstance(aux_pts, (Point, MultiPolygon)): + aux_pts = gpd.GeoSeries(aux_pts) + elif isinstance(aux_pts, gpd.GeoDataFrame): + aux_pts = aux_pts.geometry + elif not isinstance(aux_pts, gpd.GeoSeries): + raise ValueError("Wrong input type for !") + + aux_pts = aux_pts.get_coordinates().values coords = np.vstack((coords, aux_pts)) From eed2a7054f3b2f1a19f7bbbb93cb40f3aabcbaf0 Mon Sep 17 00:00:00 2001 From: "Soroosh.Mani" Date: Mon, 25 Sep 2023 08:50:00 -0400 Subject: [PATCH 17/26] Revert python ver restriction due to triangle pkg --- .github/workflows/documentation.yml | 2 -- .github/workflows/functional_test.yml | 2 +- .github/workflows/functional_test_2.yml | 2 +- .github/workflows/pylint.yml | 2 -- pyproject.toml | 2 +- 5 files changed, 3 insertions(+), 7 deletions(-) diff --git a/.github/workflows/documentation.yml b/.github/workflows/documentation.yml index bead88e8..42131d95 100644 --- a/.github/workflows/documentation.yml +++ b/.github/workflows/documentation.yml @@ -22,8 +22,6 @@ jobs: with: environment-file: ./environment.yml environment-name: ocsmesh-env - extra-specs: | - python=3.10 - name: Install dependencies shell: bash -l {0} run: | diff --git a/.github/workflows/functional_test.yml b/.github/workflows/functional_test.yml index 728c445e..661c50b4 100644 --- a/.github/workflows/functional_test.yml +++ b/.github/workflows/functional_test.yml @@ -28,7 +28,7 @@ jobs: fail-fast: false matrix: os: [ ubuntu-latest ] - python-version: [ '3.9', '3.10' ] + python-version: [ '3.9', '3.10', '3.x'] steps: - uses: actions/checkout@v3 diff --git a/.github/workflows/functional_test_2.yml b/.github/workflows/functional_test_2.yml index fe6a7182..4fc21ad1 100644 --- a/.github/workflows/functional_test_2.yml +++ b/.github/workflows/functional_test_2.yml @@ -28,7 +28,7 @@ jobs: fail-fast: false matrix: os: [ ubuntu-latest ] - python-version: [ '3.9', '3.10' ] + python-version: [ '3.9', '3.10', '3.x' ] steps: - name: Checkout diff --git a/.github/workflows/pylint.yml b/.github/workflows/pylint.yml index 12ee8bb5..3d8589f6 100644 --- a/.github/workflows/pylint.yml +++ b/.github/workflows/pylint.yml @@ -30,8 +30,6 @@ jobs: with: environment-file: ./environment.yml environment-name: ocsmesh-env - extra-specs: | - python=3.10 - name: Install dependencies shell: bash -l {0} run: | diff --git a/pyproject.toml b/pyproject.toml index d0886b8c..e68b9bdc 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -19,7 +19,7 @@ maintainers = [ description = "Package to generate computational unstructured meshes from planetary modeling." license = {file = "LICENSE"} readme = "README.md" -requires-python = '>=3.9, <3.11' # 3.8 -> scipy, 3.11 -> triangle +requires-python = '>=3.9' # 3.8 -> scipy dependencies = [ "colored-traceback", "fiona", "geopandas", "jigsawpy", "matplotlib", "netCDF4", "numba", From 575c993fb691aa87bfa1184526445979e9bfc2c8 Mon Sep 17 00:00:00 2001 From: "Soroosh.Mani" Date: Mon, 25 Sep 2023 09:50:54 -0400 Subject: [PATCH 18/26] Fix order of bbox points for hfun calc --- ocsmesh/hfun/raster.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/ocsmesh/hfun/raster.py b/ocsmesh/hfun/raster.py index 77b943ae..18430e5e 100644 --- a/ocsmesh/hfun/raster.py +++ b/ocsmesh/hfun/raster.py @@ -383,10 +383,11 @@ def msh_t( transformer = Transformer.from_crs( self.crs, utm_crs, always_xy=True) bbox = [ - *[(x, left[0]) for x in bottom], - *[(bottom[-1], y) for y in reversed(right)], - *[(x, right[-1]) for x in reversed(top)], - *[(bottom[0], y) for y in reversed(left)]] + *[(x, left[0]) for x in bottom][:-1], + *[(bottom[-1], y) for y in right][:-1], + *[(x, right[-1]) for x in reversed(top)][:-1], + *[(bottom[0], y) for y in reversed(left)][:-1] + ] geom = PolygonGeom( ops.transform(transformer.transform, Polygon(bbox)), utm_crs From 23a7e90abb9839a06d299ff7adeaf7c78659d67f Mon Sep 17 00:00:00 2001 From: "Soroosh.Mani" Date: Mon, 25 Sep 2023 10:56:13 -0400 Subject: [PATCH 19/26] Fix mamba argument for py ver --- .github/workflows/functional_test.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/functional_test.yml b/.github/workflows/functional_test.yml index 661c50b4..6152eee4 100644 --- a/.github/workflows/functional_test.yml +++ b/.github/workflows/functional_test.yml @@ -28,7 +28,7 @@ jobs: fail-fast: false matrix: os: [ ubuntu-latest ] - python-version: [ '3.9', '3.10', '3.x'] + python-version: [ '3.9', '3.10', '3.11'] steps: - uses: actions/checkout@v3 From 1dfffeb52ff8e267d05d8aa36e9aab8736f55195 Mon Sep 17 00:00:00 2001 From: "Soroosh.Mani" Date: Mon, 25 Sep 2023 10:57:03 -0400 Subject: [PATCH 20/26] Fix for pylint --- ocsmesh/cli/subset_n_combine.py | 10 ++++----- ocsmesh/geom/base.py | 3 +-- ocsmesh/utils.py | 38 +++++++++++++++++++++++++++------ 3 files changed, 36 insertions(+), 15 deletions(-) diff --git a/ocsmesh/cli/subset_n_combine.py b/ocsmesh/cli/subset_n_combine.py index b0955bb4..76d0f0a8 100755 --- a/ocsmesh/cli/subset_n_combine.py +++ b/ocsmesh/cli/subset_n_combine.py @@ -6,9 +6,7 @@ from time import time import geopandas as gpd -import jigsawpy from jigsawpy.msh_t import jigsaw_msh_t -from jigsawpy.jig_t import jigsaw_jig_t import numpy as np from pyproj import CRS, Transformer from scipy.interpolate import griddata @@ -251,7 +249,7 @@ def _calculate_mesh_size_function( buffer_crs ): - assert(buffer_crs == hires_mesh_clip.crs == lores_mesh_clip.crs) + assert buffer_crs == hires_mesh_clip.crs == lores_mesh_clip.crs # HARDCODED FOR NOW approx_elem_per_width = 3 @@ -283,7 +281,7 @@ def _calculate_mesh_size_function( hfun_cdt.size_from_mesh() hfun_cdt_sz = deepcopy(hfun_cdt.msh_t().value) / approx_elem_per_width - msht_cdt.value[:] = hfun_cdt_sz + msht_cdt.value[:] = hfun_cdt_sz hfun_rep = Hfun(Mesh(msht_cdt)) mesh_domain_rep = JigsawDriver( @@ -332,10 +330,10 @@ def _generate_mesh_for_buffer_region( ).run(sieve=0) msht_buf_apprx = deepcopy(mesh_buf_apprx.msh_t) - # If vertices are too close to buffer geom boundary, + # If vertices are too close to buffer geom boundary, # it's going to cause issues (thin elements) if msht_buf_apprx.crs != hfun_buffer.crs: - ocsmesh.utils.reproject(msht_buf_apprx, hfun_buffer.crs) + utils.reproject(msht_buf_apprx, hfun_buffer.crs) gdf_pts = gpd.GeoDataFrame( geometry=[MultiPoint(msht_buf_apprx.vert2['coord'])], crs=msht_buf_apprx.crs diff --git a/ocsmesh/geom/base.py b/ocsmesh/geom/base.py index 05091cea..581257b1 100644 --- a/ocsmesh/geom/base.py +++ b/ocsmesh/geom/base.py @@ -2,10 +2,9 @@ """ from abc import ABC, abstractmethod -from typing import List, Tuple, Any, Union +from typing import Any, Union from jigsawpy import jigsaw_msh_t -import numpy as np from pyproj import CRS, Transformer from shapely import ops from shapely.geometry import MultiPolygon diff --git a/ocsmesh/utils.py b/ocsmesh/utils.py index 40e7b9f1..5991a0b7 100644 --- a/ocsmesh/utils.py +++ b/ocsmesh/utils.py @@ -1,6 +1,6 @@ from collections import defaultdict from itertools import permutations -from typing import Union, Dict, Sequence, Tuple +from typing import Union, Dict, Sequence, Tuple, List from functools import reduce from multiprocessing import cpu_count, Pool from copy import deepcopy @@ -96,7 +96,7 @@ def cleanup_duplicates(mesh): return_index=True, return_inverse=True ) - nd_map = {e: i for e, i in enumerate(coorev)} + nd_map = dict(enumerate(coorev)) mesh.vert2 = mesh.vert2.take(cooidx, axis=0) for etype, otype in MESH_TYPES.items(): @@ -2275,7 +2275,7 @@ def shape_to_msh_t_2( ) ) - + ar_edg = np.sort(df_edg.values, axis=1) df_cnn = ( pd.DataFrame.from_records( @@ -2285,7 +2285,7 @@ def shape_to_msh_t_2( .drop_duplicates() # Remove duplicate edges .reset_index(drop=True) ) - + msht = msht_from_numpy( coordinates=df_coo[['lon', 'lat']].values, edges=df_cnn.values, @@ -2302,9 +2302,33 @@ def triangulate_polygon( aux_pts: Union[np.array, Point, MultiPoint, gpd.GeoDataFrame, gpd.GeoSeries] = None, opts='p', ) -> None: - ''' + """Triangulate the input shape, with additional points provided - ''' + Use `triangle` package to triangulate the input shape. If provided, + use the input additional points as well. The input `opts` controls + how `triangle` treats the inputs. + See the documentation for `triangle` more information + + Parameters + ---------- + shape : Polygon or MultiPolygon or GeoDataFrame or GeoSeries + Input polygon to triangulate + aux_pts : numpy array or Point or MultiPoint or GeoDataFrame or GeoSeries + Extra points to be used in the triangulation + opts : string, default='p' + Options for controlling `triangle` package + + Returns + ------- + jigsaw_msh_t + Generated triangulation + + Raises + ------ + ValueError + If input shape is invalid or the point input cannot be used + to generate point shape + """ # NOTE: Triangle can handle separate closed polygons, # so no need to have for loop @@ -2317,7 +2341,7 @@ def triangulate_polygon( if not isinstance(shape, (MultiPolygon, Polygon)): raise ValueError("Input shape must be convertible to polygon!") - + msht_shape = shape_to_msh_t_2(shape) coords = msht_shape.vert2['coord'] edges = msht_shape.edge2['index'] From be406b969604fc3b03acf3308840eb77765a9086 Mon Sep 17 00:00:00 2001 From: "Soroosh.Mani" Date: Mon, 25 Sep 2023 11:13:52 -0400 Subject: [PATCH 21/26] Fix triangulate aux point input processing --- ocsmesh/utils.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/ocsmesh/utils.py b/ocsmesh/utils.py index 5991a0b7..1e020273 100644 --- a/ocsmesh/utils.py +++ b/ocsmesh/utils.py @@ -2363,7 +2363,9 @@ def triangulate_polygon( if aux_pts is not None: - if isinstance(aux_pts, (Point, MultiPolygon)): + if isinstance(aux_pts, (np.ndarray, list, tuple)): + aux_pts = MultiPoint(aux_pts) + if isinstance(aux_pts, (Point, MultiPoint)): aux_pts = gpd.GeoSeries(aux_pts) elif isinstance(aux_pts, gpd.GeoDataFrame): aux_pts = aux_pts.geometry From 62bc5736b05339a4e6ff34c27e20d995bdefbad6 Mon Sep 17 00:00:00 2001 From: "Soroosh.Mani" Date: Mon, 25 Sep 2023 11:14:09 -0400 Subject: [PATCH 22/26] Remove numpy deprecated func --- ocsmesh/utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ocsmesh/utils.py b/ocsmesh/utils.py index 1e020273..65cefc08 100644 --- a/ocsmesh/utils.py +++ b/ocsmesh/utils.py @@ -1046,7 +1046,7 @@ def get_mesh_edges(mesh: jigsaw_msh_t, unique=True): (elm_type, np.roll(elm_type, shift=1, axis=1)), axis=2), axis=2) - edges = edges.reshape(np.product(edges.shape[0:2]), 2) + edges = edges.reshape(np.prod(edges.shape[0:2]), 2) all_edges = np.vstack((all_edges, edges)) if unique: From ab9e99d1dc017b65637695ca81bc2533b40bebc3 Mon Sep 17 00:00:00 2001 From: "Soroosh.Mani" Date: Mon, 25 Sep 2023 11:18:22 -0400 Subject: [PATCH 23/26] Update threshold for size from mesh func --- tests/api/hfun.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/api/hfun.py b/tests/api/hfun.py index 1f43f9f7..cfcb836b 100755 --- a/tests/api/hfun.py +++ b/tests/api/hfun.py @@ -483,8 +483,8 @@ def test_calculated_size(self): hfun_calc_val = hfun_calc_jig.value hfun_val_diff = self.hfun_orig_val - hfun_calc_val - # TODO: Come up with a more robust criteria - threshold = 0.2 + # TODO: Come up with a more robust criteria! + threshold = 0.21 err_value = np.max(np.abs(hfun_val_diff))/np.max(self.hfun_orig_val) self.assertTrue(err_value < threshold) From 42995a7306dc3db6d0febe8c24d0242fb83e8aba Mon Sep 17 00:00:00 2001 From: "Soroosh.Mani" Date: Mon, 25 Sep 2023 11:42:03 -0400 Subject: [PATCH 24/26] Fix for checking index ring equivalence --- tests/api/mesh.py | 24 ++++++++++++++++++++---- 1 file changed, 20 insertions(+), 4 deletions(-) diff --git a/tests/api/mesh.py b/tests/api/mesh.py index 1d4a0f68..0f2e8804 100644 --- a/tests/api/mesh.py +++ b/tests/api/mesh.py @@ -53,6 +53,18 @@ def setUp(self): self.mesh = Mesh(mesh_msht) + def _chkEqualPermutations(self, a, b): + # Make sure there are rings! + assert a[0] == a[-1] and b[0] == b[-1] + + a = np.array(a[:-1]) + b = np.array(b[:-1]) + idx = np.nonzero(a == b[0]) + if len(idx) == 0: + raise ValueError("One item in arg2 is not found in arg1") + shift = -idx[0].item() + self.assertTrue(np.all(np.roll(a, shift=shift) == np.array(b))) + def test_auto_boundary_fails_if_na_elev(self): # Set one node to nan value self.mesh.msh_t.value[-1] = np.nan @@ -79,7 +91,7 @@ def test_auto_boundary_1open_correctness(self): bdry.land().iloc[0]['index_id'], [5, 4, 3, 2, 1, 6, 11, 16, 21, 26, 27, 28, 29, 30] ) - self.assertEqual(bdry.interior().iloc[0]['index_id'], [12, 13, 18, 17, 12]) + self._chkEqualPermutations(bdry.interior().iloc[0]['index_id'], [12, 13, 18, 17, 12]) def test_auto_boundary_2open_correctness(self): @@ -100,7 +112,7 @@ def test_auto_boundary_2open_correctness(self): self.assertEqual(bdry.open().iloc[1]['index_id'], [30, 25, 20, 15, 10, 5]) self.assertEqual(bdry.land().iloc[0]['index_id'], [5, 4, 3, 2, 1]) self.assertEqual(bdry.land().iloc[1]['index_id'], [26, 27, 28, 29, 30]) - self.assertEqual(bdry.interior().iloc[0]['index_id'], [12, 13, 18, 17, 12]) + self._chkEqualPermutations(bdry.interior().iloc[0]['index_id'], [12, 13, 18, 17, 12]) def test_manual_boundary_specification_correctness(self): @@ -127,7 +139,7 @@ def test_manual_boundary_specification_correctness(self): [1, 6, 11, 16, 21, 26, 27, 28, 29, 30, 25, 20, 15, 10, 5] ) self.assertEqual(bdry.land().iloc[1]['index_id'], [2, 3, 4]) - self.assertEqual(bdry.interior().iloc[0]['index_id'], [12, 13, 18, 17, 12]) + self._chkEqualPermutations(bdry.interior().iloc[0]['index_id'], [12, 13, 18, 17, 12]) def test_manual_boundary_notaffect_interior(self): @@ -288,7 +300,7 @@ def test_auto_find_islands_only(self): bdry.find_islands() self.assertEqual(len(bdry.interior()), 1) - self.assertEqual(bdry.interior().iloc[0]['index_id'], [12, 13, 18, 17, 12]) + self._chkEqualPermutations(bdry.interior().iloc[0]['index_id'], [12, 13, 18, 17, 12]) def test_specify_boundary_on_mesh_with_no_boundary(self): @@ -303,3 +315,7 @@ def test_specify_boundary_on_mesh_with_no_boundary(self): self.assertEqual(len(bdry.interior()), 0) self.assertEqual(bdry.open().iloc[0]['index_id'], [1, 2, 3]) + + +if __name__ == '__main__': + unittest.main() From 21ccc3850433d54cea3b971a29fb2910acb6ea40 Mon Sep 17 00:00:00 2001 From: "Soroosh.Mani" Date: Mon, 25 Sep 2023 12:02:22 -0400 Subject: [PATCH 25/26] Test for values of msht_from_numpy --- ocsmesh/utils.py | 5 ++++- tests/api/utils.py | 42 ++++++++++++++++++++++++++++++++++-------- 2 files changed, 38 insertions(+), 9 deletions(-) diff --git a/ocsmesh/utils.py b/ocsmesh/utils.py index 65cefc08..896ba59f 100644 --- a/ocsmesh/utils.py +++ b/ocsmesh/utils.py @@ -2137,7 +2137,10 @@ def msht_from_numpy( np.zeros((len(mesh.vert2), 1)), dtype=jigsaw_msh_t.REALS_t ) - elif values.shape != (len(mesh.vert2), 1): + elif not isinstance(values, np.ndarray): + values = np.array(values).reshape(-1, 1) + + if values.shape != (len(mesh.vert2), 1): raise ValueError( "Input for mesh values must either be None or a" " 2D-array with row count equal to vertex count!" diff --git a/tests/api/utils.py b/tests/api/utils.py index eb737ca7..c943f438 100644 --- a/tests/api/utils.py +++ b/tests/api/utils.py @@ -592,17 +592,43 @@ def test_crs(self): self.assertEqual(out_msht_2.crs, CRS.from_user_input('esri:102008')) - def test_values(self): - # Test it has values - # Test value input of wrong size - # Test value None input - self.assertTrue(False) - #TODO: + def test_values_are_assigned(self): + out_msht = utils.msht_from_numpy( + coordinates=self.in_verts, + triangles=self.in_tria, + crs=None + ) + + self.assertTrue(len(out_msht.value) == len(self.in_verts)) + self.assertTrue(np.all(out_msht.value == 0)) + + def test_values_input_validation(self): + with self.assertRaises(ValueError) as exc_1: + utils.msht_from_numpy( + coordinates=self.in_verts, + triangles=self.in_tria, + values=[1,2,3], + crs=None + ) + + self.assertIsNotNone( + re.search( + 'values must either be', + str(exc_1.exception).lower() + ), + ) def test_kwonly_args(self): - self.assertTrue(False) - #TODO: + with self.assertRaises(Exception) as exc_1: + utils.msht_from_numpy(self.in_verts, self.in_tria) + + self.assertIsNotNone( + re.search( + 'takes 1 positional argument', + str(exc_1.exception).lower() + ), + ) class CreateRasterFromNumpy(unittest.TestCase): From e000811865c13c0a592fd2b8952c1329379f017d Mon Sep 17 00:00:00 2001 From: "Soroosh.Mani" Date: Mon, 25 Sep 2023 12:11:35 -0400 Subject: [PATCH 26/26] Add test for mesh polygon return type --- tests/api/utils.py | 15 ++++++++++++++- 1 file changed, 14 insertions(+), 1 deletion(-) diff --git a/tests/api/utils.py b/tests/api/utils.py index c943f438..9955a6ee 100644 --- a/tests/api/utils.py +++ b/tests/api/utils.py @@ -1060,9 +1060,22 @@ def test_polygons_touching_two_points_no_hole(self): class GetMeshPolygon(unittest.TestCase): def test_always_returns_multipolygon(self): - self.assertTrue(False) + poly1 = Polygon( + [[0, 0], [0, 4], [6, 4], [6, 0], [4, 2], [2, 2], [0, 0]], + ) + poly2 = Polygon( + [[0, 0], [0, -4], [6, -4], [6, 0], [4, -2], [2, -2], [0, 0]], + ) + multpoly = MultiPolygon([poly1, poly2]) + + msht_1 = utils.triangulate_polygon(poly1, opts='p') + msht_2 = utils.triangulate_polygon(multpoly, opts='p') + mesh_poly_1 = utils.get_mesh_polygons(msht_1) + mesh_poly_2 = utils.get_mesh_polygons(msht_2) + self.assertIsInstance(mesh_poly_1, MultiPolygon) + self.assertIsInstance(mesh_poly_2, MultiPolygon) if __name__ == '__main__': unittest.main()