phigre
commented
6 years ago I had a brief look at the code and it probably won't be that easy. Problem is, that the indexing is not always just equal to a slicing, but for some slices (lm, lp, km, kp), first or last indices are doubled:
One option would be to enlarge the base array by one px on each side, but then we would have to update the boundaries all the time.
The sia 2d code relies a lot on indexing arrays from 0 to n-2, 1 to n, etc. I took over the code from Alex which relied on making the indexers with
np._ix:https://github.com/OGGM/oggm/blob/97e7d371ca7d83838be1df5120fc65bd295ab7b1/oggm/core/sia2d.py#L294-L303
These indexers are static and defined out of the loop, but they are inefficient because array indexers trigger advanced indexing (i.e. returning a copy). I am quite confident that they can all be replaced with slice indexers which will trigger basic indexing (i.e. returning a view of the indexed array).
I can't be 100% sure but I have the impression that we can still win some precious computing time here.
CC @phigre