Issue running ocsmesh.Hfun function

[BahramKhazaei-NOAA]

I get an error (pasted below) running any of the code blocks that includes ocsmesh.Hfun function, for example, the following in the toturial:

hfun_obj_4 = ocsmesh.Hfun(rasters[3], hmin=1000, hmax=5000)
hfun_msh_t_4 = hfun_obj_4.msh_t()

And here is the error:

---------------------------------------------------------------------------
Exception                                 Traceback (most recent call last)
Cell In[28], line 2
      1 hfun_obj_4 = ocsmesh.Hfun(rasters[3], hmin=1000, hmax=5000)
----> 2 hfun_msh_t_4 = hfun_obj_4.msh_t()

File /scratch2/STI/coastal/save/Bahram.Khazaei/env/miniconda3/envs/myenv/lib/python3.10/site-packages/ocsmesh/hfun/raster.py:465, in HfunRaster.msh_t(self, window, marche, verbosity)
    462 if marche is True:
    463     libsaw.marche(opts, hfun)
--> 465 libsaw.jigsaw(opts, geom, window_mesh, hfun=hfun)
    467 del geom
    468 # do post processing

File /scratch2/STI/coastal/save/Bahram.Khazaei/env/miniconda3/envs/myenv/lib/python3.10/site-packages/jigsawpy/libsaw.py:746, in jigsaw(opts, geom, mesh, init, hfun)
    740 retv = JLIB.jigsaw(ct.byref(ojig),
    741                    ct.byref(gmsh),
    742                    iptr, hptr,
    743                    ct.byref(mmsh))
    745 if (retv != +0):
--> 746     raise Exception(
    747         "JIGSAW returned code: " + str(retv))
    749 #--------------------------------- copy buffers to MSH_t
    751 get_msh_t(mesh, mmsh)

Exception: JIGSAW returned code: 4

[BahramKhazaei-NOAA]

Also, I believe first line in cell 55 of the tutorial should be: geom_obj_5 = ocsmesh.Geom(rasters[1], zmax=10)

[1] is missing.

[SorooshMani-NOAA]

Thanks for testing @BahramKhazaei-NOAA. Please see my suggestions below:

ocsmesh.Hfun(rasters[3], hmin=1000, hmax=5000)

Can you please try hmin=1000.0 and hmax=5000.0 (with the added .0) and see if it works? If not it might just be a Jigsaw installation issue I guess.

Also, I believe first line in cell 55 of the tutorial should be:

No, this is correct, it is supposed to create a combination geometry from all the 4 raster files. Do you run into an issue if you don't put [1]? The warnings are OK to see. I'll need to address it on OCSMesh side at some point.

[BahramKhazaei-NOAA]

Got the same error after adding .0

Thanks for testing @BahramKhazaei-NOAA. Please see my suggestions below:

ocsmesh.Hfun(rasters[3], hmin=1000, hmax=5000)

Can you please try hmin=1000.0 and hmax=5000.0 (with the added .0) and see if it works? If not it might just be a Jigsaw installation issue I guess.

Without the [1] the code runs for hours, and after completion the next code block plots an empty figure.

Also, I believe first line in cell 55 of the tutorial should be:

No, this is correct, it is supposed to create a combination geometry from all the 4 raster files. Do you run into an issue if you don't put [1]? The warnings are OK to see. I'll need to address it on OCSMesh side at some point.

[SorooshMani-NOAA]

Got the same error after adding .0

Were you able to run any meshing commands using your environment? I think if this fails too, there must be an issue with the jigsaw library installation

Without the [1] the code runs for hours, and after completion the next code block plots an empty figure.

It is also possible sometimes that pygeos package causes an error. In your environment, if it is installed, can you try

conda remove pygeos

and also

conda install -cconda-forge "shapely<2"

Because sometimes the wrong version of these two can result in errors or delays.

[BahramKhazaei-NOAA]

I installed everything from scratch and issue with ocsmesh.Hfun is resolved. I will continue with the tutorial and if everything works we can close this issue. Thnaks!

[SorooshMani-NOAA]

OK, thanks! The whole making sure libraries work together thing is why I usually want to avoid too many different dependencies. In this case if Jigsaw was released in PyPI, we wouldn't have had this issue. Please let me know if you notice anything else

[BahramKhazaei-NOAA]

One more little thing to update in the tutorial. You provided data for NWM channel v2.1, so the path to data in the following code cell would be data/'NWM_v2.1_channel_hydrofabric_10262020/nwm_v2_1_hydrofabric.gdb:

gdf_nwm_rivers = gpd.read_file(data/'NWM_v2.0_channel_hydrofabric/nwm_v2_0_hydrofabric.gdb', layer=0) #read conus reaches
rivers_ex4 = geom_poly_5.intersection(gdf_nwm_rivers.unary_union)

[BahramKhazaei-NOAA]

In this line rivers_ex4 = geom_poly_5.intersection(gdf_nwm_rivers.unary_union) of Example 4, I got the following error:

organizePolygons() received a polygon with more than 100 parts.  The processing may be really slow.  You can skip the processing by setting METHOD=SKIP.
TopologyException: Input geom 0 is invalid: Holes are nested at -73.946092835981489 40.998176055050735
---------------------------------------------------------------------------
TopologicalError                          Traceback (most recent call last)
Cell In[48], line 2
      1 gdf_nwm_rivers = gpd.read_file(data/'NWM_v2.1_channel_hydrofabric_10262020/nwm_v2_1_hydrofabric.gdb', layer=0) #read conus reaches
----> 2 rivers_ex4 = geom_poly_5.intersection(gdf_nwm_rivers.unary_union)

File /scratch2/STI/coastal/save/Bahram.Khazaei/env/miniconda3/envs/ocsmeshenv/lib/python3.9/site-packages/shapely/geometry/base.py:695, in BaseGeometry.intersection(self, other)
    693 def intersection(self, other):
    694     """Returns the intersection of the geometries"""
--> 695     return geom_factory(self.impl['intersection'](self, other))

File /scratch2/STI/coastal/save/Bahram.Khazaei/env/miniconda3/envs/ocsmeshenv/lib/python3.9/site-packages/shapely/topology.py:73, in BinaryTopologicalOp.__call__(self, this, other, *args)
     70 if product is None:
     71     err = TopologicalError(
     72             "This operation could not be performed. Reason: unknown")
---> 73     self._check_topology(err, this, other)
     74 return product

File /scratch2/STI/coastal/save/Bahram.Khazaei/env/miniconda3/envs/ocsmeshenv/lib/python3.9/site-packages/shapely/topology.py:38, in Delegating._check_topology(self, err, *geoms)
     36 for geom in geoms:
     37     if not geom.is_valid:
---> 38         raise TopologicalError(
     39             "The operation '%s' could not be performed. "
     40             "Likely cause is invalidity of the geometry %s" % (
     41                 self.fn.__name__, repr(geom)))
     42 raise err

TopologicalError: The operation 'GEOSIntersection_r' could not be performed. Likely cause is invalidity of the geometry <shapely.geometry.multipolygon.MultiPolygon object at 0x7f56bc5e0a60>

It also took more than 2 hours for the code to crash.

[SorooshMani-NOAA]

the following code cell would be data/'NWM_v2.1_channel_hydrofabric_10262020/nwm_v2_1_hydrofabric.gdb'

Thank you. I don't know how to download version 2.0, but on my system when I ran it, I had 2.0, that's why there's this discrepancy. I can add some notes to the notebook about this. Thanks!

In this line rivers_ex4 = geom_poly_5.intersection(gdf_nwm_rivers.unary_union) of Example 4 ...

On you new environment, do you have shapely pre version 2?

[SorooshMani-NOAA]

It also took more than 2 hours for the code to crash.

This is happening on a compute node on Hera?

[BahramKhazaei-NOAA]

Here are some updates for future reference. conda remove pygeos does not work so I removed it with pip uninstall pygeos. After removing pygeos, I got the same error and long runtime for the river intersection. Also, the following code block in tutorial runs very slow (it took more than 6 hours and stopped the code):

geom_rasters = list()
for f in dem_paths:
    geom_rasters.append(ocsmesh.Raster(f))
    geom = ocsmesh.Geom(
        geom_rasters,
        base_shape=base_gs.unary_union,
        base_shape_crs=base_gs.crs,
        zmax=0)

I tried everything with a fresh installation of the package and jigsaw and no pygeos. I tried all of these on a Jupyter Notebook on Hera.

[SorooshMani-NOAA]

Thank you for testing it. I put this on my todo list to try and rerun the notebook.

[SorooshMani-NOAA]

@BahramKhazaei-NOAA I think I see what is happening with example 4. Please see https://github.com/SorooshMani-NOAA/OCSMeshTutorial/issues/4#issuecomment-1474413935

I still need to look into the second issue you mentioned above at https://github.com/SorooshMani-NOAA/OCSMeshTutorial/issues/3#issuecomment-1462201622

[SorooshMani-NOAA]

@BahramKhazaei-NOAA the other issue is related to how conda initializes PROJ environment. For some reason it sets the environment to use network, while network is not available. What you need to do to fix this is to modify the file <your_conda_installation_path>/envs/<your_env_name>/etc/conda/activate.d/proj4-activate.sh and set PROJ_NETWORK to OFF. To do that you can just comment the lines:

if [ -f "${CONDA_PREFIX}/share/proj/copyright_and_licenses.csv" ]; then
  # proj-data is installed because its license was copied over
  export PROJ_NETWORK="OFF"
else
  export PROJ_NETWORK="ON"
fi

and add

export PROJ_NETWORK="OFF"

This way every time conda is activated, it will automatically set it to OFF.

You can also manually set PROJ_NETWORK to OFF every time you activate your conda environment. Please let me know if this works for you or not.

[SorooshMani-NOAA]

@BahramKhazaei-NOAA with the above two comments, I think you should be able to go through all the OCSMesh sections of the tutorial. Please let me know if you see any other issues doing so. Thanks!

[BahramKhazaei-NOAA]

I tried to install a new venv based on the tutorial instructions, removed pygeos, installed shapely <2, and set the PROJ_NETWORK to OFF. However, still the notebook get stock at some code blocks. For example at the follwoing with a warning (or error):

geom_rasters = list()
for f in dem_paths:
    geom_rasters.append(ocsmesh.Raster(f))
    geom = ocsmesh.Geom(
        geom_rasters,
        base_shape=base_gs.unary_union,
        base_shape_crs=base_gs.crs,
        zmax=0)
print('checkpoint')
geom_poly_tidal = geom.get_multipolygon()

checkpoint CPLCreateJoinableThread() failed: Resource temporarily unavailable. CPLCreateJoinableThread() failed: Resource temporarily unavailable. CPLCreateJoinableThread() failed: Resource temporarily unavailable. CPLCreateJoinableThread() failed: Resource temporarily unavailable.

KeyboardInterrupt

[BahramKhazaei-NOAA]

@SorooshMani-NOAA do you have a working environment on Hera. I was wondering if I can use your conda env to run the tutorial?

Tagging @yunfangsun for reference here.

[SorooshMani-NOAA]

@BahramKhazaei-NOAA I have an environment that is stored at /scratch2/STI/coastal/noscrub/Soroosh.Mani/apps/miniconda3/envs/sim2 on Hera. I believe this should be what I used last time I tested the tutorial on Hera.

[BahramKhazaei-NOAA]

I was able to run the tutorial based on the sim2 virtual environment above. The following environment is the one that I installed and gives me errors:

/scratch2/STI/coastal/save/Bahram.Khazaei/env/miniconda3/envs/ocsmesh

[SorooshMani-NOAA]

Since this seems to have been an environment related issue, I'll close it. Note that now you can install ocsmesh from conda to avoid issues like this.