Sign Inversion on Read Issue with Size Function from `grd` File

Related issues: #68 and #70

Description When running Jigsaw with a single .ll mesh, crashes with JIGSAW return code 6.

m = Mesh.open('/_app/mesh_examples/hgrid.ll',crs=4326)
domain = Geom(m).get_multipolygon()
base_gs = gpd.GeoSeries(domain, crs='EPSG:4326')
g=Geom([m],base_shape=base_gs.unary_union, base_shape_crs=base_gs.crs,zmax=0)
h=Hfun([m],base_shape=base_gs.unary_union, base_shape_crs=base_gs.crs,hmin=0,hmax=2000,method='fast')
driver = JigsawDriver(geom=g, hfun=h)
mesh = driver.run()
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/ocsmesh/ocsmesh/driver.py", line 101, in run
    libsaw.jigsaw(
  File "/_venv/lib/python3.8/site-packages/jigsawpy-1.0.0-py3.8.egg/jigsawpy/libsaw.py", line 735, in jigsaw
    raise ValueError(
ValueError: JIGSAW returned code: 6

As @SorooshMani-NOAA said in #68 , this crash can be caused by negative values on the mesh file. I checked the mesh file and there were negative values. Since this mesh file is large to replace manually the negative values, I created a simple script to read its content (it is a text file) and replace the negative values by 0.0, and create a new .ll mesh file (lets call it hgrid_fix.ll). I run again the above code to open the hgrid_fix.ll mesh file, and still crashes with the same error.

I have done a further research and what I found is, by doing Hfun(m) in my original hgrid.ll mesh, I can check hmin and hmax of the mesh:

m = Mesh.open('/_app/mesh_examples/hgrid.ll',crs=4326)
h=Hfun(m)
h.hmin
-558.9873
h.hmax
50.574429

...then I tried h.size_from_mesh():

h.size_from_mesh()
h.hmin
6.4830033545312435
h.hmax
2264.144403087892

...and the values still does not seem correct. I checked my mesh file and my hmin and hmax are the opposite (-50.574429 and 558.9873 respectively). So it seems hmin and hmax are inverted for some reason.

One solution I found to fix this is by generating a new hgrid_fix2.ll file that 'inverts' all values, i.e, all positive become negative and those who are >0 will be 0. Now, If I do Hfun(m), the values will be correct:

m = Mesh.open('/_app/mesh_examples/hgrid_fix2.ll',crs=4326)
h=Hfun(m)
h.hmin
0.0
h.hmax
558.9873

With this fix, I can now run Jigsaw with no crashes, but it's now taking too long for a mesh with 95k nodes. Not sure if this is a bug, but if we are dealing hmin and hmax parameters from Hfun as positive values, there must be a way to OCSMesh assume hmin and hmax from HFunMesh as positive values too.

@rjdmartins When reading ll or grd file, ocsmesh inverts the sign. That's mostly because of how elevations are treated internally vs elevations sign needed by the solver, in specific SCHISM:

https://github.com/noaa-ocs-modeling/OCSMesh/blob/2d4b33a5e477b0b8d72389827dc36dba20772220/ocsmesh/mesh/mesh.py#L663-L665

Can you briefly explain what is it that you'd like to achieve, then I'll fix ocsmesh. The sign inversion is not a "bug" really, it's a "feature"! But I understand that this might cause confusion and break the code.

It's noteworthy to say that if you read/write your Hfun* object to the disk that should still work, since the values are inverted both on writing and on reading back. However, if you have your values (as positive) on the disk, then there's be a problem with inversion.

If what you need is a way to read and retain sign of values from an ll or grd file, then I can fix this by adding additional arguments for read and write functions.

I'll update the ticket name to better represent the issue discussed.

@SorooshMani-NOAA What I'm trying to do is to develop a script that runs OCSMesh and generates new meshes from bathymetry user inputs such as grd/ll Meshes and tif Rasters, and user parameters like a domain, hmin, hmax and z. Also, I want to use capabilities from OCSMesh such as create contours, features and other. The example above attempts to run OCSMesh to create a new mesh from an existing mesh, and using its domain.

To avoid confusion, let's forget hgrid_fix2.ll for now. The hgrid_fix.ll (https://we.tl/t-WYrnBVFoKu) I talked before is a mesh file that rectifies the negative values by setting them zero, to have only has positive values. The problem is, running that mesh only, crashes JIGSAW because of hmin and hmax. Also, running again the example and activating verbosity on JigsawDriver says "input error":

m = Mesh.open('/_app/mesh_examples/hgrid_fix.ll',crs=4326)
domain = Geom(m).get_multipolygon()
base_gs = gpd.GeoSeries(domain, crs='EPSG:4326')
g=Geom([m],base_shape=base_gs.unary_union, base_shape_crs=base_gs.crs,zmax=0)
h=Hfun([m],base_shape=base_gs.unary_union, base_shape_crs=base_gs.crs,hmin=0,hmax=2000,method='fast')
driver = JigsawDriver(geom=g, hfun=h, verbosity=1)
mesh = driver.run()

#------------------------------------------------------------
#
#   ,o, ,o,       /                                 
#    `   `  e88~88e  d88~\   /~~~8e Y88b    e    / 
#   888 888 88   88 C888         88b Y88b  d8b  /   
#   888 888 "8b_d8"  Y88b   e88~-888  Y888/Y88b/  
#   888 888  /        888D C88   888   Y8/  Y8/     
#   88P 888 Cb      \_88P   "8b_-888    Y    Y    
# \_8"       Y8""8D                             
#
#------------------------------------------------------------
# JIGSAW: an unstructured mesh generation library.  
#------------------------------------------------------------

  JIGSAW VERSION 1.0.0

  Reading CFG. data...

**input error: HFUN-HMIN = -5.59e+02
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/ocsmesh/ocsmesh/driver.py", line 101, in run
    libsaw.jigsaw(
  File "/_venv/lib/python3.8/site-packages/jigsawpy-1.0.0-py3.8.egg/jigsawpy/libsaw.py", line 735, in jigsaw
    raise ValueError(
ValueError: JIGSAW returned code: 6

The sign inversion is not a "bug" really, it's a "feature"! But I understand that this might cause confusion and break the code.

My confusion is the Hfun examples from the documentation assume that hmin and hmax parameters should be positive values. But the bathymetry mesh I'm loading has positive values, but OCSMesh assumes hmin negative.

Interestingly, bathymetry Rasters do not have any problems being processed, and they are negative values. Although, when I try to load them in Hfun(r), they do not show hmin and hmax values (no error, just nothing), which in this case I can't tell you what are their values.

r = Raster('/_app/raster_examples/gebco/gebco_2022_sub_ice_n90.0_s0.0_w-90.0_e0.0_compressed.tif')
hr = Hfun(r)
hr.hmin
hr.hmax

Hi @rjdmartins thank you for your feedback. Some of the internals of OCSMesh is not documented in full details in the manual. So please let me explain them here, and hopefully we'll find a solution for the OCSMesh issues you're dealing with.

First of all let me stress the following "features" of ocsmesh:

When writing or reading grd or ll files, the sign of the values, e.g. elevation values is inverted. This is true for both HfunMesh or Mesh objects. This means "values" you see in the mesh object is negative of what you see on the file on disk.
For 2dm read/write this is not the case and sign is retained.
All size function objects pass the size information to jigsawpy in the form of an unstructured grid (through the return value of msh_t() method), whether the input to the size function object is raster (like HfunRaster) or mesh (like HfunMesh).

These are crucial to navigating some of the issues when reading/writing hfunmesh or mesh from disk. In fact when writing hfunmesh_obj to disk you need to call hfunmesh_obj.mesh.write(...).

To see the sign inversion for yourself, please look at the mesh values after you open it in the above code:

m = Mesh.open('/_app/mesh_examples/hgrid_fix.ll',crs=4326)
m.msh_t.value  # <-- Look at these values

and you see all of them are inverted compared to what you had in the file.

As for the following:

Interestingly, bathymetry Rasters do not have any problems being processed, and they are negative values. ...

Note HfunMesh and HfunRaster process their inputs in fundamentally different ways. HfunCollector follows the same different logics for processing rasters vs mesh inputs in the list you provide to it. So I'll just explain what happens for HfunMesh vs HfunRaster, and the same logic goes for mesh vs raster inputs for HfunCollector.

When an object of type HfunRaster is created, the raster values (e.g. elevations) are used just as guide to find contours or apply bathymetry based sizes. The actual values that are used for element size specification are not well defined until you apply a control, such as add_contour or add_patch, etc. Also note that for HfunRaster, the minimum and maximum sizes (hmin and hmax) are not "found", they are inputs that the user has to provide. These input values will constraint the hfun object when generating the final element-size unstructured-grid hfun_raster.msh_t(). So if they are not provide (no constraint) they will simply be None, which is what you see.

When an object from HfunMesh is created, the Mesh object that was passed as input is the actual object that is used for storing the data. So all of the values of Mesh object is used directly as sizes of HfunMesh element size. Any modification of these values through the HfunMesh object is applied to the underlying mesh (e.g. if you call .size_from_mesh()). Also note that in this case the hmin and hmax are not inputs, but are taken from the underlying mesh object and nothing is "assumed". So if the mesh values are negative, hmin and hmax are going to be negative.

I hope this helps clearing some of the confusion around why ocsmesh behaves a certain way.

As for your script:

m = Mesh.open('/_app/mesh_examples/hgrid_fix.ll',crs=4326)
domain = Geom(m).get_multipolygon()
base_gs = gpd.GeoSeries(domain, crs='EPSG:4326')
g=Geom([m],base_shape=base_gs.unary_union, base_shape_crs=base_gs.crs,zmax=0)
h=Hfun([m],base_shape=base_gs.unary_union, base_shape_crs=base_gs.crs,hmin=0,hmax=2000,method='fast')
driver = JigsawDriver(geom=g, hfun=h, verbosity=1)
mesh = driver.run()

Depending on what is it exactly that you want to achieve I'd do something like:

m = Mesh.open('/_app/mesh_examples/hgrid_fix.ll',crs=4326)
+ m.size_from_mesh()
domain = Geom(m).get_multipolygon()
base_gs = gpd.GeoSeries(domain, crs='EPSG:4326')
g=Geom([m],base_shape=base_gs.unary_union, base_shape_crs=base_gs.crs,zmax=0)
h=Hfun([m],base_shape=base_gs.unary_union, base_shape_crs=base_gs.crs,hmin=0,hmax=2000,method='fast')
driver = JigsawDriver(geom=g, hfun=h, verbosity=1)
mesh = driver.run()

or use a bunch of controls to define mesh sizes:

m = Mesh.open('/_app/mesh_examples/hgrid_fix.ll',crs=4326)
domain = Geom(m).get_multipolygon()
+ m.add_feature(feature=my_feature_line, expansion_rate=my_exapnsion_value, target_size=fine_element_size)
+ ...
base_gs = gpd.GeoSeries(domain, crs='EPSG:4326')
g=Geom([m],base_shape=base_gs.unary_union, base_shape_crs=base_gs.crs,zmax=0)
h=Hfun([m],base_shape=base_gs.unary_union, base_shape_crs=base_gs.crs,hmin=0,hmax=2000,method='fast')
driver = JigsawDriver(geom=g, hfun=h, verbosity=1)
mesh = driver.run()

Also if you only want the input mesh to be automatically treated as boundary, you can pass it as base_mesh argument.

By the way, if you need to use the mesh object stored in m for both mesh size and geometry (two different purposes) I suggest you actually make a copy:

+ from copy import deepcopy
m = Mesh.open('/_app/mesh_examples/hgrid_fix.ll',crs=4326)
+ m_for_geom = deepcopy(m)
domain = Geom(m).get_multipolygon()
+ ...
base_gs = gpd.GeoSeries(domain, crs='EPSG:4326')
- g=Geom([m],base_shape=base_gs.unary_union, base_shape_crs=base_gs.crs,zmax=0)
+ g=Geom([m_for_geom],base_shape=base_gs.unary_union, base_shape_crs=base_gs.crs,zmax=0)
h=Hfun([m],base_shape=base_gs.unary_union, base_shape_crs=base_gs.crs,hmin=0,hmax=2000,method='fast')
driver = JigsawDriver(geom=g, hfun=h, verbosity=1)
mesh = driver.run()

Please let me know if you still see any issues.

Hi @SorooshMani-NOAA I think you need to rectify your solution, Mesh does not have the ability to call size_from_mesh and add_feature, these are from Hfun class. In any case, the only way is by doing deepcopy of my mesh (to avoid corruption) and calling it as Hfun, then h.size_from_mesh(). Literally more or less this:

m = Mesh.open('/_app/mesh_examples/hgrid.ll',crs=4326)
h1= Hfun(deepcopy(m))
h1.size_from_mesh()
domain = Geom(m).get_multipolygon()
base_gs = gpd.GeoSeries(domain, crs='EPSG:4326')
g=Geom([m],base_shape=base_gs.unary_union, base_shape_crs=base_gs.crs,zmax=0)
h=Hfun([h1],base_shape=base_gs.unary_union, base_shape_crs=base_gs.crs,hmin=0,hmax=2000,method='fast')
(...do some add_feature stuff here...)
driver = JigsawDriver(geom=g, hfun=h, verbosity=1)
mesh = driver.run()

@rjdmartins you're right, I mixed up mesh and hfun-mesh when I was writing the comment. Excuse me for the confusion. Also for the corruption issue you've got the right solution; as I mentioned the hfun mesh modifies the values within the mesh object (no copy) so if you want to keep the values, you'd need to deepcopy it as you said.

noaa-ocs-modeling / OCSMesh

Sign Inversion on Read Issue with Size Function from `grd` File #71