Selecting with slice returns nodes out of the bounding rectangle

[janfer95]

Hi Xugrid-Team,

Thanks for this awesome package!

I'm using it quite frequently to extract subgrids from a larger study area by making use of uda.ugrid.sel(x=slice(x1, x2), y=slice(y1, y2)). Usually this works well, but in certain cases this results in nodes (and thus triangles) that are just out of the bounding rectangle.

Looking into the source code I saw that this is due to the fact that the locate_bounding_box method of the Ugrid2d class uses the face centroids, not the face coordinates. So it happens that a face centroid might still be in the bounding rectangle, while one of the nodes might not.

Is there a way to choose the face coordinates instead? Is there a reason to prefer face centroids over face coordinates?

Thank you!

[Huite]

Hi @janfer95,

Happy to hear xugrid is useful to you!

Good question regarding the selection. I think there's two reasons that I originally chose this approach:

• It kinda matches what xarray generally does: it treats coordinates as midpoints. So if you have raster data, .sel(x=slice(x1, x2) will return a result that technically includes vertices that are outside of the selection window.
• It's cheaper too: we only need to check two coordinates (of the centroid) rather than all the node coordinates.

However, at the least, there should be an option to easily use the node coordinates instead. We'll have to do some thinking about it.

I'll get back to you and post a snippet how you can (relatively easily) select based on node coordinates instead.

[janfer95]

Thanks for the quick answer, @Huite!

Matching xarray functionalities and faster computations are indeed some good reasons to use centroids as a default.

As for an option to easily use the node coordinates instead, this should probably go into a separate method, right? It would be possible (and straightforward) to use an extra argument in the .sel() method, but I'm not sure if this does not deviate too much from "typical xarray" style.

The drawback of a separate method is that it would copy a lot of code and would be more difficult to maintain, since the only thing that really changes is the face_index that is passed to the topology_subset method.

For now I can always do the following to create a strict subgrid (and if needed extend this to get the actual data on the grid):

from xugrid import Ugrid2d

def sel_nodes_in_box(
    grid: Ugrid2d,
    xmin: float,
    xmax: float,
    ymin: float,
    ymax: float
) -> Ugrid2d:
    face_node_coordinates = grid.face_node_coordinates
    face_node_x = face_node_coordinates[:, :, 0]
    face_node_y = face_node_coordinates[:, :, 1]

    face_index = np.nonzero(
        (
            (face_node_x >= xmin)
            & (face_node_x < xmax)
            & (face_node_y >= ymin)
            & (face_node_y < ymax)
        ).all(axis=1)
    )[0]

    return grid.topology_subset(face_index)

In the codebase the only thing that really would be added / modified would be locate_bounding_box, that could possibly take an additional argument to compute either centroid or node face_indexes.

Let me know what you think is the best approach to implement this. If it's easier for you I can also add this myself then and create a pull request.

[janfer95]

Hi @Huite ! Are there any updates on this issue? If you tell me which approach to implement you prefer / is best I can implement it and open a pull request myself.