Should flopy.plot.PlotCrossSection include partial cells?

[javgs-bd]

The following code shows that regardless of a very small min_segment_length argument, flopy plots the cross section including only those cells that are fully contained by the defining line. It also changes the x axis limits of the xs plot to those defined by the included model cells. I see no change if I modify the min_segment_length argument (except the code crashes when it is higher than a cell length)

Xsection.ipynb

If I'm not wrong, and the intended behavior of flopy.plot.PlotCrossSection is to always include the full extent of intersected cells above a tolerace length (min_segment_length), the code below should result in at least two additional cells being included in the xs plot (the ones partially touched by the line at the extremes)

However, beyond the role of min_segment_length, I wonder if flopy.plot.PlotCrossSection should always respect the length of the provided line that defines the xsection, and only decide which cells to include or not. This is crucial when coding with flopy for numerical against conceptual xsections comparison (using HGU elevation rasters, for instance).

While it is true that there is an argument called "extent" available, this is not straightforward to use when the xsections are not aligned with the grid orthogonal directions. This occurs because flopy can change what goes into x-axis of the xsection plot (sometimes it takes x coords, sometimes y-coords, depending on the line angle).

So, I wonder if i) there is a bug with min_segment_length and ii) if the default extent of the xsection plot returned by flopy should always respect the length of the line passed to define it. If not, then this should be warned in the description of the "line" argument.

Also, it would be ideal that the docs (https://flopy.readthedocs.io/en/latest/source/flopy.plot.crosssection.html) could indicate when flopy will use ycoords in the x-axis to deploy the xsection plot. Or maybe, make it an optional arg.

[jlarsen-usgs]

@javgs-bd

The long-term behavior of the PlotCrossSection class has been to return cells with two points of intersection with a line or multiple line segments. The reason for the behavior, is intersection with a single vertex (cell corner) is a fairly common case we've seen in the line-based intersection and we need to filter that out for plotting purposes. Code improvements could be made for representing the ends of a line or collection of line segments to check if the line ends within a cell that the line intersects with add to add the vertex to it. I'll take a look at this and try to update the intersection routine to consider end points.

As far as the projection direction is concerned when working lines (x-coords or y-coords on the x-axis) this is determined by the maximum difference in the line coordinates between the x and y direction. We could add in an optional method to force the projection direction fairly easily. We can also update the documentation to clarify how the code decides on the projection direction when supplying a line.

[javgs-bd]

Thank you @jlarsen-usgs. If the user would have control on what coordinates will end up in the x-axis of the cross-section plot, I believe that the extent argument could be used to enforce the x-axis limits to correspond with the corresponding coordinates of the original line array. It would also help that the x-axis Xsection plot direction go always from the first coordinate to the second. At least for our usage case, that would allow completely rigorous comparison of conceptual vs numerical cross-sections.

From your reply I'm not sure I'm understanding the role of min_segment_length. If a line only crosses one cell side, then that cell is never included in the plot regardless of the arg value?.

[jlarsen-usgs]

In the current version of the code, the min_segment_length does not apply to cells that only have one point of intersection, because they are discarded from plotting prior to that evaluation. Once the code is updated to include cells on the ends of the line segment, the min_segment_length parameter would apply.

[jlarsen-usgs]

@ javgs-bd

Pull Request #2782 should add in the additional functionality that you've described here.

Here's a rundown of what the code updates do:

Cross-sectional lines where the end-points end in a model cell now include those endpoints in the cross-section. For example, this cross-sectional line:

will now plot as this instead of truncating any "incomplete" cells.

I also added an optional view= parameter. view="auto", which is the default maintains the original behavior of PlotCrossSection. Passing in view="x" or view="y" allows users to specify the frame of reference. For example, with the same cross section if view="x" is passed in, we get:

xc = flopy.plot.PlotCrossSection(
    modelgrid=modelgrid,
    line={"line": line},
    geographic_coords=True,
    view="x"
)