Mesh mask error

[georgebv]

Hi Rob, thank you for writing this great tool. I am having a problem with the dhitools.mesh.Mesh() object mask() method error. When I run it for my mesh I am getting this error:

Traceback (most recent call last):
  File "E:\anaconda3\envs\py36\lib\site-packages\numpy\core\fromnumeric.py", line 51, in _wrapfunc
    return getattr(obj, method)(*args, **kwds)
AttributeError: 'list' object has no attribute 'argmax'
During handling of the above exception, another exception occurred:
Traceback (most recent call last):
  File "<input>", line 5, in <module>
  File "C:\Users\grbh\Desktop\Repositories\dhitools\dhitools\mesh.py", line 546, in mask
    poly_mask = _polygon_mask(all_poly_coords)
  File "C:\Users\grbh\Desktop\Repositories\dhitools\dhitools\mesh.py", line 902, in _polygon_mask
    max_area_idx = np.argmax(poly_areas)
  File "E:\anaconda3\envs\py36\lib\site-packages\numpy\core\fromnumeric.py", line 1037, in argmax
    return _wrapfunc(a, 'argmax', axis=axis, out=out)
  File "E:\anaconda3\envs\py36\lib\site-packages\numpy\core\fromnumeric.py", line 61, in _wrapfunc
    return _wrapit(obj, method, *args, **kwds)
  File "E:\anaconda3\envs\py36\lib\site-packages\numpy\core\fromnumeric.py", line 41, in _wrapit
    result = getattr(asarray(obj), method)(*args, **kwds)
ValueError: attempt to get argmax of an empty sequence

Funny enough, if I rerun the following 2 lines several times (5-6 times, random every time) eventually it works splendidly. I can send you my mesh if this would help with debugging.

m = dhitools.mesh.Mesh('path\to\mesh.mesh')
mask = m.mask()

Also had to change config.py since my MIKE is not on drive C. You can actually use the built-in os module to check all drives one by one until you find the MIKE SDK (or throw an exception if you do not) and set the variable dynamically during import.

[georgebv]

Got an update on the error (another side of it): when I execute these commands

m = dhitools.mesh.Mesh('path\to\mesh.mesh')
mask = m.mask()

separately (one by one) I am getting the following error trace:

Traceback (most recent call last):
  File "<input>", line 1, in <module>
  File "C:\Users\grbh\Desktop\Repositories\dhitools\dhitools\mesh.py", line 540, in mask
    all_polygons = _extract_all_polygons(be_idx)
  File "C:\Users\grbh\Desktop\Repositories\dhitools\dhitools\mesh.py", line 872, in _extract_all_polygons
    node_next = next_edge[next_edge != node_start][0]
IndexError: index 0 is out of bounds for axis 0 with size 0

I was also looking at the boolean_mask method and it looks like you can speed it up significantly by using gdal raster clipper - this operation is very fast and produces the exact same grid as your code is making at this moment. Added advantage of using gdal is that you would have easy access to gdal raster tools such as exporting to tiff.

[robjameswall]

Hi Georgii, great to hear that you're finding the tool useful. Is the mesh you're using subset of any prior mesh? I've noticed during some earlier testing that mesh subsets can cause issues with the boundary codes. That would be great if you could share your mesh file.

Thanks for the suggestion for the MIKE SDK path and the boolean_mask method. The current mask implementation is very slow so I'll try out the gdal clipper.

[georgebv]

Hi Rob, the mesh I was using is not a subset. Unfortunately I would not be able to send it to you, but I think I found a better solution - one which doesn't require using mesh at all. It goes as follows:

When you extract nodes and elements (tri or quad polygons) from the DFSU file, you can create a list of shapely.geometry.Polygon objects and then run shapely.ops.cascaded_union(polygons) on this list to get a domain polygon - it will have domain (exterior) and islands (interior) polygons in it. You can then export this polygon to shapefile using fiona and then use it to clip a raster file you generate using gdal. This operation takes seconds even for very large meshes and doesn't require the original mesh file to work. The advantage of using polygons is that you can also display the mesh using matplotlib the same way its shown in MIKE (even better since I am not a fan of MIKE's colormaps) - after all, every element is just a polygon with 3 or 4 vertices.

I had to write my own code for DFSU files interface to complete a project, but I would be interested in using your code in the future as its open source and has neat object-oriented structure.

[robjameswall]

Hi Georgii

That sounds like a great way to obtain a mesh mask, especially given the time improvement. I will try and implement this over the next few weeks. Alternatively, I would be happy to accept a pull request if you'd like to contribute.

Glad that you're finding the library useful!