Schism plotting fails on open boundary greater than 1

[benjaminleighton]

The below code generates a synthetic mesh that cause a failure on plotting the grid.plot. We are seeing this problem on real meshes. A fix is proposed in the branch schism_plotting_fix but I'm not sure if it removes existing functionality and is as yet untested.

def create_sample_gr3(filename='sample.gr3'):
    # Simpler mesh for testing
    nx, ny = 5, 4
    x = np.linspace(0, 100, nx)
    y = np.linspace(0, 80, ny)
    X, Y = np.meshgrid(x, y)
    
    nodes_x = X.flatten()
    nodes_y = Y.flatten()
    depths = np.ones_like(nodes_x) * 10
    
    # Create elements
    elements = []
    for j in range(ny-1):
        for i in range(nx-1):
            node1 = j*nx + i
            node2 = j*nx + i + 1
            node3 = (j+1)*nx + i
            node4 = (j+1)*nx + i + 1
            
            elements.append([node1+1, node2+1, node3+1])
            elements.append([node2+1, node4+1, node3+1])
    
    logger.debug(f"Creating mesh with {len(nodes_x)} nodes and {len(elements)} elements")
    
    # Store the content first to debug
    content = []
    
    # Header
    content.append('simple_mesh')
    
    # Number of elements and nodes
    content.append(f'{len(elements)} {len(nodes_x)}')
    
    # Nodes
    for i in range(len(nodes_x)):
        content.append(f'{i+1} {nodes_x[i]:.1f} {nodes_y[i]:.1f} {depths[i]:.1f}')
    
    # Elements
    for i, elem in enumerate(elements):
        content.append(f'{i+1} 3 {elem[0]} {elem[1]} {elem[2]}')
    
   # Store the content first to debug
    content = []
    
    # Header
    content.append('simple_mesh')
    
    # Number of elements and nodes
    content.append(f'{len(elements)} {len(nodes_x)}')
    
    # Nodes
    for i in range(len(nodes_x)):
        content.append(f'{i+1} {nodes_x[i]:.1f} {nodes_y[i]:.1f} {depths[i]:.1f}')
    
    # Elements
    for i, elem in enumerate(elements):
        content.append(f'{i+1} 3 {elem[0]} {elem[1]} {elem[2]}')
   # Open boundaries
    content.append('2')  # Number of open boundaries
    content.append('')   # Empty line that will be skipped
    
    # Left boundary
    left_nodes = list(range(1, ny+1))
    content.append(f'{len(left_nodes)} 2')  # number of nodes and boundary type
    for node in left_nodes:
        content.append(f'{node} 2')
        
    # Right boundary
    right_nodes = list(range(nx*(ny-1)+1, nx*ny+1))
    content.append(f'{len(right_nodes)} 2')
    for node in right_nodes:
        content.append(f'{node} 2')
        
    # Land boundaries
    content.append('0')  # Number of land boundaries
    content.append('')   # Empty line for consistency
    
    # Write content to file - ensure no empty lines at the end
    with open(filename, 'w') as f:
        f.write('\n'.join(content))  # Join with newlines
        f.write('\n')  # Single final newline
            
    return pathlib.Path(filename).absolute()

# After creating the file, let's also examine it directly
def examine_file_detailed(filename):
    """Examine file contents in detail"""
    print("\nDetailed file examination:")
    with open(filename, 'rb') as f:
        content = f.read()
        print(f"Total file size: {len(content)} bytes")
        lines = content.split(b'\n')
        print(f"Number of lines: {len(lines)}")
        for i, line in enumerate(lines):
            print(f"Line {i+1:3d}: Length={len(line):3d} Hex={line.hex()} Text='{line.decode()}'")

# Create and examine the mesh
mesh_file = create_sample_gr3()
examine_file_detailed(str(mesh_file))

# Try to read with pyschism
print("\nAttempting to read with pyschism...")
try:
    grid = Hgrid.open(mesh_file)
    print("Successfully read grid!")
except Exception as e:
    print(f"\nError reading grid: {str(e)}")
    # Find where it failed
    with open(mesh_file, 'r') as f:
        lines = f.readlines()
        header = lines[0]
        ne, np = map(int, lines[1].split())
        current_line = 2 + np + ne  # Skip header, size line, nodes, and elements
        print(f"\nFailed while trying to parse line {current_line + 1}:")
        print(f"Line content (hex): {lines[current_line].encode().hex()}")
        print(f"Line content (text): '{lines[current_line].strip()}'")

[tomdurrant]

Sorry to miss this @benjaminleighton, looks good to me. Can you merge the latest from main, and then sent me a PR and I'll double check?