Add support for signed distance for nonuniform grids

[xylar]

This requires interpolating to "pixel" or "index" space and is only effective if the geometry is sampled at a resolution close to the finest on the grid.

To support this, some refactoring has been done so that geometry is subdivided first before computing distance and mask (then signed distance).

[xylar]

@dengwirda, could you review this at your convenience? I included an example I played around with in my testing in the documentation. The resulting distance, mask and signed-distance functions look like this:

[dengwirda]

Thanks @xylar, I checked this via an ICoM config., and all works as advertised.
In terms of computational expense though, can I check expectations? For me, forming a distance function (for a high-res case) seems to take an hour or so?

[xylar]

@dengwirda, that's disappointing. Let's chat about alternative approaches (for both ICoM and COMPASS).

Would it be hard to add some print statements to find out which part of the process is slow? If it's transforming the geometry, I could look into that. But I assume the speed issue is from finding the distance, not the sign or transforming the geometry into index space. FLANN (what we use to find nearest-neighbor distance) has been pretty performant for the tests I've done but it still gets slower the more points you have in the shapes you're using or in your reference grid.

[dengwirda]

@xylar, okay, I've tried to dig into this a little further, and am seeing the following:

• Essentially all of the time is spent determining nearest neighbour distances.
• Using FLANN this takes ~2500 sec.
• Using spatial.cKDtree (the c is the "fast" Cython version) this takes ~1200 sec.
• Using spatial.cKDtree and distance, _ = tree.query(pts, eps=0.025) (the eps=0.025 means dist is approx to within 2.5%) this takes ~200 sec.

This is for quite a big problem (|lon-lat| = 3819 x 5401) with the distance taken to a high-res USGS watershed polygon.

Since the FLANN vs spatial.cKDtree comparison is opposite to expectations, I'm worried that I've done something wrong here... Is there anything to keep in mind when installing FLANN?

Distance map (in km) from watershed polygons (white):

[dengwirda]

As an additional aside, I typically use jigsaw's marche solver to compute distance functions, but as (a PDE-based approximation to) the geodesic distance, rather than the Euclidean one.

• This takes ~40 sec for this case.

I think there's no problem using Euclidean distances, but there are some interesting differences compared to the geodesic results:

dist_3d - dist_geodesic (km)

dist_3d - dist_geodesic (clipped colourbar)

[xylar]

@dengwirda, thanks a lot of looking into this so extensively.

I don't think we are particularly attached to the methods we're using. If jigsaw can be used in place of the tools we have, it would be worth rewriting the signed-distance code to use it instead. I believe @sbrus89 used Euclidean distance, rather than geodesic distance, out of convenience, not because of any particular preference for that metric.

We could certainly add support for using cKDtree, including setting eps as an input argument if that's the easier/better way to go.

[dengwirda]

@xylar, no worries --- I'll approve this now, and put together a different PR (when this is merged) to add the geodesic option to signed_distance.py.

[xylar]

Thanks @dengwirda!