diff --git a/.github/workflows/documentation.yml b/.github/workflows/documentation.yml index 98b5d66c..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.* - name: Install dependencies shell: bash -l {0} run: | diff --git a/.github/workflows/functional_test.yml b/.github/workflows/functional_test.yml index 6a4f753f..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.11' ] + python-version: [ '3.9', '3.10', '3.11'] steps: - uses: actions/checkout@v3 diff --git a/.github/workflows/functional_test_2.yml b/.github/workflows/functional_test_2.yml index 92153e6d..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', '3.11' ] + python-version: [ '3.9', '3.10', '3.x' ] steps: - name: Checkout diff --git a/.github/workflows/pylint.yml b/.github/workflows/pylint.yml index 3a8b0e20..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.* - name: Install dependencies shell: bash -l {0} run: | diff --git a/ocsmesh/cli/subset_n_combine.py b/ocsmesh/cli/subset_n_combine.py index be863e56..76d0f0a8 100755 --- a/ocsmesh/cli/subset_n_combine.py +++ b/ocsmesh/cli/subset_n_combine.py @@ -6,16 +6,15 @@ 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 from scipy.spatial import cKDTree -from shapely.geometry import MultiPolygon, Polygon, GeometryCollection -from shapely.ops import polygonize, unary_union +from shapely.geometry import MultiPolygon, Polygon, GeometryCollection, MultiPoint +from shapely.ops import polygonize, unary_union, transform -from ocsmesh import Raster, Geom, Mesh, Hfun, utils +from ocsmesh import Raster, Geom, Mesh, Hfun, utils, JigsawDriver logging.basicConfig( stream=sys.stdout, @@ -242,96 +241,118 @@ 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 + + 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) - 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 + # 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 = utm - return jig_hfun +# utils.reproject(msht_cdt, utm) + hfun_cdt = Hfun(Mesh(msht_cdt)) + hfun_cdt.size_from_mesh() + 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)) - def _generate_mesh_for_buffer_region( - self, buffer_polygon, jig_hfun, crs): + mesh_domain_rep = JigsawDriver( + geom=Geom(buffer_domain, crs=utm), + hfun=hfun_rep, + initial_mesh=False + ).run(sieve=0) - seam = Geom(buffer_polygon, crs=crs) + msht_domain_rep = deepcopy(mesh_domain_rep.msh_t) + utils.reproject(msht_domain_rep, utm) - jig_geom = seam.msh_t() + 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)] - # IMPORTANT: Setting these to -1 results in NON CONFORMAL boundary -# jig_geom.vert2['IDtag'][:] = -1 -# jig_geom.edge2['IDtag'][:] = -1 + msht_domain_rep.value[:] = domain_sz + hfun_interp = Hfun(Mesh(msht_domain_rep)) - jig_init = deepcopy(seam.msh_t()) - jig_init.vert2['IDtag'][:] = -1 - jig_init.edge2['IDtag'][:] = -1 + return hfun_interp - # 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() - # TODO: Use marche to calculate mesh size in the area between - # the two regions? + def _generate_mesh_for_buffer_region( + self, buffer_polygon, hfun_buffer, buffer_crs): - _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 + utm = utils.estimate_bounds_utm( + buffer_polygon.bounds, buffer_crs ) - idx = gdf_elems.representative_point().sindex.query( - seam.get_multipolygon(), predicate='intersects' + transformer = Transformer.from_crs(buffer_crs, utm, always_xy=True) + buffer_polygon = transform(transformer.transform, buffer_polygon) + + mesh_buf_apprx = JigsawDriver( + geom=Geom(buffer_polygon, crs=utm), + hfun=hfun_buffer, + initial_mesh=False + ).run(sieve=0) + + 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: + 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=gdf_aux_pts, opts='p' ) - flag = np.zeros(len(gdf_elems), dtype=bool) - flag[idx] = True - jig_mesh.tria3 = jig_mesh.tria3[flag] + msht_buffer.crs = utm + + utils.reproject(msht_buffer, buffer_crs) - return jig_mesh + return msht_buffer def _transform_mesh(self, mesh, out_crs): @@ -623,9 +644,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..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 @@ -144,42 +143,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/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 diff --git a/ocsmesh/utils.py b/ocsmesh/utils.py index 685fe14a..896ba59f 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 @@ -13,6 +13,7 @@ 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) @@ -23,6 +24,7 @@ LineString, LinearRing) from shapely.ops import polygonize, linemerge, unary_union import geopandas as gpd +import pandas as pd import utm @@ -94,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(): @@ -191,86 +193,18 @@ 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] + ) + poly = unary_union(elm_polys) + if isinstance(poly, Polygon): + poly = MultiPolygon([poly]) - # 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 poly def repartition_features(linestring, max_verts): @@ -1112,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: @@ -2164,8 +2098,11 @@ 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) ): mesh = jigsaw_msh_t() @@ -2179,14 +2116,282 @@ 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 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!" + ) + + mesh.value = values 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]) + + 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( + 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 shape_to_msh_t_2( + shape: Union[Polygon, MultiPolygon, gpd.GeoDataFrame, gpd.GeoSeries] +) -> jigsaw_msh_t: + gdf_shape = shape + 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 + .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_seg = ( + df_lonlat.join( + df_lonlat.groupby('index').transform(np.roll, 1, axis=0), + lsuffix='_1', + rsuffix='_2' + ).dropna() + .set_index('index') + ) + + df_nodes = ( + 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']) + ) + + # 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 = msht_from_numpy( + coordinates=df_coo[['lon', 'lat']].values, + edges=df_cnn.values, + values=np.zeros((len(df_coo) ,1)), + crs=gdf_shape.crs + ) + + 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: + """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 + + 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_2(shape) + 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 neg_shape.geoms: + holes.append(np.array(poly.representative_point().xy).ravel()) + holes = np.array(holes) + + + if aux_pts is not None: + 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 + 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)) + + input_dict = {'vertices': coords, 'segments': edges} + if len(holes): + input_dict['holes'] = holes + cdt = tr.triangulate(input_dict, opts=opts) + + msht_tri = msht_from_numpy( + coordinates=cdt['vertices'], + triangles=cdt['triangles'], + 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/pyproject.toml b/pyproject.toml index bcaafe11..e68b9bdc 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -19,13 +19,13 @@ 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.8 -> scipy dependencies = [ "colored-traceback", "fiona", "geopandas", "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"] diff --git a/tests/api/hfun.py b/tests/api/hfun.py index 3570545b..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) @@ -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/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() diff --git a/tests/api/utils.py b/tests/api/utils.py index d555b066..9955a6ee 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 @@ -67,7 +69,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) @@ -588,6 +592,44 @@ def test_crs(self): self.assertEqual(out_msht_2.crs, CRS.from_user_input('esri:102008')) + 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): + 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): @@ -653,5 +695,387 @@ def test_data_masking(self): self.assertEqual(rast.src.nodata, fill_value) +class ShapeToMeshT(unittest.TestCase): + + 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_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) + + 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_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) + + 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): + 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): + 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 TriangulatePolygon(unittest.TestCase): + + + 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_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) + + 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) + + # 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) + ) + + + def test_aux_points(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! + ) + 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)) + + + 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 GetMeshPolygon(unittest.TestCase): + def test_always_returns_multipolygon(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_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()