Fix interp_points + add assertion in test

[rhugonnet]

Changed interp_points() to scipy.ndimage.map_coordinates instead of griddata. Commented the rest of the code that might allow to write a mutiprocessing loop if needed in the future.
Added assertions to test_interp().
All tests passing.
Resolves #100
Resolves #69

One current limitation: we need to know if the pixel coordinates refers to center or corner (is_area or is_point). Right now I shift x/y by 0.5 pixel assuming it's the center, which works for our test Raster.
For robustness, we'd need to add the is_area or is_point info from GDAL into the Raster class. This is necessary to further improve the inter_points method.

[rhugonnet]

EDIT: adding "Resolves" in comment doesn't seem to link to the PR.

[erikmannerfelt]

I only have some performance suggestions.

[rhugonnet]

No wonder I couldn't get my head around those tests.
I've discovered a new interesting behaviour of the rio.ds.index() function:

r.ds.index(x,y)
Out[33]: (75, 301)
r.ds.index(x,y,op=np.float32)
Out[34]: (75.0, 302.0)

Just by changing operator, you shift y by 1. Right now test_interp() works the way it is intended with np.float32. But, when comparing to value_at_coords(), it only works if the operator is math.floor (lower y index).

How should we go about this?

[erikmannerfelt]

No wonder I couldn't get my head around those tests.
I've discovered a new interesting behaviour of the rio.ds.index() function:

r.ds.index(x,y)
Out[33]: (75, 301)
r.ds.index(x,y,op=np.float32)
Out[34]: (75.0, 302.0)

Just by changing operator, you shift y by 1. Right now test_interp() works the way it is intended with np.float32. But, when comparing to value_at_coords(), it only works if the operator is math.floor (lower y index).

How should we go about this?

Hahaha what!? That should probably be added as an issue on rasterio, no? A quick search in their issues shows no result corresponding to this phenomenon.

Could this be circumvented for now, and then hopefully be fixed by rasterio in time?

[rhugonnet]

Wait for the next level:

r.ds.index(x,y,op=np.float32,precision=10)
Out[20]: (138.0, 504.0)
r.ds.index(x,y,op=np.float32)
Out[21]: (138.0, 504.0)
r.ds.index(x,y)
Out[22]: (138, 503)
r.ds.index(x,y,precision=0)
Out[23]: (137, 504)
r.ds.index(x,y,precision=1)
Out[24]: (137, 504)
r.ds.index(x,y,precision=10)
Out[25]: (138, 504)

What's problematic is that the "math.floor" argument is by default in the underlying function, and is the one showing all those changes:

def rowcol(transform, xs, ys, op=math.floor, precision=None):

which is wrapped by index

[rhugonnet]

Need to test this in more details, but it looks like the "right" solution is the np.float32 and the default setting with math.floor can be falsely rounding (with default precision=None or other) when you provide the exact coordinates:

r.bounds
Out[27]: BoundingBox(left=478000.0, bottom=3088490.0, right=502000.0, top=3108140.0)
xmin = 478000
yax=3108140
ymax=3108140
r.ds.index(xmin,ymax)
Out[31]: (0, 0)
r.ds.index(xmin,ymax,op=np.float32)
Out[32]: (0.0, 0.0)
r.ds.index(xmin,ymax,op=np.float32,precision=10)
Out[33]: (0.0, 3.637979e-12)
r.ds.index(xmin,ymax,precision=10)
Out[34]: (0, 0)
r.ds.index(xmin,ymax,precision=1)
Out[35]: (-1, 0)
r.ds.index(xmin,ymax,precision=0)
Out[36]: (-1, 0)
r.ds.index(xmin,ymax)
Out[37]: (0, 0)

[erikmannerfelt]

Wait for the next level:

r.ds.index(x,y,op=np.float32,precision=10)
Out[20]: (138.0, 504.0)
r.ds.index(x,y,op=np.float32)
Out[21]: (138.0, 504.0)
r.ds.index(x,y)
Out[22]: (138, 503)
r.ds.index(x,y,precision=0)
Out[23]: (137, 504)
r.ds.index(x,y,precision=1)
Out[24]: (137, 504)
r.ds.index(x,y,precision=10)
Out[25]: (138, 504)

What's problematic is that the "math.floor" argument is by default in the underlying function, and is the one showing all those changes:

def rowcol(transform, xs, ys, op=math.floor, precision=None):

which is wrapped by index

Call an exorcist!

This should definitely be an issue in rasterio, no?

Also, why do they need math.floor() in the first place? Maybe it's new (and rio still supports py=2.7), but round(1.1) returns 1 (type(round(1.1)) == int).

[adehecq]

I thought we had discussed this before, but I cannot find any mention to it. But I agree, the use of floor by default is a bit strange. Maybe float32 should be the default. We could add an optional argument opt to xy2ij too.

[adehecq]

@rhugonnet what is the status of this??

[adehecq]

@rhugonnet Tomorrow, your first challenge will be to merge this PR! :-P

[rhugonnet]

@erikmannerfelt @adehecq
Defaulted to np.float32 as operator. Raises warning when other operator is used.
Modified area_or_point as optional input (None, 'Area', 'Point') that defaults to GDAL metadata AREA_OR_POINT fetched through rasterio when the value is None.
Updated all tests, passing.

[erikmannerfelt]

There seems to be a merge conflict. Otherwise, it looks very good!

[erikmannerfelt]

There should be an Iterable import here:
https://github.com/GlacioHack/GeoUtils/blob/5e2b975a1a281fb38e1695e58916b4eadddeaeb1/geoutils/georaster.py#L9

[rhugonnet]

I've fixed some additional bugs after merging.
Now only need approval to merge the PR!

[erikmannerfelt]

Great job!!