From 1f9a52c1863527babd6a7dda928d663f329e4b5e Mon Sep 17 00:00:00 2001 From: Keith Roberts Date: Sat, 27 May 2023 10:14:05 -0400 Subject: [PATCH] Major revision (#63) * IMPRV: adding fort.14 writer * IMPRV: adding t3s writer * IMPRV: adding in t3s writer * Forgot to commit * Bugfix: fixing reading of DEM * IMPRV: clipping DEMs and new method total_bounds for Region * Major revision of DEM class, slope edgefx, new distance from point edgefx, impr. vis. of edgefx, new helper functions in several places * Get polygon from various input types * IMPRV: adding in resolution from line function * BUGFIX: handle multipolygons. IMPRV: begin work on labeling boundaries * IMPRV: adding automatic detection of ocean boundaries * IMPRV: more refinement for the ocean boundary * IMPRV: general improvements * BUGFIX: wl mesh size function in geographic coords * Fix code style issues with Black * FIX: formatting * FIX: addressing flake8 * Fix code style issues with Black --------- Co-authored-by: oceanmesh <40673418+CHLNDDEV@users.noreply.github.com> Co-authored-by: Lint Action --- oceanmesh/__init__.py | 62 ++++++-- oceanmesh/boundary.py | 106 +++++++++++++ oceanmesh/clean.py | 57 ++++++- oceanmesh/edgefx.py | 288 +++++++++++++++++++++++++++++++----- oceanmesh/fix_mesh.py | 13 +- oceanmesh/geodata.py | 239 ++++++++++++++++++++---------- oceanmesh/grid.py | 85 +++++------ oceanmesh/mesh_generator.py | 134 ++++++++++++++--- oceanmesh/region.py | 12 ++ 9 files changed, 790 insertions(+), 206 deletions(-) create mode 100644 oceanmesh/boundary.py diff --git a/oceanmesh/__init__.py b/oceanmesh/__init__.py index 21b0f2c..5440f5e 100644 --- a/oceanmesh/__init__.py +++ b/oceanmesh/__init__.py @@ -8,31 +8,54 @@ ), "The environment variable CGAL_BIN must be set." os.add_dll_directory(os.environ["CGAL_BIN"]) -from oceanmesh.clean import (delete_boundary_faces, delete_exterior_faces, - delete_faces_connected_to_one_face, - delete_interior_faces, laplacian2, - make_mesh_boundaries_traversable) -from oceanmesh.edgefx import (bathymetric_gradient_sizing_function, - distance_sizing_function, enforce_mesh_gradation, - enforce_mesh_size_bounds_elevation, - feature_sizing_function, - multiscale_sizing_function, - wavelength_sizing_function) +from oceanmesh.boundary import identify_ocean_boundary_sections +from oceanmesh.clean import ( + delete_boundary_faces, + delete_exterior_faces, + delete_faces_connected_to_one_face, + delete_interior_faces, + laplacian2, + make_mesh_boundaries_traversable, + mesh_clean, +) +from oceanmesh.edgefx import ( + bathymetric_gradient_sizing_function, + distance_sizing_from_line_function, + distance_sizing_from_point_function, + distance_sizing_function, + enforce_mesh_gradation, + enforce_mesh_size_bounds_elevation, + feature_sizing_function, + multiscale_sizing_function, + wavelength_sizing_function, +) from oceanmesh.edges import draw_edges, get_poly_edges from oceanmesh.filterfx import filt2 -from oceanmesh.geodata import DEM, Shoreline +from oceanmesh.geodata import ( + DEM, + Shoreline, + create_circle_coords, + get_polygon_coordinates, +) from oceanmesh.grid import Grid, compute_minimum from oceanmesh.region import Region, warp_coordinates from oceanmesh.signed_distance_function import ( - Difference, Domain, Intersection, Union, create_bbox, create_circle, - multiscale_signed_distance_function, signed_distance_function) + Difference, + Domain, + Intersection, + Union, + create_bbox, + create_circle, + multiscale_signed_distance_function, + signed_distance_function, +) from .fix_mesh import fix_mesh, simp_vol - from .mesh_generator import ( generate_mesh, generate_multiscale_mesh, - plot_mesh, + plot_mesh_bathy, + plot_mesh_connectivity, write_to_fort14, write_to_t3s, ) @@ -41,16 +64,21 @@ "create_bbox", "Region", "compute_minimum", + "create_circle_coords", "bathymetric_gradient_sizing_function", "multiscale_sizing_function", "delete_boundary_faces", "delete_faces_connected_to_one_face", - "plot_mesh", + "distance_sizing_from_point_function", + "distance_sizing_from_line_function", + "plot_mesh_connectivity", + "plot_mesh_bathy", "make_mesh_boundaries_traversable", "enforce_mesh_size_bounds_elevation", "laplacian2", "delete_interior_faces", "delete_exterior_faces", + "mesh_clean", "SizeFunction", "Grid", "DEM", @@ -60,8 +88,10 @@ "Union", "Difference", "Intersection", + "identify_ocean_boundary_sections", "Shoreline", "generate_multiscale_mesh", + "get_polygon_coordinates", "distance_sizing_function", "feature_sizing_function", "enforce_mesh_gradation", diff --git a/oceanmesh/boundary.py b/oceanmesh/boundary.py new file mode 100644 index 0000000..ea3824a --- /dev/null +++ b/oceanmesh/boundary.py @@ -0,0 +1,106 @@ +import numpy as np + +from .edges import get_winded_boundary_edges +import matplotlib.pyplot as plt + +__all__ = ["identify_ocean_boundary_sections"] + + +def identify_ocean_boundary_sections( + points, + cells, + topobathymetry, + depth_threshold=-50.0, + min_nodes_threshold=10, + plot=False, +): + """Identify the contiguous sections on the ocean boundary based on depth + that could be forced in a numerical model as ocean-type boundaries (e.g., elevation-specified) + + Parameters + ---------- + points: numpy.ndarray + Array of points (x,y) + cells : numpy.ndarray + Array of cells + topobathymetry : numpy.ndarray + Array of topobathymetry values (depth below datum negative) + depth_threshold : float, optional + Depth threshold to identify ocean boundary nodes, by default -50 m below the datum + min_nodes_threshold : int, optional + Minimum number of nodes to be considered a boundary section, by default 10 + plot : bool, optional + Plot the mesh and the identified boundary sections, by default False + + Returns + -------- + boundary_sections : list + List of tuples of the nodes that define the ocean boundary sections + Note these map back into the points array. + + """ + # Identify the nodes on the boundary of the mesh + boundary_edges = get_winded_boundary_edges(cells) + boundary_edges = boundary_edges.flatten() + unique_indexes = np.unique(boundary_edges, return_index=True)[1] + boundary_nodes_unmasked = [ + boundary_edges[unique_index] for unique_index in sorted(unique_indexes) + ] + # Define a boolean array of valid nodes + bathymetry_on_boundary = topobathymetry[boundary_nodes_unmasked] + # Append a NaN value to the array to align with the original + bathymetry_on_boundary = np.append(bathymetry_on_boundary, np.nan) + stops = np.nonzero(bathymetry_on_boundary <= depth_threshold)[0] + + # Plot the mesh + if plot: + fig, ax = plt.subplots() + ax.triplot(points[:, 0], points[:, 1], cells, color="k", lw=0.1) + + first = True + boundary_sections = [] + start_node = None + end_node = None + for idx, (s1, s2) in enumerate(zip(stops[:-1], stops[1:])): + if s2 - s1 < min_nodes_threshold: + if first: + start_node = s1 + first = False + # We've reached the end of the list + elif idx == len(stops) - 2: + # Append the start and end nodes to the boundary sections list + end_node = s2 + boundary_sections.append([start_node, end_node]) + # Its not the end and we haven't found a section yet + else: + end_node = s2 + elif s2 - s1 >= min_nodes_threshold and not first: + # Append the start and end nodes to the boundary sections list + boundary_sections.append([start_node, end_node]) + # Reset the start node, the last node didn't satisfy the threshold + # and it appears we have a new section + start_node = s1 + first = True + # We've reached the end of the list + elif idx == len(stops) - 2: + # Save the end node + end_node = s2 + # Append the start and end nodes to the boundary sections list and finish + boundary_sections.append([start_node, end_node]) + if plot: + for s1, s2 in boundary_sections: + ax.scatter( + points[boundary_nodes_unmasked[s1:s2], 0], + points[boundary_nodes_unmasked[s1:s2], 1], + 5, + c="r", + ) + ax.set_title("Identified ocean boundary sections") + plt.show() + + # Map back to the original node indices associated with the points array + boundary_sections = [ + (boundary_nodes_unmasked[s1], boundary_nodes_unmasked[s2]) + for s1, s2 in boundary_sections + ] + return boundary_sections diff --git a/oceanmesh/clean.py b/oceanmesh/clean.py index 177b501..3335432 100644 --- a/oceanmesh/clean.py +++ b/oceanmesh/clean.py @@ -16,9 +16,61 @@ "delete_faces_connected_to_one_face", "delete_boundary_faces", "laplacian2", + "mesh_clean", ] +def mesh_clean( + points, + cells, + min_element_qual=0.01, + min_percent_disconnected_area=0.05, + max_iter=20, + tol=0.01, + pfix=None, +): + """Clean a mesh by removing bad quality elements and boundary faces. + + Parameters + ---------- + points : array-like + Mesh vertices. + cells : array-like + Mesh connectivity table. + min_element_qual : float, optional + Enforce Mmnimum quality of elements through deletion. The default is 0.01. + min_percent_disconnected_area : float, optional + Delete disconnected areas smaller than this threshold. The default is 0.05 x + the total area of the mesh. + max_iter : int, optional + Maximum number of iterations for the Laplacian smoothing. The default is 20. + tol : float, optional + Tolerance for the Laplacian smoothing. The default is 0.01. Laplacian terminates + after max_iter or when the maximum change in any vertex is less than tol. + pfix : array-like, optional + Indices of points to fix during the Laplacian smoothing. The default is None. + + Returns + ------- + points : array-like + Mesh vertices cleaned. + cells : array-like + Mesh connectivity table cleaned. + + """ + points, cells = make_mesh_boundaries_traversable( + points, cells, min_disconnected_area=min_percent_disconnected_area + ) + points, cells = delete_boundary_faces(points, cells, min_qual=min_element_qual) + points, cells = delete_faces_connected_to_one_face(points, cells) + points, cells = laplacian2(points, cells, max_iter=max_iter, tol=tol, pfix=pfix) + points, cells = make_mesh_boundaries_traversable( + points, cells, min_disconnected_area=min_percent_disconnected_area + ) + points, cells, _ = fix_mesh(points, cells, delete_unused=True) + return points, cells + + def _arg_sortrows(arr): """Before a multi column sort like MATLAB's sortrows""" i = arr[:, 1].argsort() # First sort doesn't need to be stable. @@ -293,7 +345,7 @@ def _depth_first_search(faces): return nflag -def delete_faces_connected_to_one_face(vertices, faces, max_iter=5): +def delete_faces_connected_to_one_face(vertices, faces, max_iter=np.inf): """Iteratively deletes faces connected to one face. Parameters @@ -304,7 +356,8 @@ def delete_faces_connected_to_one_face(vertices, faces, max_iter=5): The "uncleaned" mesh connectivity. max_iter: float, optional The number of iterations to repeatedly delete faces connected to one face - + If the mesh is well-formed, this will converge in a finite + number of iterations. Default is np.inf. Returns ------- diff --git a/oceanmesh/edgefx.py b/oceanmesh/edgefx.py index dd94353..fd70331 100644 --- a/oceanmesh/edgefx.py +++ b/oceanmesh/edgefx.py @@ -1,12 +1,17 @@ import logging +import math +import time +import geopandas as gpd +import matplotlib.pyplot as plt import numpy as np import scipy.spatial import skfmm +from _HamiltonJacobi import gradient_limit from inpoly import inpoly2 +from shapely.geometry import LineString from skimage.morphology import medial_axis -from _HamiltonJacobi import gradient_limit from oceanmesh.filterfx import filt2 from . import edges @@ -18,6 +23,8 @@ "enforce_mesh_gradation", "enforce_mesh_size_bounds_elevation", "distance_sizing_function", + "distance_sizing_from_point_function", + "distance_sizing_from_line_function", "wavelength_sizing_function", "multiscale_sizing_function", "feature_sizing_function", @@ -80,7 +87,7 @@ def enforce_mesh_size_bounds_elevation(grid, dem, bounds): return grid -def enforce_mesh_gradation(grid, gradation=0.15, crs=4326): +def enforce_mesh_gradation(grid, gradation=0.15, crs="EPSG:4326"): """Enforce a mesh size gradation bound `gradation` on a :class:`grid` Parameters @@ -127,12 +134,190 @@ def enforce_mesh_gradation(grid, gradation=0.15, crs=4326): return grid_limited +def _line_to_points_array(line): + """Convert a shapely LineString to a numpy array of points""" + return np.array(line.coords) + + +def _resample_line(row, min_edge_length): + """Resample a line to a minimum edge length""" + line = row["geometry"] + resampled_points = [] + distance = 0 + while distance < line.length: + resampled_points.append(line.interpolate(distance)) + distance += min_edge_length / 2 + resampled_line = LineString(resampled_points) + row["geometry"] = resampled_line + return row + + +def distance_sizing_from_line_function( + line_file, + bbox, + min_edge_length, + rate=0.15, + max_edge_length=None, + coarsen=1, + crs="EPSG:4326", +): + """Mesh sizes that vary linearly at `rate` from a line or lines + + Parameters + ---------- + line_file: str + Path to a vector file containing LineString(s) + bbox: list or tuple + A list or tuple of the form [xmin, xmax, ymin, ymax] denoting the bounding box of the + domain + min_edge_length: float + The minimum edge length of the mesh + rate: float + The decimal percent mesh expansion rate from the line(s) + coarsen: int + The coarsening factor of the mesh + crs: A Python int, dict, or str, optional + The coordinate reference system + max_edge_length: float, optional + The maximum edge length of the mesh + + Returns + ------- + grid: class:`Grid` + A grid ojbect with its values field populated with distance sizing + """ + logger.info("Building a distance sizing from point function...") + line_geodataframe = gpd.read_file(line_file) + assert ( + line_geodataframe.crs == crs + ), "The crs of the point geodataframe must match the crs of the grid" + # check all the geometries are points + assert all( + line_geodataframe.geometry.geom_type == "LineString" + ), "All geometries must be linestrings" + + # Resample the spacing along the lines so that the minimum edge length is met + line_geodataframe = line_geodataframe.apply( + _resample_line, axis=1, min_edge_length=min_edge_length + ) + + # Get the coordinates of the linestrings from the geodataframe + # Convert all the LineStrings in the dataframe to arrays of points + points_list = [ + _line_to_points_array(line) for line in line_geodataframe["geometry"] + ] + points = np.concatenate(points_list) + + # Create a mesh size function grid + grid = Grid( + bbox=bbox, + dx=min_edge_length * coarsen, + hmin=min_edge_length, + extrapolate=True, + values=0.0, + crs=crs, + ) + # create phi (-1 where point(s) intersect grid points -1 elsewhere 0) + phi = np.ones(shape=(grid.nx, grid.ny)) + lon, lat = grid.create_grid() + # find location of points on grid + indices = grid.find_indices(points, lon, lat) + phi[indices] = -1.0 + try: + dis = np.abs(skfmm.distance(phi, [grid.dx, grid.dy])) + except ValueError: + logger.info("0-level set not found in domain or grid malformed") + dis = np.zeros((grid.nx, grid.ny)) + 999 + tmp = min_edge_length + dis * rate + if max_edge_length is not None: + tmp[tmp > max_edge_length] = max_edge_length + grid.values = np.ma.array(tmp) + grid.build_interpolant() + return grid + + +def distance_sizing_from_point_function( + point_file, + bbox, + min_edge_length, + rate=0.15, + max_edge_length=None, + coarsen=1, + crs="EPSG:4326", +): + '''Mesh sizes that vary linearly at `rate` from a point or points + contained within a geopandas dataframe. + + Parameters + ---------- + point_geodataframe: str + Path to a vector file containing Points + bbox: list or tuple + A list or tuple of the form [xmin, xmax, ymin, ymax] denoting the bounding box of the + domain + min_edge_length: float + The minimum edge length of the mesh + rate: float + The decimal percent mesh expansion rate from the point(s) + coarsen: int + The coarsening factor of the background grid + crs: A Python int, dict, or str, optional + The coordinate reference system + + Returns + ------- + grid: class:`Grid` + A grid ojbect with its values field gradient limited + + """ + + + ''' + logger.info("Building a distance sizing from point function...") + point_geodataframe = gpd.read_file(point_file) + assert ( + point_geodataframe.crs == crs + ), "The crs of the point geodataframe must match the crs of the grid" + # check all the geometries are points + assert all( + point_geodataframe.geometry.geom_type == "Point" + ), "All geometries must be points" + # Get the coordinates of the points from the geodataframe + points = np.array(point_geodataframe.geometry.apply(lambda x: (x.x, x.y)).tolist()) + # Create a mesh size function grid + grid = Grid( + bbox=bbox, + dx=min_edge_length * coarsen, + hmin=min_edge_length, + extrapolate=True, + values=0.0, + crs=crs, + ) + # create phi (-1 where point(s) intersect grid points -1 elsewhere 0) + phi = np.ones(shape=(grid.nx, grid.ny)) + lon, lat = grid.create_grid() + # find location of points on grid + indices = grid.find_indices(points, lon, lat) + phi[indices] = -1.0 + try: + dis = np.abs(skfmm.distance(phi, [grid.dx, grid.dy])) + except ValueError: + logger.info("0-level set not found in domain or grid malformed") + dis = np.zeros((grid.nx, grid.ny)) + 999 + tmp = min_edge_length + dis * rate + if max_edge_length is not None: + tmp[tmp > max_edge_length] = max_edge_length + grid.values = np.ma.array(tmp) + grid.build_interpolant() + return grid + + def distance_sizing_function( shoreline, rate=0.15, max_edge_length=None, - coarsen=1, - crs=4326, + coarsen=1.0, + crs="EPSG:4326", ): """Mesh sizes that vary linearly at `rate` from coordinates in `obj`:Shoreline Parameters @@ -190,12 +375,11 @@ def distance_sizing_function( try: dis = np.abs(skfmm.distance(phi, [grid.dx, grid.dy])) except ValueError: - logger.info("0-level set not found in domain") + logger.info("0-level set not found in domain or grid malformed") dis = np.zeros((grid.nx, grid.ny)) + 999 tmp = shoreline.h0 + dis * rate if max_edge_length is not None: tmp[tmp > max_edge_length] = max_edge_length - grid.values = np.ma.array(tmp, mask=mask) grid.build_interpolant() return grid @@ -207,10 +391,10 @@ def bathymetric_gradient_sizing_function( filter_quotient=50, min_edge_length=None, max_edge_length=None, - min_elevation_cutoff=50, - type_of_filter="barotropic", - filter_cutoffs=None, - crs=4326, + min_elevation_cutoff=-50.0, + type_of_filter="lowpass", + filter_cutoffs=1000, + crs="EPSG:4326", ): """Mesh sizes that vary proportional to the bathymetryic gradient. Bathymetry is filtered by default using a fraction of the @@ -226,9 +410,9 @@ def bathymetric_gradient_sizing_function( slope_parameter: integer, optional The number of nodes to resolve bathymetryic gradients min_edge_length: float, optional - The minimum allowable edge length in meters in the domain. + The minimum allowable edge length in CRS units in the domain. max_edge_length: float, optional - The maximum allowable edge length in meters in the domain. + The maximum allowable edge length in CRS units in the domain. min_elevation_cutoff: float, optional abs(elevation) < this value the sizing function is not calculated. type_of_filter: str, optional @@ -250,40 +434,60 @@ def bathymetric_gradient_sizing_function( logger.info("Building a slope length sizing function...") - x, y = dem.create_grid() - tmpz = dem.eval((x, y)) + xg, yg = dem.create_grid() + tmpz = dem.eval((xg, yg)) + grid = Grid( bbox=dem.bbox, - dx=dem.dx, - dy=dem.dy, + dx=min_edge_length / 2.0, + dy=min_edge_length / 2.0, extrapolate=True, values=0.0, crs=crs, - hmin=dem.dx, + hmin=min_edge_length, ) logger.info(f"Enforcing a minimum elevation cutoff of {min_elevation_cutoff}") - tmpz[np.abs(tmpz) < min_elevation_cutoff] = min_elevation_cutoff + tmpz[tmpz >= min_elevation_cutoff] = min_elevation_cutoff dx, dy = dem.dx, dem.dy # for gradient function nx, ny = dem.nx, dem.ny - grid_details = (nx, ny, dx, dy) - xg, yg = dem.create_grid() coords = (xg, yg) + if crs == "EPSG:4326": + mean_latitude = np.mean(dem.bbox[2:]) + meters_per_degree = ( + 111132.92 + - 559.82 * np.cos(2 * mean_latitude) + + 1.175 * np.cos(4 * mean_latitude) + - 0.0023 * np.cos(6 * mean_latitude) + ) + dy *= meters_per_degree + dx *= meters_per_degree + + grid_details = (nx, ny, dx, dy) + if type_of_filter == "barotropic" and filter_quotient > 0: + # Temporary fix for barotropic filter not implemented correctly yet + assert ValueError("Barotropic filter not implemented yet") + logger.info("Baroptropic Rossby radius calculation...") bs, time_taken = rossby_radius_filter( tmpz, dem.bbox, grid_details, coords, filter_quotient, True ) elif type_of_filter == "baroclinic" and filter_quotient > 0: + # Temporary fix for barotropic filter not implemented correctly yet + assert ValueError("Baroclinic filter not implemented yet") + logger.info("Baroclinic Rossby radius calculation...") bs, time_taken = rossby_radius_filter( tmpz, dem.bbox, grid_details, coords, filter_quotient, False ) elif "pass" in type_of_filter: - logger.info("Using a {type_of_filter} filter...") - bs = filt2(tmpz, dy, filter_cutoffs, type_of_filter) + logger.info(f"Using a {type_of_filter} filter...") + tmpzs = filt2(tmpz, dy, filter_cutoffs, type_of_filter) + by, bx = _earth_gradient(tmpzs, dy, dx) + bs = np.sqrt(bx**2 + by**2) # get overall slope else: msg = f"The type_of_filter {type_of_filter} is not known and remains off" logger.info(msg) @@ -292,8 +496,12 @@ def bathymetric_gradient_sizing_function( # Calculating the slope function eps = 1e-10 # small number to approximate derivative - dp = np.maximum(1, tmpz) - grid.values = (2 * np.pi / slope_parameter) * dp / (bs + eps) + dp = np.clip(tmpz, None, -1) + grid.values = (2 * np.pi / slope_parameter) * np.abs(dp) / (bs + eps) + + # Convert back to degrees from meters (if geographic) + if crs == "EPSG:4326": + grid.values /= meters_per_degree if max_edge_length is not None: grid.values[grid.values > max_edge_length] = max_edge_length @@ -339,8 +547,6 @@ def rossby_radius_filter(tmpz, bbox, grid_details, coords, rbfilt, barot): the time taken to prform the filtering process. """ - import math - import time x0, xN, y0, yN = bbox @@ -462,7 +668,7 @@ def feature_sizing_function( min_edge_length=None, max_edge_length=None, plot=False, - crs=4326, + crs="EPSG:4326", ): """Mesh sizes vary proportional to the width or "thickness" of the shoreline @@ -484,7 +690,6 @@ def feature_sizing_function( crs: A Python int, dict, or str, optional The coordinate reference system - Returns ------- :class:`Grid` object @@ -534,8 +739,6 @@ def feature_sizing_function( dis = np.abs(skfmm.distance(phi2, [grid_calc.dx, grid_calc.dy])) if plot: - import matplotlib.pyplot as plt - fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4)) ax1.pcolor(lon, lat, skel, cmap=plt.cm.gray) ax1.axis("off") @@ -579,7 +782,7 @@ def wavelength_sizing_function( max_edge_length=None, period=12.42 * 3600, # M2 period in seconds gravity=9.81, # m/s^2 - crs=4326, + crs="EPSG:4326", ): """Mesh sizes that vary proportional to an estimate of the wavelength of a period (default M2-period) @@ -615,15 +818,32 @@ def wavelength_sizing_function( lon, lat = dem.create_grid() tmpz = dem.eval((lon, lat)) - grav = 9.807 + if crs == "EPSG:4326": + mean_latitude = np.mean(dem.bbox[2:]) + meters_per_degree = ( + 111132.92 + - 559.82 * np.cos(2 * mean_latitude) + + 1.175 * np.cos(4 * mean_latitude) + - 0.0023 * np.cos(6 * mean_latitude) + ) + grid = Grid( bbox=dem.bbox, dx=dem.dx, dy=dem.dy, extrapolate=True, values=0.0, crs=crs ) tmpz[np.abs(tmpz) < 1] = 1 - grid.values = period * np.sqrt(grav * np.abs(tmpz)) / wl + grid.values = period * np.sqrt(gravity * np.abs(tmpz)) / wl + + # Convert back to degrees from meters (if geographic) + if crs == "EPSG:4326": + grid.values /= meters_per_degree + if min_edgelength is None: min_edgelength = np.amin(grid.values) + else: + grid.values[grid.values < min_edgelength] = min_edgelength + grid.hmin = min_edgelength + if max_edge_length is not None: grid.values[grid.values > max_edge_length] = max_edge_length @@ -712,7 +932,7 @@ def func(qpts): return func, new_list_of_grids -def _earth_gradient(F, dy, dx): +def _earth_gradient(F, dx, dy): """ earth_gradient(F,HX,HY), where F is 2-D, uses the spacing specified by HX and HY. HX and HY can either be scalars to specify diff --git a/oceanmesh/fix_mesh.py b/oceanmesh/fix_mesh.py index 6b3c36d..525b677 100644 --- a/oceanmesh/fix_mesh.py +++ b/oceanmesh/fix_mesh.py @@ -1,4 +1,5 @@ import numpy as np +import warnings def simp_qual(p, t): @@ -18,9 +19,15 @@ def length(p1): a = length(p[t[:, 1]] - p[t[:, 0]]) b = length(p[t[:, 2]] - p[t[:, 0]]) c = length(p[t[:, 2]] - p[t[:, 1]]) - r = 0.5 * np.sqrt((b + c - a) * (c + a - b) * (a + b - c) / (a + b + c)) - R = a * b * c / np.sqrt((a + b + c) * (b + c - a) * (c + a - b) * (a + b - c)) - return 2 * r / R + # Suppress Runtime warnings here because we know that mult1/denom1 can be negative + # as the mesh is being cleaned + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + mult1 = (b + c - a) * (c + a - b) * (a + b - c) / (a + b + c) + denom1 = np.sqrt((a + b + c) * (b + c - a) * (c + a - b) * (a + b - c)) + r = 0.5 * mult1 + R = a * b * c / denom1 + return 2 * r / R def fix_mesh(p, t, ptol=2e-13, dim=2, delete_unused=False): diff --git a/oceanmesh/geodata.py b/oceanmesh/geodata.py index a56e0ec..27d74ce 100644 --- a/oceanmesh/geodata.py +++ b/oceanmesh/geodata.py @@ -1,16 +1,21 @@ import errno import logging import os +from pathlib import Path import fiona import geopandas as gpd import matplotlib.path as mpltPath +import matplotlib.pyplot as plt import numpy as np import numpy.linalg -import rioxarray as rxr +import rasterio +import rasterio.crs +import rasterio.warp import shapely.geometry import shapely.validation from pyproj import CRS +from rasterio.windows import from_bounds from .grid import Grid from .region import Region @@ -20,7 +25,48 @@ logger = logging.getLogger(__name__) -__all__ = ["Shoreline", "DEM"] +__all__ = ["Shoreline", "DEM", "get_polygon_coordinates", "create_circle_coords"] + + +def create_circle_coords(radius, center, arc_res): + """ + Given a radius and a center point, creates a numpy array of coordinates + defining a circle in a CCW direction with a given arc resolution. + + Parameters: + radius (float): the radius of the circle + center (tuple): the (x,y) coordinates of the center point + arc_res (float): the arc resolution of the circle in degrees + + Returns: + numpy.ndarray: an array of (x,y) coordinates defining the circle + """ + # Define the angle array with the given arc resolution + angles = np.arange(0, 360 + arc_res, arc_res) * np.pi / 180 + + # Calculate the (x,y) coordinates of the circle points + x_coords = center[0] + radius * np.cos(angles) + y_coords = center[1] + radius * np.sin(angles) + + # Combine the (x,y) coordinates into a single array + coords = np.column_stack((x_coords, y_coords)) + + return coords + + +def get_polygon_coordinates(vector_file): + """Get the coordinates of a polygon from a vector file or plain csv file""" + # detect if file is a shapefile or a geojson or geopackage + if ( + vector_file.endswith(".shp") + or vector_file.endswith(".geojson") + or vector_file.endswith(".gpkg") + ): + gdf = gpd.read_file(vector_file) + polygon = np.array(gdf.iloc[0].geometry.exterior.coords.xy).T + elif vector_file.endswith(".csv"): + polygon = np.loadtxt(vector_file, delimiter=",") + return polygon def _convert_to_array(lst): @@ -435,7 +481,29 @@ class Shoreline(Region): """ The shoreline class extends :class:`Region` to store data that is later used to create signed distance functions to - represent irregular shoreline geometries. + represent irregular shoreline geometries. This data + is also involved in developing mesh sizing functions. + + Parameters + ---------- + shp : str or pathlib.Path + Path to shapefile containing shoreline data. + bbox : tuple + Bounding box of the region of interest. The format is + (xmin, xmax, ymin, ymax). + h0 : float + Minimum grid spacing. + crs : str, optional + Coordinate reference system of the shapefile. Default is + 'EPSG:4326'. + refinements : int, optional + Number of refinements to apply to the shoreline. Default is 1. + minimum_area_mult : float, optional + Minimum area multiplier. Default is 4.0. + Note that features with area less than h0*minimum_area_mult + are removed. + smooth_shoreline : bool, optional + Smooth the shoreline. Default is True. """ def __init__( @@ -443,11 +511,14 @@ def __init__( shp, bbox, h0, - crs=4326, + crs="EPSG:4326", refinements=1, minimum_area_mult=4.0, smooth_shoreline=True, ): + if isinstance(shp, str): + shp = Path(shp) + if isinstance(bbox, tuple): _boubox = np.asarray(_create_boubox(bbox)) else: @@ -559,6 +630,9 @@ def _read(self): # transform if necessary s = self.transform_to(gpd.read_file(self.shp), self.crs) + # Explode to remove multipolygons or multi-linestrings (if present) + s = s.explode(index_parts=True) + polys = [] # store polygons delimiter = np.empty((1, 2)) @@ -569,10 +643,13 @@ def _read(self): # extent of geometry bbox2 = [g.bounds[r] for r in re] if _is_overlapping(_bbox, bbox2): - if g.type == "LineString": + if g.geom_type == "LineString": poly = np.asarray(g.coords) - else: # a polygon + elif g.geom_type == "Polygon": # a polygon poly = np.asarray(g.exterior.coords.xy).T + else: + raise ValueError(f"Unsupported geometry type: {g.geom_type}") + poly = remove_dup(poly) polys.append(np.row_stack((poly, delimiter))) @@ -595,8 +672,6 @@ def plot( ylim=None, ): """Visualize the content in the shp field of Shoreline""" - import matplotlib.pyplot as plt - flg1, flg2 = False, False if ax is None: @@ -654,84 +729,90 @@ def plot( class DEM(Grid): """ Digitial elevation model read in from a tif or NetCDF file - parent class is a :class:`Grid` if dem input is a string. If dem is an array, - it assumes a uniform grid-spacing, while if a function is given, the grid will consist of - 1001 points in both horizontal directions. """ - def __init__(self, dem, crs=4326, bbox=None, nnodes=1000): - if type(dem) == str: - basename, ext = os.path.splitext(dem) - if ext.lower() in [".nc"] or [".tif"]: - topobathy, reso, bbox = self._read(dem, bbox, crs) - topobathy = topobathy.astype(float) - - else: - raise ValueError(f"DEM file {dem} has unknown format {ext[1:]}.") + def __init__(self, dem, crs="EPSG:4326", bbox=None, extrapolate=False): + """Read in a DEM from a tif or NetCDF file for later use + in developing mesh sizing functions. + + Parameters + ---------- + dem : str or pathlib.Path + Path to the DEM file + crs : str, optional + Coordinate reference system of the DEM, by default 'EPSG:4326' + bbox : oceanmesh.Region class + Bounding box of the DEM, by default None. + Note that if none, it will read in the entire DEM. + extrapolate : bool, optional + Extrapolate the DEM outside the bounding box, by default False + """ - self.dem = dem + if isinstance(dem, str): + dem = Path(dem) + + if bbox is not None: + assert isinstance(bbox, Region), "bbox must be a Region class object" + # Extract the total bounds from the extent + bbox = bbox.total_bounds + + if dem.exists(): + msg = f"Reading in {dem}" + logger.info(msg) + # Open the raster file using rasterio + with rasterio.open(dem) as src: + nodata_value = src.nodata + self.meta = src.meta + # entire DEM is read in + if bbox is None: + bbox = src.bounds + topobathy = src.read(1) + # then clip the DEM to the box + else: + # + _bbox = (bbox[0], bbox[2], bbox[1], bbox[3]) + window = from_bounds(*_bbox, transform=src.transform) + topobathy = src.read(1, window=window, masked=True) + topobathy = np.transpose(topobathy, (1, 0)) + # Ensure its a floating point array + topobathy = topobathy.astype(np.float64) + topobathy[ + topobathy == nodata_value + ] = np.nan # set the no-data value to nan + elif not dem.exists(): + raise FileNotFoundError(f"File {dem} could not be located.") - elif type(dem) == numpy.ndarray: - topobathy = dem.astype(float) - reso = ( - (bbox[1] - bbox[0]) / topobathy.shape[0], - (bbox[3] - bbox[2]) / topobathy.shape[1], - ) - self.dem = "input" - - elif callable(dem): # if input is a function - dx = (bbox[1] - bbox[0]) / nnodes - lon, lat = ( - np.arange(bbox[0], bbox[1] + dx, dx), - np.arange(bbox[2], bbox[3] + dx, dx), - ) - reso = (dx, dx) - lon, lat = np.meshgrid(lon, lat) - topobathy = np.rot90(dem(lon, lat), 2) - self.dem = "function" - - topobathy[abs(topobathy) > 1e5] = np.NaN super().__init__( bbox=bbox, crs=crs, - dx=abs(reso[0]), - dy=abs(reso[1]), - values=np.rot90(topobathy, 3), + dx=self.meta["transform"][0], + dy=abs( + self.meta["transform"][4] + ), # Note: grid spacing in y-direction is negative. + values=topobathy, + extrapolate=extrapolate, # user-specified potentially "dangerous" option ) super().build_interpolant() - @staticmethod - def transform_to(xry, dst_crs): - """Transform xarray ``xry`` representing - a raster to dst_crs - """ - dst_crs = CRS.from_user_input(dst_crs) - if not xry.crs.equals(dst_crs): - msg = f"Reprojecting raster from {xry.crs} to {dst_crs}" - logger.debug(msg) - xry = xry.rio.reproject(dst_crs) - return xry - - def _read(self, filename, bbox, crs): - """Read in a digitial elevation model from a NetCDF or GeoTif file""" - logger.debug("Entering: DEM._read") - - msg = f"Reading in {filename}" - logger.info(msg) - - with rxr.open_rasterio(filename) as r: - # warp/reproject it if necessary - r = self.transform_to(r, crs) - # entire DEM is read in - if bbox is None: - bnds = r.rio.bounds() - bbox = (bnds[0], bnds[2], bnds[1], bnds[3]) - else: - # then we clip the DEM to the box - r = r.rio.clip_box( - minx=bbox[0], miny=bbox[2], maxx=bbox[1], maxy=bbox[3] - ) - topobathy = r.data[0] - - logger.debug("Exiting: DEM._read") - return topobathy, r.rio.resolution(), bbox + def flip(self): + """Flip the DEM upside down""" + self.values = -self.values + super().build_interpolant() + return self + + def plot(self, coarsen=1, holding=False, **kwargs): + """Visualize the DEM""" + fig, ax, pc = super().plot( + coarsen=coarsen, + holding=True, + cmap="terrain", + **kwargs, + ) + ax.set_xlabel("X-coordinate") + ax.set_ylabel("Y-coordinate") + ax.set_aspect("equal") + cbar = fig.colorbar(pc) + cbar.set_label("Topobathymetric depth (m)") + if not holding: + plt.show() + return fig, ax diff --git a/oceanmesh/grid.py b/oceanmesh/grid.py index 4a237f0..670103b 100644 --- a/oceanmesh/grid.py +++ b/oceanmesh/grid.py @@ -80,13 +80,20 @@ class Grid(Region): """ def __init__( - self, bbox, dx, dy=None, crs=4326, hmin=None, values=None, extrapolate=False + self, + bbox, + dx, + dy=None, + crs="EPSG:4326", + hmin=None, + values=None, + extrapolate=False, ): super().__init__(bbox, crs) if dy is None: - dy = dx + dy = dx # equidistant grid in both x and y dirs if not passed self.bbox = bbox - self.x0y0 = (bbox[0], bbox[2]) # bottom left corner coordinates + self.x0y0 = (bbox[0], bbox[2]) # bottom left corner coordinates (x,y) self.dx = dx self.dy = dy self.nx = None # int((self.bbox[1] - self.bbox[0]) // self.dx) + 1 @@ -102,8 +109,8 @@ def dx(self): @dx.setter def dx(self, value): - if value < 0: - raise ValueError("Grid spacing (dx) must be > 0.0") + if value <= 0: + raise ValueError("Grid spacing (dx) must be >= 0.0") self.__dx = value @property @@ -112,8 +119,8 @@ def dy(self): @dy.setter def dy(self, value): - if value < 0: - raise ValueError("Grid spacing (dy) must be > 0.0") + if value <= 0: + raise ValueError("Grid spacing (dy) must be >= 0.0") self.__dy = value @property @@ -124,7 +131,7 @@ def values(self): def values(self, data): if np.isscalar(data): self.nx = int((self.bbox[1] - self.bbox[0]) // self.dx) + 1 - self.ny = int((self.bbox[3] - self.bbox[2]) // self.dy) + 1 + self.ny = abs(int((self.bbox[3] - self.bbox[2]) // self.dy) + 1) data = np.tile(data, (self.nx, self.ny)) self.__values = data self.nx, self.ny = data.shape @@ -151,8 +158,9 @@ def create_vectors(self): 1D array contain data with `float` type of y-coordinates. """ - x = self.x0y0[0] + np.arange(0, self.nx) * self.dx - y = self.x0y0[1] + np.arange(0, self.ny) * self.dy + x = self.x0y0[0] + np.arange(0, self.nx) * self.dx # ascending monotonically + y = self.x0y0[1] + np.arange(0, self.ny) * abs(self.dy) + y = y[::-1] # descending monotonically return x, y def create_grid(self): @@ -171,7 +179,8 @@ def create_grid(self): """ x, y = self.create_vectors() - return np.meshgrid(x, y, sparse=False, indexing="ij") + xg, yg = np.meshgrid(x, y, indexing="ij") + return xg, yg def find_indices(self, points, lon, lat, tree=None, k=1): """Find linear indices `indices` into a 2D array such that they @@ -334,63 +343,39 @@ def blend_into(self, coarse, blend_width=10, p=1, nnear=6, eps=0.0): def plot( self, - hold=False, - show=True, - vmin=None, - vmax=None, + holding=False, coarsen=1, - xlabel=None, - ylabel=None, - title=None, - cbarlabel=None, - filename=None, - xlim=None, - ylim=None, + plot_colorbar=False, + **kwargs, ): """Visualize the values in :obj:`Grid` Parameters ---------- - hold: boolean, optional + holding: boolean, optional Whether to create a new plot axis. Returns ------- + fig: ax: handle to axis of plot handle to axis of plot. """ - - x, y = self.create_grid() - + _xg, _yg = self.create_grid() fig, ax = plt.subplots() ax.axis("equal") - c = ax.pcolormesh( - x[:-1:coarsen, :-1:coarsen], - y[:-1:coarsen, :-1:coarsen], - self.values[:-1:coarsen, :-1:coarsen], - vmin=vmin, - vmax=vmax, - shading="auto", + pc = ax.pcolor( + _xg[::coarsen, ::coarsen], + _yg[::coarsen, ::coarsen], + self.values[::coarsen, ::coarsen], + **kwargs, ) - cbar = plt.colorbar(c) - if cbarlabel is not None: - cbar.set_label(cbarlabel) - if xlabel is not None: - ax.set_xlabel(xlabel) - if ylabel is not None: - ax.set_ylabel(ylabel) - if title is not None: - ax.set_title(title) - if xlim is not None: - ax.set_xlim(xlim) - if ylim is not None: - ax.set_ylim(ylim) - if hold is False and show: + if plot_colorbar: + fig.colorbar(pc) + if holding is False: plt.show() - if filename is not None: - plt.savefig(filename) - return ax + return fig, ax, pc def build_interpolant(self): """Construct a RegularGriddedInterpolant sizing function stores it as diff --git a/oceanmesh/mesh_generator.py b/oceanmesh/mesh_generator.py index f6695ed..ceba4d6 100644 --- a/oceanmesh/mesh_generator.py +++ b/oceanmesh/mesh_generator.py @@ -14,31 +14,43 @@ from .edgefx import multiscale_sizing_function from .fix_mesh import fix_mesh from .grid import Grid -from .signed_distance_function import (Domain, - multiscale_signed_distance_function) +from .signed_distance_function import Domain, multiscale_signed_distance_function logger = logging.getLogger(__name__) __all__ = [ "generate_mesh", "generate_multiscale_mesh", - "plot_mesh", + "plot_mesh_connectivity", + "plot_mesh_bathy", "write_to_fort14", "write_to_t3s", ] -def write_to_fort14(points, cells, filepath, project_name="Created with pyoceanmesh"): +def write_to_fort14( + points, + cells, + filepath, + topobathymetry=None, + project_name="Created with oceanmesh", + flip_bathymetry=False, +): """ - Write mesh data to a fort.14 file. - - Parameters: + Parameters + ----------- points (numpy.ndarray): An array of shape (np, 2) containing the x, y coordinates of the mesh nodes. cells (numpy.ndarray): An array of shape (ne, 3) containing the indices of the nodes that form each mesh element. filepath (str): The file path to write the fort.14 file to. + topobathymetry (numpy.ndarray): An array of shape (np, 1) containing the topobathymetry values at each node. + project_name (str): The name of the project to be written to the fort.14 file. + flip_bathymetry (bool): If True, the bathymetry values will be multiplied by -1. Returns: - Message indicating file written. + -------- + points (numpy.ndarray): An array of shape (np, 2) containing the x, y coordinates of the mesh nodes. + cells (numpy.ndarray): An array of shape (ne, 3) containing the indices of the nodes that form each mesh element. + filepath (str): The file path to write the fort.14 file to. """ logger.info("Exporting mesh to fort.14 file...") @@ -46,13 +58,26 @@ def write_to_fort14(points, cells, filepath, project_name="Created with pyoceanm npoints = np.size(points, 0) nelements = np.size(cells, 0) + if topobathymetry is not None: + assert ( + len(topobathymetry) == npoints + ), "topobathymetry must be the same length as points" + else: + topobathymetry = np.zeros((npoints, 1)) + + if flip_bathymetry: + topobathymetry *= -1 + # Shift cell indices by 1 (fort.14 uses 1-based indexing) cells += 1 # Open file for writing with open(filepath, "w") as f_id: # Write mesh name - f_id.write(f"{project_name} \n") + if flip_bathymetry: + f_id.write(f"{project_name} (bathymetry flipped) \n") + else: + f_id.write(f"{project_name} \n") # Write number of nodes and elements np.savetxt( @@ -67,7 +92,7 @@ def write_to_fort14(points, cells, filepath, project_name="Created with pyoceanm for k in range(npoints): np.savetxt( f_id, - np.column_stack((k + 1, points[k][0], points[k][1], 0.0)), + np.column_stack((k + 1, points[k][0], points[k][1], topobathymetry[k])), delimiter=" ", fmt="%i %f %f %f", newline="\n", @@ -89,6 +114,7 @@ def write_to_fort14(points, cells, filepath, project_name="Created with pyoceanm return f"Wrote the mesh to {filepath}..." + def write_to_t3s(points, cells, filepath): """ Write mesh data to a t3s file. @@ -167,16 +193,79 @@ def write_to_t3s(points, cells, filepath): return f"Wrote the mesh to {filepath}..." +def plot_mesh_connectivity(points, cells, show_plot=True): + """Plot the mesh connectivity using matplotlib's triplot function. + Parameters + ---------- + + points : numpy.ndarray + A 2D array containing the x and y coordinates of the points. + cells : numpy.ndarray + A 2D array containing the connectivity information for the triangles. + show_plot : bool, optional + Whether to show the plot or not. The default is True. -def plot_mesh(points, cells, count=0, show=True, pause=999): + Returns + ------- + ax : matplotlib.axes.Axes + The axes object containing the plot. + """ triang = tri.Triangulation(points[:, 0], points[:, 1], cells) - plt.triplot(triang) - plt.gca().set_aspect("equal", adjustable="box") - plt.title(f"Iteration {count}") - if show: + fig, ax = plt.subplots(figsize=(10, 10)) + ax.triplot(triang) + ax.set_aspect("equal", adjustable="box") + ax.set_title("Mesh connectivity") + if show_plot: plt.show(block=False) - plt.pause(pause) - plt.close() + return ax + + +def plot_mesh_bathy(points, bathymetry, connectivity, show_plot=True): + """ + Create a tricontourf plot of the bathymetry data associated with the points, + using the triangle connectivity information to plot the contours. + + Parameters + ---------- + points : numpy.ndarray + A 2D array containing the x and y coordinates of the points. + bathymetry : numpy.ndarray + A 1D array containing the bathymetry values associated with each point. + connectivity : numpy.ndarray + A 2D array containing the connectivity information for the triangles. + show_plot : bool, optional + Whether or not to display the plot. Default is True. + + Returns + ------- + matplotlib.axes._subplots.AxesSubplot + The axis handle of the plot. + + """ + # Create a Triangulation object using the points and connectivity table + triangulation = tri.Triangulation(points[:, 0], points[:, 1], connectivity) + + # Create a figure and axis object + fig, ax = plt.subplots(figsize=(10, 10)) + + # Plot the tricontourf + tricontourf = ax.tricontourf(triangulation, bathymetry, cmap="jet") + + # Add colorbar + plt.colorbar(tricontourf) + + # Set axis labels + ax.set_xlabel("Longitude") + ax.set_ylabel("Latitude") + + # Set title + ax.set_title("Mesh Topobathymetry") + + # Show the plot if requested + if show_plot: + plt.show() + + return ax def _parse_kwargs(kwargs): @@ -197,6 +286,7 @@ def _parse_kwargs(kwargs): "blend_max_iter", "blend_nnear", "lock_boundary", + "pseudo_dt", }: pass else: @@ -317,6 +407,8 @@ def generate_mesh(domain, edge_length, **kwargs): The minimum element size in the domain. REQUIRED IF NOT USING :class:`edge_length` * *plot* (``int``) -- The mesh is visualized every `plot` meshing iterations. + * *pseudo_dt* (``float``) -- + The pseudo time step for the meshing algorithm. (default==0.2) Returns ------- @@ -335,6 +427,7 @@ def generate_mesh(domain, edge_length, **kwargs): "min_edge_length": None, "plot": 999999, "lock_boundary": False, + "pseudo_dt": 0.2, } opts.update(kwargs) _parse_kwargs(kwargs) @@ -353,7 +446,7 @@ def generate_mesh(domain, edge_length, **kwargs): np.random.seed(opts["seed"]) L0mult = 1 + 0.4 / 2 ** (_DIM - 1) - delta_t = 0.20 + delta_t = opts["pseudo_dt"] geps = 1e-3 * np.amin(min_edge_length) deps = np.sqrt(np.finfo(np.double).eps) # * np.amin(min_edge_length) @@ -395,6 +488,7 @@ def generate_mesh(domain, edge_length, **kwargs): _, bpts = _external_topology(p, t) for fix in bpts: ifix.append(_closest_node(fix, p)) + nfix = len(ifix) # Find where pfix went if nfix > 0: @@ -406,10 +500,6 @@ def generate_mesh(domain, edge_length, **kwargs): # Remove points outside the domain t = _remove_triangles_outside(p, t, fd, geps) - # plot the mesh every count iterations - if count % opts["plot"] == 0 and count != 0: - plot_mesh(p, t, count=count, pause=5.0) - # Number of iterations reached, stop. if count == (max_iter - 1): p, t, _ = fix_mesh(p, t, dim=_DIM, delete_unused=True) diff --git a/oceanmesh/region.py b/oceanmesh/region.py index 15a47bb..668dd8b 100644 --- a/oceanmesh/region.py +++ b/oceanmesh/region.py @@ -25,6 +25,18 @@ def crs(self): def bbox(self): return self.__bbox + @property + def total_bounds(self): + if isinstance(self.bbox, tuple): + return self.bbox + else: + return ( + self.bbox[:, 0].min(), + self.bbox[:, 0].max(), + self.bbox[:, 1].min(), + self.bbox[:, 1].max(), + ) + @bbox.setter def bbox(self, value): if isinstance(value, tuple):