search

mdbartos / pysheds

:earth_americas: Simple and fast watershed delineation in python.
GNU General Public License v3.0
717 stars 196 forks source link

nearest_cell should find index of the cell that the pour point falls in #110

Closed debboutr closed 4 years ago

debboutr commented 4 years ago

Hey, loving this package, thank you. I have a use-case where the cell that gets picked on the flow direction raster is very specific and am realizing that the nearest_cell method will find the nearest top/left cell point and not produce the cell index of the one that the point falls within. Here's an image that will help demonstrate. pour_point_edit

np.around(~affine * (x, y))
Out[210]: array([ 757., 1823.])

np.floor(~affine * (x, y))
Out[211]: array([ 756., 1822.])

I would think that this is what is needed in this method, but I am only just becoming familiar with this package. A possible fix to this could be to just switch np.around to np.floor

# line 580 of grid.py
col, row = np.floor(~affine * (x, y)).astype(int)

maybe it could be an arg to the method to determine which numpy method to use? I would think that using floor would always be the preferred method, but there are always so many use cases with this that I'd imagine that I'm missing. Let me know if this makes sense. Thanks!

mdbartos commented 4 years ago

Glad to hear this package was useful for you, and thanks for taking such a detailed look at this problem.

I think it would make sense to use some sort of keyword argument to specify behavior. Perhaps instead of changing the function, the argument could control whether to use the cell corner or centroid. So maybe something like:

np.around(~affine * (x + c, y - c))

If you're up for it, this could be a good pull request.

Thanks! MDB

debboutr commented 4 years ago

I'll give it a shot. I don't think what you've described above will work out too well though. I like the concept, but we don't know where x and y lie in the projection so negative vs. positive values would have different results. Mind if I push a branch up here? or should I fork and then make a PR that way?

mdbartos commented 4 years ago

Greetings,

This tutorial gives a good overview of the affine module: https://www.perrygeo.com/python-affine-transforms.html

I think the problem is that you want the cell reference to be the center of the pixel, but affine is using the upper-left corner. You can get around this by adding an offset as shown in the example above. Would need to test though.

To contribute, you can fork and then submit a PR to the development branch (I just pulled the most recent changes into development).