How to get ridoff elements outside the roi?

[ivicajan]

Dear all,
I'm building a mesh for a SCHISM wet-dry simulation and got almost there ;(
The problem I am facing now is that I still have elements where I don't want them (i.e. outside the roi I defined as a exclusion zone).
Possibly I am doing some stupid mistake and would appreciate if you could have a look, have some idea and save me time.

I put all inside the jupyter notebook with all other needed files to be able to run the example. Files are on my github:
https://github.com/ivicajan/OCSmesh_sandpit/tree/main

Thanks in advance,
Regards
Ivica

[SorooshMani-NOAA]

@ivicajan, the problem is that the only the exterior of base_shape input to the collector geometry type is used. During the implementation of the collector geometry, we were just looking for a way to limit the extent of the domain, not to enforce its internal shape.

I'm not sure what exactly you're trying to do, but if you want to enforce internal meshing boundary for your domain, you can first use the collector object to extract the polygon from multiple DEMs and then do boolean operations on the extracted shape and the shape you'd like to enforce and then pass the resulting (multi)polygon as a shape geometry to the jigsaw driver.

[ivicajan]

Hi,
I am trying to make a mesh with dry nodes inside the green polygon and wet within the blue. In a fact, it is doesn't matter for the meshing but I would like to keep elements within the green boundary i.e. not extending outside the red curve. Currently, it is using a raster to define boundary..
Those are the polygons I would like to limit domain:

Right now Hfun have elements ~500m (blue color) and in the region where I don't want them (i.e. outside the red polygon):

If you can see what is wrong, code is inside the notebook:
https://github.com/ivicajan/OCSmesh_sandpit/blob/main/sandpit.ipynb

Thanks!
Ivica

[SorooshMani-NOAA]

Yes, I glanced over the notebook, but I was not sure what boundary you're trying to enforce.

Right now Hfun have elements ~500m (blue color) and in the region where I don't want them (i.e. outside the red polygon):

In OCSMesh, the size function doesn't follow the boundary specification. You shouldn't worry about the hfun elements. However I noticed in your notebook you also have elements where you don't seem to want:

Playing with the notebook, this is the domain of meshing, which looks fine based on your stated intent:

But as you can see there are many small "islands" etc. in the domain which when converted to Jigsaw input might throw it off course and cause it to fail to find the right in vs out of domain regions, resulting in the mesh you get. I would suggest cleaning up the shape for such domains. Something as basic as:

gdf_polys = gpd.GeoDataFrame(geometry=[domain_poly], crs=geom.crs).explode(ignore_index=True)
gdf_main = gpd.GeoDataFrame(geometry=gdf_polys.iloc[gdf_polys.area.argmax()])

and then creating a shape geometry from it:

geom_main = Geom(gdf_main.unary_union, crs=gdf_main.crs)
...
driver = JigsawDriver(geom_main, hfun, initial_mesh=False)
mesh = driver.run()

will result in a mesh like:

Note that you'll still get the elements out of the ROI in the Hfun, but don't worry about those! I hope this resolves the issue.

[ivicajan]

Thanks a lot for your time and tips!
I am still missing something in my notebook to get to the mesh you just plot. Would it be possibly to send me your notebook to compare with what I am doing?
Thanks in advance!
Ivica

[SorooshMani-NOAA]

No problem! Please see https://github.com/SorooshMani-NOAA/OCSmesh_sandpit_support/blob/main/sandpit.ipynb

[ivicajan]

Thanks for your time!
I just cloned your OCSmesh_sandpit_support and went through the notebook but still can't get the grid as you did (without elements outside roi). I believe we have the same hfun (at least from the figure what I can see) but the mesh generation produces a different elements/mesh

Could it be because we use possibly different versions of jigsaw?
in my case JIGSAW VERSION 0.9.14

[ivicajan]

After I removed jigsawpy that I installed via conda as in your README.md and installed the latest version 1.0.0 from the github (https://github.com/dengwirda/jigsaw-python/tree/master) I got the same result at sou did.
Turns out it was jigsawpy version at the end.

Thanks a lot for your effort!

[SorooshMani-NOAA]

I see, I didn't really try out your notebook without making changes to the calculated geometry object. Good to know that the issue was actually with jigsaw version.

In the tutorial I tried to keep it simple by installing from conda. But I know there's a version lag and some missing bugfixes in the conda version compare to the GitHub version.

[ivicajan]

While you are on the topic, what would be an efficient way to get mesh coordinates back to lon/lat projection?
RIght now it is in UTM centered at the mean of the mesh?

Cheers,
Ivica

[ivicajan]

Just figured that after saving in file it is already in WGS.

[SorooshMani-NOAA]

Right now OCSMesh only writes meshes in epsg:4326 all the time, no matter what the mesh object CRS is. There's actually a ticket #57 to write in different CRS too, but for now if all you need is 4326, you should be covered when writing to disk.

To reproject a mesh, you need to get its internal representation (call mesh_obj.msh_t to get an object of type jigsaw_msh_t) and reproject it to a different crs using ocsmesh.utils.reproject(msht_obj, dest_crs). Then after reprojection you can create a new Mesh object from it.