Intersecting lines with GridIntersect returns vertices as separate numpy x, y arrays if grid is offset

[mkennard-aquaveo]

• flopy version 3.3.3

In t065_test_gridintersect.py, test_rect_grid_linestring_in_and_out_of_cell can be modified to show the resulting vertices, which all works fine:

def test_rect_grid_linestring_in_and_out_of_cell():
    # avoid test fail when shapely not available
    try:
        import shapely
    except:
        return
    gr = get_rect_grid()
    ix = GridIntersect(gr, method="structured")
    result = ix.intersect(
        LineString([(5., 9), (15., 5.), (5., 1.)]))
    assert len(result) == 2
    assert result.cellids[0] == (1, 0)
    assert result.cellids[1] == (1, 1)
    assert result.vertices[0] == [(5.0, 9.0), (10.0, 7.0), (10.0, 3.0), (5.0, 1.0)]
    assert result.vertices[1] == [(10.0, 7.0), (15.0, 5.0), (10.0, 3.0)]
    assert np.allclose(result.lengths.sum(), 21.540659228538015)
    return result

But if I create a new test with the only difference being that the grid has an xyoffset, the resulting vertices are numpy arrays with the Xs all grouped together and Ys all grouped together:

def test_rect_grid_linestring_in_and_out_of_cell_xyoffset():
    # avoid test fail when shapely not available
    try:
        import shapely
    except:
        return
    gr = get_rect_grid(xyoffset=1.0)  # <- The only thing different is this xy offset
    ix = GridIntersect(gr, method="structured")
    result = ix.intersect(
        LineString([(5., 9), (15., 5.), (5., 1.)]))
    assert len(result) == 2
    assert result.cellids[0] == (1, 0)
    assert result.cellids[1] == (1, 1)
    assert result.vertices[0] == [(5.0, 9.0), (11.0, 6.6), (11.0, 3.4), (5.0, 1.0)]  # Fails. [array([5., 11., 11., 5.]), array([9., 6.6, 3.4, 1.])]
    assert result.vertices[1] == [(11.0, 6.6), (15.0, 5.0), (11.0, 3.4)]  # Fails. [array([11., 15., 11.]), array([6.6, 5., 3.4])]
    assert np.allclose(result.lengths.sum(), 21.540659228538015)
    return result

I think the problem is in this code which shows up in a couple of places:

            v_realworld = []
            for pt in vertices:
                pt = np.array(pt)
                rx, ry = transform(
                    pt[:, 0],
                    pt[:, 1],
                    self.mfgrid.xoffset,
                    self.mfgrid.yoffset,
                    self.mfgrid.angrot_radians,
                    inverse=False,
                )
                v_realworld.append([rx, ry])
            vertices = v_realworld

Maybe "v_realworld.append([rx, ry])" should be "v_realworld.append([zip(rx, ry)])"?

[dbrakenhoff]

Thanks for submitting this issue. I think you're right, in this case it should be:

v_realworld.append(list(zip(rx, ry)))

It also works (the vertices are returned correctly) when keepzerolengths=True, and the linestring ends on the boundary between two grid cells. So I'll submit the change.

It might be a good idea to take a look at that bit of code that handles tranforming and rotating the results. It is repeated several places and moving it into a function is probably smart? I don't have time for that this week, but I can take a look soon.