imposing boundary to small islands or structure features based on the shape "roi" file (and not DEM)

[ivicajan]

Hi,
Possibly someone already solved the issue and willing to share ideas.
I am wondering if there is an easy way to impose hard constrain to have exact shape of small islands/features (with no elements over the islands, acting as a perfect wall and not getting wet) with exact boundary based on the shape file rather than using (coarse) DEM.
I am having small islands/features that should be left as islands - empty with predefined boundary (see fig 1). I can create perfect ROI shape file (region of interest) capturing exactly what I would like to have in the mesh as hard coast, but during the meshing it is still filled with triangles and not having coastline as should.

I am using something like this (domainm is my ROI in meters):
geom = Geom(
deepcopy(rasters),
zmax=5,
base_shape=domainm.unary_union,
base_shape_crs=domainm.crs,
)
domain_poly = geom.get_multipolygon()

Tried reducing raster zmax values but doesn't help. The resulting geom (domain_poly) looks like:

And the mesh is not what should be.

Any hints are welcome ;)

Thanks,
Ivica

p.s. I can put together working notebook example.

[felicio93]

Hi @ivicajan,

If I understand correctly you don't want to use the DEM as a constrain but instead you want to use a predefined shapefile.
In that case, you should not use a zmax nor pass a list of Rasters, as that will define the geom object based on your DEMs.

Here is a comparison I made trying to replicate your issue:
If I pass a domain shapefile without passing a zmax and list of rasters:

geom = Geom(
[],
base_shape=domain.union_all(),
base_shape_crs=domain.crs,
)

my final domain_poly will look like (it will keep the original domain shapefile boundaries):

if I pass a zmax and the rasters:
geom = Geom(
deepcopy(geom_rast_list),
base_shape=domain.union_all(),
base_shape_crs=domain.crs,
zmax=15
)

my final domain_poly will look like:

Please let us know if you have any questions or encounter any issues, we will try to get back to you ASAP.

[ivicajan]

Perfect, it is exactly what I was looking for.
Tried to follow that from approach using working examples, but was stuck as I was missing the first line with the [] in the geom:
geom = Geom(
[],
base_shape=domainm.unary_union,
base_shape_crs=domainm.crs,
)

Thanks!

[ivicajan]

Closing as resolved

[felicio93]

For the sake of documentation. The definition of the Geom object alone will not ensure that the fixed constrained will be enforced.

If I understand your issue correctly, to actually do what you want would require directly using jigsawpy to create the mesh, instead of the:

driver = JigsawDriver(geom, hfun)
mesh = driver.run()

as you might be doing if you are following on the OCSMesh tech report.

instead you would be creating your mesh via:
jigsawpy.lib.jigsaw(opts,geom_m,fixed_mesh,init=init,hfun=hfun_m)

where:
opt controls the mesh gen process within jigsaw, more details on these options here:

opts = jigsawpy.jigsaw_jig_t()
opts.mesh_dims = +2
opts.hfun_hmin = np.min(hfun_m.value)/100000 / ### change accordingly, this is just a dummy value for this test
opts.hfun_hmax = np.min(hfun_m.value)/50000  ### change accordingly, this is just a dummy value for this test

fixed mesh is the empty mesh where your final mesh will be written on:

fixed_mesh = jigsawpy.jigsaw_msh_t()
fixed_mesh.mshID = 'euclidean-mesh'
fixed_mesh.ndims = +2
fixed_mesh.crs = geom_m.crs

geom_m and hfun_m are simply:

geom_m = geom.msh_t()
hfun_m = hfun.msh_t()

and init is a copy of geom_m, with 'IDtag' = -1, this will ensure that the vertices of your fixed constrains are treated as gaps:

init = deepcopy(geom_m)
init.vert2['IDtag'] = -1

My geom domain looks like:

If I zoom in to the fixed constrains area in red, the resulting mesh looks like:

Note that if my resolution is too coarse, it won't guarantee the gaps are respected, e.g.:

[ivicajan]

Hi Felicio,
thanks for the tips. I am aware of the issue if not having fine enough resolution along/inside the problematic region. To resolve that I've created an additional lines that I am using to define an increase resolution and to be able to fulfill constrain of small triangles. In that sense my weighting function looks a bit patchy:

and the grid looks fine, possibly having a bit over killed resolution.

The approach of calling jigsaw directly seems powerful but at the same time a bit complicated. If you could make notebook with the example that would be great.

Cheers,
Ivica

[felicio93]

What you have looks really good, I am glad you found a solution within OCSMesh (without having to call jigsaw directly)!

The fact your hfun is patchy should not be a problem, especially if you are creating a mesh for SCHISM.

Sure I can provide a notebook, but based on what you have it should be very straight forward for you to run it directly.
Since you already have a geom and hfun you would just have to:

import jigsawpy:

import jigsawpy
from jigsawpy import jigsaw_msh_t 

convert your hfun and geom to msh_t objects:

geom_m = geom.msh_t()
hfun_m = hfun.msh_t()

create an init:

init = deepcopy(geom_m)
init.vert2['IDtag'] = -1

an empty mesh:

fixed_mesh = jigsawpy.jigsaw_msh_t()
fixed_mesh.mshID = 'euclidean-mesh'
fixed_mesh.ndims = +2
fixed_mesh.crs = geom_m.crs

and specify your opts:

opts = jigsawpy.jigsaw_jig_t()
opts.mesh_dims = +2
opts.hfun_hmin = np.min(hfun_m.value)/100000 / ### change accordingly, this is just a dummy value for this test
opts.hfun_hmax = np.min(hfun_m.value)/50000  ### change accordingly, this is just a dummy value for this test

and finally run jigsawpy directly:

jigsawpy.lib.jigsaw(opts,geom_m,fixed_mesh,init=init,hfun=hfun_m)

to save the mesh:

ocsmesh.Mesh(fixed_mesh).write(PATH/mesh_name.2dm", format='2dm', overwrite=True)

Please let me know if that does not work.