Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

precompute finite value mask for kdtree query #2201

Merged
merged 9 commits into from
Jul 14, 2023
12 changes: 10 additions & 2 deletions lib/cartopy/img_transform.py
Original file line number Diff line number Diff line change
Expand Up @@ -268,16 +268,24 @@ def regrid(array, source_x_coords, source_y_coords, source_proj, target_proj,
target_x_points.flatten(),
target_y_points.flatten())

# Find mask of valid points before querying kdtree: scipy >= 1.11 errors
# when querying nan points, might as well use for pykdtree too.
indices = np.zeros(target_xyz.shape[0], dtype=int)
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

was I right in having the indices array default to 0? I based this on the indices[mask]=0 line after this, wasn't entirely certain.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this is fine to use 0 here as you noted because I think it is just settting it to some valid index. But, I think we still need to add these non-finite locations you found here into the mask below as well? Something like mask = ~finite_xyz | (indices >= len(xyz))

finite_xyz = np.all(np.isfinite(target_xyz), axis=-1)

if _is_pykdtree:
kdtree = pykdtree.kdtree.KDTree(xyz)
# Use sqr_dists=True because we don't care about distances,
# and it saves a sqrt.
_, indices = kdtree.query(target_xyz, k=1, sqr_dists=True)
_, indices[finite_xyz] = kdtree.query(target_xyz[finite_xyz, :],
k=1,
sqr_dists=True)
greglucas marked this conversation as resolved.
Show resolved Hide resolved
else:
# Versions of scipy >= v0.16 added the balanced_tree argument,
# which caused the KDTree to hang with this input.
kdtree = scipy.spatial.cKDTree(xyz, balanced_tree=False)
_, indices = kdtree.query(target_xyz, k=1)
_, indices[finite_xyz] = kdtree.query(target_xyz[finite_xyz, :], k=1)

mask = indices >= len(xyz)
chrishavlin marked this conversation as resolved.
Show resolved Hide resolved
indices[mask] = 0

Expand Down