Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Coreg.fit_pts(..., mask_high_curv=True) does not seem to work. #404

Closed
erikmannerfelt opened this issue Aug 15, 2023 · 4 comments · Fixed by #480
Closed

Coreg.fit_pts(..., mask_high_curv=True) does not seem to work. #404

erikmannerfelt opened this issue Aug 15, 2023 · 4 comments · Fixed by #480
Labels
bug Something isn't working

Comments

@erikmannerfelt
Copy link
Contributor

I'm getting an error when enabling the mask_high_curv=True flag (after adhering to the warning: "Warning: There is no curvature in dataframe. Set mask_high_curv=True for more robust results").

File /nix/store/kyrwb500hfcyyabb66msvq2whpvhmh0d-python3-3.10.10-env/lib/python3.10/site-packages/xdem/coreg/base.py:907, in Coreg.fit_pts(self, reference_dem, dem_to_be_aligned, inlier_mask, transform, samples, subsample, verbose, mask_high_curv, order, z_name, weights)
    903 ref_dem = reference_dem[ref_valid]
    905 if mask_high_curv:
    906     maxc = np.maximum(
--> 907         np.abs(get_terrain_attribute(tba_dem, attribute=["planform_curvature", "profile_curvature"])), axis=0
    908     )
    909     # Mask very high curvatures to avoid resolution biases
    910     mask_hc = maxc.data > 5.0

TypeError: bad operand type for abs(): 'Raster'

The code in question is here:

xdem/xdem/coreg/base.py

Lines 905 to 910 in 02902c0

if mask_high_curv:
maxc = np.maximum(
np.abs(get_terrain_attribute(tba_dem, attribute=["planform_curvature", "profile_curvature"])), axis=0
)
# Mask very high curvatures to avoid resolution biases
mask_hc = maxc.data > 5.0

I don't know if this piece of code was ever tested, or if something new has broken it! Either way, it doesn't seem to work right now, as np.abs(Raster) doesn't work it seems. @adehecq or @rhugonnet, do you know if this has ever worked? I see three potential solutions:

  • Make np.abs(Raster) work in GeoUtils.
  • Change the problematic line to: np.max(np.abs(get_terrain_attribute(tba_dem, attribute=[...]).data), axis=0) (note the .data addition to make it a masked array instead of a raster).
  • Remove this functionality here altogether and make it a Filter

Also, I don't think it's great that a maxc limit of 5.0 is hardcoded. Should that be made a keyword argument?

What do you think @rhugonnet?

@erikmannerfelt erikmannerfelt added the bug Something isn't working label Aug 15, 2023
@rhugonnet
Copy link
Member

Yes the tests are still poor for these new point functions (I intended to work a bit on that in the next PR on coregistration, but I'm not there yet). We didn't insist too much on tests in #346 as we knew we were about to rework some of the module.
I have no idea if it has ever worked, never used it.

I agree, actually I think both directions should be done:

To keep new features cleanly tested in the future (and ensure they are added to API, we update the package history, etc...), we could add an automated "PR" checklist in .github? Something like this: https://github.com/pyproj4/pyproj/blob/main/.github/PULL_REQUEST_TEMPLATE.md.
What do you think @erikmannerfelt @adehecq?

@rhugonnet
Copy link
Member

I added np.abs and np.absolute (aliases) in GlacioHack/geoutils#393, it was just a couple lines!

I remember now: the reason that line of code might have worked before is because the Raster class had the __array_interface__. Unfortunately we had to deactivate it for now because it created infinite loops when called by np.ma functions, and messed up some priorities for arithmetic functions. They are working on fixing that in NumPy by adding an interface for masked arrays (...eventually!) 😅

@erikmannerfelt
Copy link
Contributor Author

Yes the tests are still poor for these new point functions (I intended to work a bit on that in the next PR on coregistration, but I'm not there yet). We didn't insist too much on tests in #346 as we knew we were about to rework some of the module. I have no idea if it has ever worked, never used it.

I agree, actually I think both directions should be done:

* We should definitely move this to `coreg/filters` in time, with the limit a keyword argument,

* It'd be great if we added `np.abs` to GeoUtils in single-input handled functions: https://github.com/GlacioHack/geoutils/blob/main/geoutils/raster/raster.py#L68.

To keep new features cleanly tested in the future (and ensure they are added to API, we update the package history, etc...), we could add an automated "PR" checklist in .github? Something like this: https://github.com/pyproj4/pyproj/blob/main/.github/PULL_REQUEST_TEMPLATE.md. What do you think @erikmannerfelt @adehecq?

Awesome @rhugonnet, thanks for the input! Yeah a checklist would be great. Even though there's no check that the listed parts are actually implemented, I think it could work in our favour in the long term.

Also, thanks a lot for fixing the problem on the GU side!!

@erikmannerfelt
Copy link
Contributor Author

Oh, and perhaps internally we could have a routine for PRs like #346. I was the one who pushed it despite it not being absolutely finished, so I'm mostly to blame. But I suspect there may be times when this could happen again, such as when @adehecq needed to quickly implement some functionality for his workshop a year ago (?).

It would be good to have an "express train" routine to make sure that the essentials are merged, and cleanup happens soon thereafter!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
bug Something isn't working
Projects
None yet
Development

Successfully merging a pull request may close this issue.

2 participants