Implement ind_from_latlon using a BallTree

[leuty]

I tried using ind_from_latlon on with the new ICON-1 domain of MCH and it was kinda slow (20 Min). I implemented a new version using a BallTree.

I also added functionality to return the n closest points. That is super useful if one want to compare against instruments and avoid double penalties by slight location offsets.

[leuty]

Note, an version for spherical earth can be found here:

https://github.com/Unidata/python-workshop/blob/fall-2016/notebooks/netcdf-by-coordinates.ipynb

[victoria-cherkas]

Breaking changes (ind_from_latlon is used in scripts in icon-vis):

• verbose argument is removed
• returns a List[int] rather than int.

Also I think its nice to have an example usage section in the docstring.

@AnnikaLau what do you think about the above? Perhaps its nicer to return a List[int] and then in our icon-vis scripts we can just select ind[0].
And I personally thing we should keep the verbose option. But let me know what you think.

[AnnikaLau]

Thanks for the update/improvement! I'm still going through it but one thing I'm not sure of yet, would ind[0] be the closest of the closest points?

[leuty]

@AnnikaLau
I think the distinction is made by either passing k=n or k=[n] to cKDTree.query, see here.

Would it help if I re-implemented the verbose option for ind_from_latlon ? Then you can test.

[AnnikaLau]

@AnnikaLau I think the distinction is made by either passing k=n or k=[n] to cKDTree.query, see here.

Would it help if I re-implemented the verbose option for ind_from_latlon ? Then you can test.

Yes I think the verbose option is useful for the case k=1, which we use for the icon-vis examples. So that the users are aware of that the location of the nearest point might be quite different. So it would be great if you could re-implement it for that case.

[leuty]

@AnnikaLau
So that the users are aware of that the location of the nearest point might be quite different.

I dont understand. Does it yield a different location for your case?

[AnnikaLau]

@AnnikaLau So that the users are aware of that the location of the nearest point might be quite different.

I dont understand. Does it yield a different location for your case?

It takes the nearest point to the coordinates given. So in the case the grid cells are big, this can be quite a bit off.

[AnnikaLau]

I created a corresponding merge request for icon-vis: C2SM/icon-vis#39

[AnnikaLau]

Suggested change

	>>> lats = np.rad2deg(ds.clat.values[:])
	>>> lons = np.rad2deg(ds.clon.values[:])
	>>> lats = xr.DataArray(np.rad2deg(data.clat), coords=data.clat.coords, dims=data.clat.dims)
	>>> lons = xr.DataArray(np.rad2deg(data.clon), coords=data.clon.coords, dims=data.clon.dims)

[AnnikaLau]

The old version doesn't work anymore with the changes to the script.

[leuty]

I dont understand why exactly it stopped, but fine with me.

[AnnikaLau]

I get errors when executing the next line otherwise

[leuty]

is this still present? Cant reproduce

[leuty]

@AnnikaLau
I realized that my approach did not work with global domains. I now use another Tree Data Structure and Haversines to compute the distances. The downside is that it adds another dependency: skikit-learn.

[AnnikaLau]

Radian and degree are now mixed up for input and output. I would only give degree as input and output as it is much more intuitive.

[AnnikaLau]

Apart from one last thing, which should be fixed, I think it looks good and should be merged together with C2SM/icon-vis#39.
However, the tests are currently failing on icon-vis, so we should wait until this is fixed before we can merge it...

[AnnikaLau]

Waiting for approval of C2SM/icon-vis#39