Cannot use RiverMapper & pyDEM in the same script

[SorooshMani-NOAA]

@feiye-vims while you're working on the latest generated mesh to see if it's useful or not, I started using RiverMapper and pyDEM to learn a bit more about them. To start I tried both serial and parallel sample for each of them successfully.

My next step was to try to combine the two serial scripts so that I end up with a single script that gets DEM as input and outputs river polygons. This is where I ran into an issue.

After some digging I realized that the issue stems from this line:

RiverMeshTools/RiverMapper/RiverMapper/make_river_map.py

Line 35 in 4ee6cd1

np.seterr(all='raise')

Basically setting this error handling in RiverMapper affects the results of pyDEM. So when at the top of myscript I import:

from pyDEM.dem import dem
from RiverMapper.make_river_map import make_river_map

I run into an issue during the calculation of the DEM thalwegs:

        dem_obj.compute_river(acc_limit=catchment_area)
        gdf = dem_obj.write_shapefile(npt_smooth=thalweg_smoothing)
        gdf.to_file(out_stem)

So I was wondering if there's any reason why this error handling is set for numpy errors?

[SorooshMani-NOAA]

I'm not sure if this just means we shouldn't set the error handling or we should fix the pyDEM's failure as a result of setting this error handling. If I comment that line in make_river_map.py I get a rivermap successfully in the end.

Update
This comment has been updated to reflect the results after fixing the script issue

[feiye-vims]

Hi Soroosh, I'm glad you are testing the tools.

I set "np.seterr(all='raise') " to catch any error during numpy calculation such as division by zero or other invalid values. It is werid it can cause a different output from pyDEM without actually throwing an error in pyDEM. I can ask Linlin to check what caused the problem. @cuill

In the meantime, I'll just comment this line out, since make_river_map never throws this error in my previous tests.

Update: I've committed the change.

[SorooshMani-NOAA]

Hi, no problem! OK, I'll also let you know if I figured out why pyDEM is doing this.

In the meantime, I'll just comment this line out

Thank you, I'll let you know if during my tests I notice anything strange in RiverMapper as a results of commenting this out.

@cuill in my case the easiest way to see the issue is by adding

+ import numpy as np
+ np.seterr(all='raise') 

after importing all other packages in the pyDEM_Samples/Serial/run_serial.py example case.

[cuill]

@SorooshMani-NOAA Thanks for the information. I'll take a look.

[cuill]

It looks like the issue is caused by the difference between the python built-in function and numpy's function. For example, python's built-in max(0, 1) will return 1, while np.max(0, 1) will raise an error. The correct format is np.max([0, 1]). I am not sure how pyDEM was triggered to call numpy's functions in this case. I'll modify the code to use numpy's functions to avoid this error.

[SorooshMani-NOAA]

Thank you for investigating this. Please let me know when I can test on my side to close the ticket.

[cuill]

@SorooshMani-NOAA You can pull the latest changes and test now. The issue with NaN values was also resolved. It should only have feature_id not reach_id.

[SorooshMani-NOAA]

@cuill I'm using the latest code from the repo, but I still get the same error when I set np.seterr(all='raise') at the top of my script. The error happens during writing the dataframe to file due to invalid coordinates:

Traceback (most recent call last):
  File "/mnt/c/Users/Soroosh.Mani/WorkArea/Development/RiverMapper/own_test/dem_to_rivermap_parallel.py", line 315, in <module>
    main()
  File "/mnt/c/Users/Soroosh.Mani/WorkArea/Development/RiverMapper/own_test/dem_to_rivermap_parallel.py", line 270, in main
    shp_files = extract_thalwegs_mpi_wrap(
  File "/mnt/c/Users/Soroosh.Mani/WorkArea/Development/RiverMapper/own_test/dem_to_rivermap_parallel.py", line 75, in extract_thalwegs_mpi_wrap
    dist_gen_files = extract_thalwegs(
  File "/mnt/c/Users/Soroosh.Mani/WorkArea/Development/RiverMapper/own_test/dem_to_rivermap_parallel.py", line 130, in extract_thalwegs
    gdf.to_file(out_stem)
  File "/home/smani/miniconda3/envs/meshdev/lib/python3.10/site-packages/geopandas/geodataframe.py", line 1263, in to_file
    _to_file(self, filename, driver, schema, index, **kwargs)
  File "/home/smani/miniconda3/envs/meshdev/lib/python3.10/site-packages/geopandas/io/file.py", line 539, in _to_file
    _to_file_fiona(df, filename, driver, schema, crs, mode, **kwargs)
  File "/home/smani/miniconda3/envs/meshdev/lib/python3.10/site-packages/geopandas/io/file.py", line 568, in _to_file_fiona
    colxn.writerecords(df.iterfeatures())
  File "/home/smani/miniconda3/envs/meshdev/lib/python3.10/site-packages/fiona/collection.py", line 549, in writerecords
    self.session.writerecs(records, self)
  File "fiona/ogrext.pyx", line 1413, in fiona.ogrext.WritingSession.writerecs
RuntimeError: GDAL Error: Coordinates with non-finite values are not allowed. Failed to write record: <fiona.model.Feature object at 0x7fa3aca5c4c0>

[cuill]

@SorooshMani-NOAA Yes, I just got this error, too. In the previous test, I accidentally comment out np.seterr(all='raise') which did not raise an error. I'll look into this again when I have time.

[cuill]

OSGeo/gdal#3542