Fix hashtable.query by checking for both roots in index_search

[erikvansebille]

Description

While working on the Delft3D with NaNs for lon//lat, I realised that the issue wasn't the creation of the hashtable itself (that works fine with NaNs for lon/lat points), but that the point-in-cell test always returned False, even if a point was clearly in a cell.

After some further digging, I discovered that the issue was in index_search._bilinear_inverse(), which did not return the correct barycentric coordinates

With some help from an LLM (see prompt below), I then found the issue and a solution:

The function produces two quadratic roots for eta, but the current implementation always selects the
(−bb+det)/(2a) solution. It should of course also check the (−bb-det)/(2a) solution.

This PR implements that. And now the hashtable search works on the Delft3D grid!

Checklist

• Closes None
• Tests added
• This PR targets the correct branch (main for normal development, v3-support for v3 support)

AI Disclosure

• This PR contains AI-generated content.
  • I have tested any AI-generated content in my PR.
  • I take responsibility for any AI-generated content in my PR.
  • Describe how you used it (e.g., by pasting your prompt):

"I use the _bilinear_inverse function here to find out the baricentric coordinates of a point_in_cell check. However, it does not seem to be working because if I input px: [[39878.2421875 39849.01953125 39819.69921875 39790.2734375 ]
[40011.484375 39982.2109375 39952.8359375 39923.3515625 ]
[39982.2109375 39952.8359375 39923.3515625 39893.75 ]
[39849.01953125 39819.69921875 39790.2734375 39760.71875 ]], py: [[450116.90625 450081.46875 450045.90625 450010.3125 ]
[450014.3125 449978.875 449943.3125 449907.6875 ]
[449978.875 449943.3125 449907.6875 449871.90625]
[450081.46875 450045.90625 450010.3125 449974.5625 ]], xq: [39900. 39900. 39900. 39900.], yq: [450000. 450000. 450000. 450000.]
then I get xsi: [-1481.18252686 -1380.11892455 -762.69672299 -868.69661282], eta: [ 4291.41378809 3983.77629991 6657.93389165 10363.58362922], but none of these are between 0 and 1, as I would expect because the point lays within one of the cells"

[erikvansebille]

UPDATE: Actually, this fix turns out not to be needed. I had accidentally swapped the definition of the X and Y axes (called M and N in Delft3D). If I use

ds["N"].attrs["axis"] = "Y"
ds["M"].attrs["axis"] = "X"

then the index_search works also without the fix in this PR

So the question is then: should we keep this fix? Or better to not merge it into main? Note that the code does become slower. Running with this PR, a simulation took 65 minutes, while it took 41 minutes on main with the axes-swapping