-
Notifications
You must be signed in to change notification settings - Fork 4
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
prevent high orthogonality after coastline cutting #150
Labels
enhancement
New feature or request
Comments
88 tasks
Simplified testcase: import matplotlib.pyplot as plt
plt.close('all')
# import dfm_tools as dfmt
# import contextily as ctx
import numpy as np
from meshkernel import (ProjectionType, MakeGridParameters, MeshKernel,
MeshRefinementParameters, GriddedSamples, RefinementType,
GeometryList, DeleteMeshOption,
)
# input params
lon_min, lon_max, lat_min, lat_max = 147.75, 147.9, -40.4, -40.25
dxy = 0.05 #0.05
crs = 'EPSG:4326'
# grid generation and refinement with GEBCO bathymetry
# create base grid
projection = ProjectionType.SPHERICAL
make_grid_parameters = MakeGridParameters(angle=0,
origin_x=lon_min,
origin_y=lat_min,
upper_right_x=lon_max,
upper_right_y=lat_max,
block_size_x=dxy,
block_size_y=dxy)
mk = MeshKernel(projection=projection)
mk.curvilinear_compute_rectangular_grid_on_extension(make_grid_parameters)
mk.curvilinear_convert_to_mesh2d() #convert to ugrid/mesh2d
#refine
min_edge_size = 500 #in meters
print('Refining basegrid')
gebco_elev = np.array([[-37, -37, -36, -36],
[-41, -39, -38, -34],
[-42, -32, -27, -6],
[-38, -6, -2, 2],
[-43, -42, -31, -20]],dtype=np.int16)
lat_np = np.array([-40.397917, -40.364583, -40.33125 , -40.297917, -40.264583])
lon_np = np.array([147.76875 , 147.802083, 147.835417, 147.86875 ])
values_np = gebco_elev.flatten()
gridded_samples = GriddedSamples(x_coordinates=lon_np,y_coordinates=lat_np,values=values_np)
#refinement
mesh_refinement_parameters = MeshRefinementParameters(min_edge_size=min_edge_size, #always in meters
refinement_type=RefinementType(1), #Wavecourant/1,
connect_hanging_nodes=True, #set to False to do multiple refinement steps (e.g. for multiple regions)
smoothing_iterations=2,
max_courant_time=120)
mk.mesh2d_refine_based_on_gridded_samples(gridded_samples=gridded_samples,
mesh_refinement_params=mesh_refinement_parameters,
use_nodal_refinement=True) #TODO: what does this do?
print("max ortho after refining", mk.mesh2d_get_orthogonality().values.max()) #TODO: couple back to uds, currently ordering mismatch: https://github.com/Deltares/MeshKernelPy/issues/72
#cutcells
print('Cutting cells')
geometry_separator = -999
xx = np.array([147.83625 , 147.839556, 147.855833, 147.877528, 147.904139,
147.911222, 147.885861, 147.849111, 147.85625 , 147.83625 ])
yy = np.array([-40.305028, -40.301667, -40.293778, -40.30625 , -40.291222,
-40.301639, -40.326278, -40.324611, -40.309222, -40.305028])
xx = np.concatenate([xx-0.075, [geometry_separator], xx])
yy = np.concatenate([yy, [geometry_separator], yy])
delete_pol_geom = GeometryList(x_coordinates=xx, y_coordinates=yy, geometry_separator=geometry_separator) #TODO: .copy()/to_numpy() makes the array contiguous in memory, which is necessary for meshkernel.mesh2d_delete()
mk.mesh2d_delete(geometry_list=delete_pol_geom,
delete_option=DeleteMeshOption.INSIDE_NOT_INTERSECTED,
invert_deletion=False)
ortho = mk.mesh2d_get_orthogonality().values
# print('mk_object.mesh2d_get_obtuse_triangles_mass_centers()')
# print(mk_object.mesh2d_get_obtuse_triangles_mass_centers().values)
# print('mk_object.mesh2d_get_orthogonality()')
print("max ortho after cutting",ortho.max()) #TODO: couple back to uds, currently ordering mismatch: https://github.com/Deltares/MeshKernelPy/issues/72
# print('mk_object.mesh2d_get_hanging_edges()')
# print(mk_object.mesh2d_get_hanging_edges())
# mk_object.mesh2d_delete_hanging_edges()
m2d = mk.mesh2d_get()
bool_ortho_high = ortho>0.4
sel_x = m2d.edge_x[bool_ortho_high]
sel_y = m2d.edge_y[bool_ortho_high]
sel_o = ortho[bool_ortho_high]
fig, ax = plt.subplots()
m2d.plot_edges(ax=ax,linewidth=1)
ax.scatter(sel_x,sel_y,5,c=sel_o, zorder=5)
# ctx.add_basemap(ax=ax, crs=crs, attribution=False)
# dfmt.plot_coastlines(ax=ax, res='i', min_area=5, crs=crs) Results in:
|
closed with latest meshkernelpy pr |
Adding to the code in #150 (comment), one can now use the new illegalcells_geom = mk.mesh2d_get_face_polygons(num_edges=6)
illegalcells_np = np.c_[illegalcells_geom.x_coordinates, illegalcells_geom.y_coordinates]
illegalcells_np[illegalcells_np==illegalcells_geom.geometry_separator] = np.nan
fig, ax = plt.subplots()
m2d.plot_edges(ax=ax,linewidth=1)
# ax.scatter(sel_x,sel_y,5,c=sel_o, zorder=5)
ax.plot(illegalcells_np[:,0], illegalcells_np[:,1], "r-") |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Is your feature request related to a problem? Please describe.
In some rare cases, removing grid with coastline polygons deletes only one edge which result in high orthogonality. This is a known issue in workflows from the past, but this is the first time I encounter it in a grid generated with meshkernel.
The below example generates a grid with an orthogonality of max 0.0262, but after cutting it the max orthogonality is 0.418
Example (can probably be simplified to speed it up, but is also sensitive to change):
Describe the solution you'd like
Prevent high orthogonality with cutting grid.
Describe alternatives you've considered
The text was updated successfully, but these errors were encountered: