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

In cube.intersection, find split cells using a tolerant equality check #4220

Merged
merged 4 commits into from
Jul 6, 2021

Conversation

lbdreyer
Copy link
Member

@lbdreyer lbdreyer commented Jul 2, 2021

When performing a cube.intersection, Iris sometimes has to wrap the coord bounds. As an example when the coord on the original cube is between (0,360), but someone requests an intersection of (-60, 60). This can result in split cells that straddle the wrap point. For example, your original bounds may be:

array([[299.7, 299.8],
       [299.8, 299.9],
       [299.9, 300. ],
       [300. , 300.1],
       [300.1, 300.2]])

But your wrapped bounds, with wrapping point of -60, will look like:

array([[299.7, 299.8],
       [299.8, 299.9],
       [299.9, -60. ], # straddle point!
       [-60. , -59.9],
       [-59.9, -59.8]])

Iris calculates the diff between the upper and lower bounds, for bounds inside the requested intersection, resulting in the variables pre_wrap_delta and post_wrap_delta, and then it compares these to find any straddle points.
If there are any differences, it will then use these to determine where the new wrapping point (which avoids the straddle points) should be.

Notably, Iris first checks whether there are any diffs between pre_wrap_delta and post_wrap_delta with np.allclose(). If there are differences, it then determines where they are. This np.allclose check was introduced in #1641 as we had found previous issues with numerical tolerances (post_wrap_delta was almost like pre_wrap_delta).
However, if np.allclose returns False as there is a definite straddle point, in order to find where the straddle points are Iris uses np.where(pre_wrap_delta!= post_wrap_delta). But this again causes problems with numerical tolerances, and so when Iris tries to work out where the new wrapping point should be, it gets this wrong.

I've outlined an example below :

Using the bounds and wrapped bounds we have above, the pre_wrap_delta looks like:

array([[0.1],
        [0.1],
        [0.1],
        [0.1]]),

(note the intersections doesn't include the first cell (299.8,299.9) as that is outside our intersection area of (-60, 60)
and the post_wrap_delta looks like

array([[ 1.000e-01],
       [-3.599e+02],
       [ 1.000e-01],
       [ 1.000e-01]])

Those almost match but not quite:

>>> pre_wrap_delta != post_wrap_delta
array([[ True],
       [ True],
       [ True],
       [ True]])
>>> pre_wrap_delta[-1][0]
0.0999999999999659
>>> post_wrap_delta[-1][0]
0.10000000000000142

So Iris thinks they are all straddle points, and so ends up choosing a wrong new wrap point.

This PR changes the search for straddle points to instead use a tolerant equality check:

split_cell_indices, _ = np.where(
                    ~np.isclose(pre_wrap_delta, post_wrap_delta)
                )

Which allows for these floating point problems.

@lbdreyer lbdreyer requested a review from bjlittle July 2, 2021 02:04
@rcomer
Copy link
Member

rcomer commented Jul 2, 2021

So the logic introduced at #4059 assumed there would be only one split cell, and that split cell would be around the requested minimum. This PR ensures that that assumption is valid, so looks right to me.

Looking at this again makes me realise #4059 should have been slightly different though. If the check had used maximum - modulus instead of maximum % modulus, the recalculation would have happened, since -300 is not in cells. I've opened #4221 to follow up on that at some point (I don't think it's urgent).

Copy link
Member

@bjlittle bjlittle left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'd forgotten just how much Cube.intersection hurts my head 🤯

Awesome fix 🚀 💯

Just a couple of minor changes to bolster the tests and I think we're there 👍

lib/iris/tests/unit/cube/test_Cube.py Outdated Show resolved Hide resolved
@lbdreyer
Copy link
Member Author

lbdreyer commented Jul 5, 2021

Docs are broken 😩
I'll try fixing tomorrow morning!

Copy link
Member

@bjlittle bjlittle left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@lbdreyer Awesome. Just some minor changes to the whatsnew. I think my suggestions should fix the docs building.

Also, could you:

  • add a reference to the 3.0.3.rst at the top of the docs/iris/src/whatsnew/index.rst
  • bump the iris.__version__ to be 3.0.3 (might as well do this now 👍)

docs/iris/src/whatsnew/3.0.3.rst Outdated Show resolved Hide resolved
docs/iris/src/whatsnew/3.0.3.rst Outdated Show resolved Hide resolved
docs/iris/src/whatsnew/3.0.3.rst Outdated Show resolved Hide resolved
docs/iris/src/whatsnew/3.0.3.rst Outdated Show resolved Hide resolved
docs/iris/src/whatsnew/3.0.3.rst Outdated Show resolved Hide resolved
docs/iris/src/whatsnew/3.0.3.rst Outdated Show resolved Hide resolved
Copy link
Member

@bjlittle bjlittle left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@lbdreyer Fantastic, thanks 🍻

@bjlittle bjlittle merged commit 51fbe69 into SciTools:v3.0.x Jul 6, 2021
abooton pushed a commit that referenced this pull request Aug 10, 2021
* Add release highlights and pin rc version (#3898)

* Add release highlights and pin rc version

* review actions

* reorder release highlights (#3899)

Tweak release highlights

* Add whatsnew announcement (#3900)

* Fix spelling (#3903)

* Fix unit label handling (#3902)

* Add failing test of plotting

* Implement fix to pass test

* Update idiff to ignore irrelevant hyphens in path

* Update imagerepo (following docs)

* Update after review by @trexfeathers

* Add whatsnew entries

* Move whatsnew entries into correct file

* Release Docs Improvements (#3895)

* Minor phrasing change in 'Release candidate'.

* Before release deprecations.

* Whatsnew highlights section.

* Relax setup.py setup requirements (#3909)

* Updated CF saver version in User Guide and docstring (#3925)

* Updated CF saver version in User Guide and docstring

* Remove references to CF version of the loader in docstrings

* Added whatsnew

* Pin cftime<1.3.0

* Migrate to cirrus-ci (#3928)

* migrate from travis-ci to cirrus-ci

* added whatsnew entries

* ignore url for doc link check (#3929)

* whatsnew for coord default units (#3924)

* Cube._summary_coord_extra: efficiency and bugfix (#3922)

* Add Documentation Title Case Capitalization (#3940)

* Use Title Case Capitalisation for Documentation

* add whatsnew enter

* CI requirements drop pip packages (#3939)

* requirements pip to conda

* use pip install over develop

* default PY_VER to python versions

* update links (#3942)

* update links

* added s to http

* Add support for 1-d weights in collapse. (#3943)

* Remove warning for convert_units on lazy data (#3951)

* drop stickler references in docs (#3953)

* drop stickler references in docs

* remove sticker from common links

* update docs for travis-ci to cirrus-ci (#3954)

* update docs for travis-ci to cirrus-ci

* add 'travis-ci' reference locally to whatsnew

* update whatsnew comment

* docs for nox (#3955)

* docs for nox

* add titles, notices and additional detail

* review actions

* Resolve test coverage (#3947)

* test coverage for __init__ and __call__

* test coverage for metadata resolve and coverage

* partial test coverage for metadata mapping

* python 3.6 workaround for deepcopy of mock.sentinel

* test coverage for Resolve._free_mapping

* test coverage for Resolve convenience methods

* add test stub for Resolve._metadata_mapping

* fix Test__tgt_cube_position

* test coverage for shape

* test coverage for _as_compatible_cubes

* test coverage for Resolve._metadata_mapping

* test coverage for Resolve._prepare_common_dim_payload

* test coverage for Resolve._prepare_common_aux_payload

* test coverage for Resolve._prepare_points_and_bounds

* test coverage for Resolve._create_prepared_item

* test coverage for Resolve._prepare_local_payload_dim

* test coverage for Resolve._prepare_local_payload_aux

* test coverage for Resolve._prepare_local_payload_scalar + docs URL skip

* test coverage for Resolve._prepare_local_payload

* test coverage for Resolve._metadata_prepare

* added docs URL linkcheck skip

* test coverage for Resolve._prepare_factory_payload

* test coverage for Resolve._get_prepared_item

* review actions

* test coverage for Resolve.cube

* pin v3.0.0 version and whatnew date (#3956)

* update github ci checks image (#3957)

* Promote unknown units to dimensionless in aux factories (#3965)

* promote unknown to dimensionless units in aux factories

* patch aux factories to promote unknown to dimensionless units for formula terms

* add whatnew PR for entry

* Release branch prepare for v3.0.2 (#4044)

* update intersphinx mapping and matplotlib urls (#4003)

* update intersphinx mapping and matplotlib urls

* use matplotlib intersphinx where possible

* review actions

* review actions

* cirrus-ci compute credits (#4007)

* cirrus-ci conditional tasks (#4019)

* cirrus-ci conditional tasks

* use bc for bash arithmetic

* revert back to sed

* use expr

* reword

* minor documentation changes

* review actions

* prepare v3.0.2 release

* Fix test_incompatible_dimensions test (#3977)

* test_incompatible_dimensions used a ragged array for the test, which has been deprecated in numpy, and now fails if dtype is anything other than object.  This test appears to be checking that the addition of a [2x4] masked array to a [2x3] masked cube should raise a ValueError. This commit fixes the creation of `data3` object to be a [2x4] non-ragged array.

* Added entry to what's new

* Added name to core developer list :)

* Update latest.rst

Fixed space in PR macro call

* update whatsnew v3.0.2

Co-authored-by: James Penn <james@jamespenn.co.uk>

* um_stash_source attribute improved handling (#4035)

* Modified pyke rule

* Tests added

* Black and whatsnew

* Include PR number

* Remove latest.rst

* Add what's new

* Support for py38 and Cartopy 0.19 (#4130)

* mpl 3.4.1 updates (#4087)

* replace most recent hashes (#4112)

* Corrected plot_anomaly_log_colouring for new Matplotlib linscale rules. (#4115)

* Cartopy 0.19 updates (#4128)

* Use assertArrayAllClose for sqrt test (#4118)

* using AllClose for sqrt test

* Omitting the checksum from test cml

* use ArrayAllClose (rebase reset it?)

* Iris py38 (#3976)

* support for py38

* update CI and noxfile

* enforce alphabetical xml element attribute order

* full tests for py38 + fix docs-tests

* add whatsnew entry

* update doc-strings + review actions

* Alternate xml handling routine (#29)

* all xml tests pass for nox tests-3.8

* restored docstrings

* move sort_xml_attrs

* make sort_xml_attrs a classmethod

* update sort_xml_attr doc-string

Co-authored-by: Bill Little <bill.james.little@gmail.com>

* add jamesp to whatsnew + minor tweak

Co-authored-by: James Penn <james@jamespenn.co.uk>

* reinstate black and nox links

* Linkcheck update (#4104)

* Updated links

* Added login remark

* Removed extra space

* change to kick cirrus

* kick cirrus

* test verbose on cirrus

* Removed test settings.

Co-authored-by: Bill Little <bill.james.little@gmail.com>
Co-authored-by: Martin Yeo <40734014+trexfeathers@users.noreply.github.com>
Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>
Co-authored-by: James Penn <james.penn@metoffice.gov.uk>
Co-authored-by: James Penn <james@jamespenn.co.uk>
Co-authored-by: tkknight <2108488+tkknight@users.noreply.github.com>

* Intersection bounds fix (replacement PR) (#4059)

* fix intersection out of bounds point

* fix intersection out of bounds point take 2

* add referencing comment

* bootstrap ci for 3.0.x (#4154)

* add cirrus-ci docker support

* wip

* wip

* wip

* wip

* wip

* add pip

* revert docker

* fix add_weekday (#32)

* fix add_weekday

* fix black link

* fix link

Co-authored-by: Ruth Comer <10599679+rcomer@users.noreply.github.com>

* Tweak to speed up dask wrapping of netcdf variables (#4135)

* Use 'meta' in da.from_array to stop it sampling netcdf variables, which is quite slow.

* Fix PR number.

* Fix test.

* Review changes.

* Update docs/iris/src/whatsnew/3.0.2.rst

Co-authored-by: lbdreyer <lbdreyer@users.noreply.github.com>

* Fb fix cube coord arithmetic (#4159)

* fix coord with cube arithmetic

* add coord arithmetic test coverage

* add a whatsnew entry

* Pp daskfix (#4141)

* Remove workaround when dask-wrapping PP data, obsoleted by #4135.

* Remove old slice testing

* Add whats new

* move whitespace?

* missing line

Co-authored-by: lbdreyer <laura.dreyer@metoffice.gov.uk>

* add release date to v3.0.2 whatsnew (#4160)

* update readme logo img src and href (#4006) (#4216)

* In cube.intersection, find split cells using a tolerant equality check (#4220)

* In cube.intersection, find split cells using a tolerant equality check

* remove unused var

* Add bounds check to tests; add what's new

* Fix whats new; update version number

* Update 3.0.3.rst

Update release date in `3.0.3.rst`

* Wide cubestr fix 3v0vx v2 (#4233)

* Widen cube printout for long ancil or cell-measure names.

* Adjust result for fixed cube-units printout.

* Added whatsnew.

* Include newest whatsnew in index.

* Review: fix whatsnew structure.

* Fix initial sections display.

* Fix some typos in the cube maths docs (#4248)

* (More of) Wide cubestr fix 3v0vx v2 (#4238)

* Fix whatsnew for #4233.

* Move 'empty slicings' whatsnew entry to 3.0.2 section.

* unpin cftime (#4222)

* Test wrangling vs v3.0.x (#4249)

* test changes

* add awol init

* rename system test cases

* fix PartialDateTime tests

* fix typo

* add whatsnew

* Fix mergeback PR #4035

* Fix mergeback PR #4035 tests

* Update mergeback nox conda-lock files

Co-authored-by: tkknight <2108488+tkknight@users.noreply.github.com>
Co-authored-by: Zeb Nicholls <zebedee.nicholls@climate-energy-college.org>
Co-authored-by: Martin Yeo <40734014+trexfeathers@users.noreply.github.com>
Co-authored-by: Jon Seddon <17068361+jonseddon@users.noreply.github.com>
Co-authored-by: Ruth Comer <ruth.comer@metoffice.gov.uk>
Co-authored-by: Patrick Peglar <patrick.peglar@metoffice.gov.uk>
Co-authored-by: James Penn <james@jamespenn.co.uk>
Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>
Co-authored-by: James Penn <james.penn@metoffice.gov.uk>
Co-authored-by: Ruth Comer <10599679+rcomer@users.noreply.github.com>
Co-authored-by: lbdreyer <lbdreyer@users.noreply.github.com>
Co-authored-by: lbdreyer <laura.dreyer@metoffice.gov.uk>
@lbdreyer lbdreyer deleted the fix_intersection branch February 21, 2023 11:19
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants