Elements created outside of the Polygon defined as domain.

[simonweppe]

Hi there @SorooshMani-NOAA , and thanks a lot for your tool.

Im trying to mesh a domain defined as a Polygon, but I am getting some elements created outside of it ..playing with expansion rate seems to have an effect of the amount of "unwanted" elements , but I haven't found a solution to strictly enforce that mesh extent.

Any recommendations or hints ??
Thanks

domain='domain_v1.shp' # actual shoreline (with no islands)
base_gs = gpd.read_file(domain)

# load the list of rasters
dem_paths = [
            'tin_file_30m_v4_negdown.tif', # depth needs to be negative down
            ]
# load rasters file(s)
geom_rasters = [Raster(p) for p in dem_paths]

geom1 = PolygonGeom(polygon = base_gs.unary_union,crs = base_gs.crs)

hfun_rasters = [Raster(p) for p in dem_paths]
# for testing
size_mesh_min = 100
size_mesh_max = 1000

hfun = Hfun(
    deepcopy(hfun_rasters),
    base_shape=base_gs.unary_union,
    base_shape_crs=base_gs.crs,
    hmin=size_mesh_min, 
    hmax=size_mesh_max,
    #method='fast'
    )

hfun.add_constant_value(value=size_mesh_min, lower_bound=-2, upper_bound=2.0) # constant element size for size function between depth +5.0 and -2
hfun.add_contour(level=-2, expansion_rate=0.003, target_size=size_mesh_min) # element size starts at 30m at depth = -2, then expands

driver = JigsawDriver(geom=geom1, hfun=hfun)
mesh = driver.run()
mesh.write('gridv0.gr3', format='grd', overwrite=True)

[SorooshMani-NOAA]

This looks strange. Something like this shouldn't happen. The best guess I have is that during the conversion and crs-transformation of shape into PSLG for Jigsaw, something happens to the edges that results in holes in the shape.

Can you please plot the output of the geometry multipolygon?

import geopandas as gpd

shape = geom1.get_multipolygon()
gpd.GeoDataFrame(geometry=[shape]).to_file('your/path.shp')

By the way is there any specific reason you're using collector size function instead of raster based?

[simonweppe]

Here is the plot of the Polygon im using. I have attached the <.shp> file as well.
geom1.zip

By the way is there any specific reason you're using collector size function instead of raster based?
You mean using hfun = Hfun(.. rather than hfun = HfunRaster(.. ? Because I'll have probably have more rasters to use in the final one.

[SorooshMani-NOAA]

I see, it makes sense to use collector in that case. Thanks for the domain shape; I can't promise when I'll be able to get to testing it though. One suggestion I have is to simplify the shape and try again to see if the issue is resolved or not; try different tolerances and both True and False for preserve_topology when simplifying.

[simonweppe]

Thanks @SorooshMani-NOAA . Actually, something really odd was going on.

Besides that issue with the polygon I had noticed that my mesh resolution was all weird and not following the specs I was specifying at all. I had installed ocsmesh with conda (as recommended), so I created a new conda env and this time installed using pip ... Now it seems to work as expected i.e. the mesh nicely follows the contour I'm specifying. I have no idea why ...

I ran into another issue though, for a more complex case geometry case where I used

    geom = Geom(
        deepcopy(geom_rasters),
        base_shape=base_gs.unary_union,
        base_shape_crs=base_gs.crs,
        zmax=2.0) 

instead of
geom1 = PolygonGeom(polygon = base_gs.unary_union,crs = base_gs.crs)

Running the meshing I got that error, ever got something like this ? I'll try to use

xcb] Unknown sequence number while processing queue
[xcb] Most likely this is a multi-threaded client and XInitThreads has not been called
[xcb] Aborting, sorry about that.

I'll try to define a multi polygon with the islands in there, and use that as geometry, to see if that could be a workaround.

[SorooshMani-NOAA]

I'm glad that your issue is resolved. The error seems to be related to x-server. OCSMesh uses matplotlib to calculate contours, so it internally creates figures, axes, etc. and this can result in issues sometimes, see #131. I suggest that you try the workaround in that ticket, by calling matplotlib.use('agg') at the beginning of your script. If you're using Jupyter notebook this causes the plotting to fail, but at least the meshing works!!

[simonweppe]

Oh thanks ! yes indeed this solves the other issue I had ! thanks for your reactivity - greatly appreciated.