Skip to content

Fix boundaries losing precision for accurate polygon calculations#677

Merged
djhoese merged 1 commit intopytroll:mainfrom
djhoese:bugfix-poly90
Aug 19, 2025
Merged

Fix boundaries losing precision for accurate polygon calculations#677
djhoese merged 1 commit intopytroll:mainfrom
djhoese:bugfix-poly90

Conversation

@djhoese
Copy link
Member

@djhoese djhoese commented Aug 18, 2025

In a real world case identified by @spruceboy, if a 32-bit coordinate is given to AreaBoundary (or any boundary class) that is exactly at 90 degrees north, the conversion to radians and then validity checks in the spherical polygon class will raise an error.

  File "/opt/cspp/polar2grid_v_3_1/libexec/python_runtime/lib/python3.11/site-packages/pyresample/spherical.py", line 571, in aedges
    yield Arc(SCoordinate(lon_start, lat_start),
              ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/opt/cspp/polar2grid_v_3_1/libexec/python_runtime/lib/python3.11/site-packages/pyresample/spherical.py", line 131, in __init__
    lon, lat = _check_lon_lat(lon, lat)
               ^^^^^^^^^^^^^^^^^^^^^^^^
  File "/opt/cspp/polar2grid_v_3_1/libexec/python_runtime/lib/python3.11/site-packages/pyresample/spherical.py", line 119, in _check_lon_lat
    _check_lat_validity(lat)
  File "/opt/cspp/polar2grid_v_3_1/libexec/python_runtime/lib/python3.11/site-packages/pyresample/spherical.py", line 111, in _check_lat_validity
    raise ValueError("Latitude values must range between [-pi/2, pi/2].")

What I discovered is that due to precision issues you end up with the input coordinate being 4.371139006309477e-08 greater than np.pi / 2 and failing the validity check. See ssec/polar2grid#766 for more info.

The alternative solution could be allowing an epsilon greater/less than 90/-90 in the validity check but then the question becomes how to handle that coordinate? Should it be clipped? Should it be left alone and hope that the math is "good enough". Maybe it is just better to be more accurate with this PR's fix.

CC @ghiggi

@spruceboy, feel free to hack your Polar2Grid installation and add the dtype=np.float64 shown in the changes here for pyresample/boundary/legacy_boundary.py.

  • Closes #xxxx
  • Tests added
  • Tests passed
  • Fully documented

@codecov
Copy link

codecov bot commented Aug 18, 2025

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 94.01%. Comparing base (e589c27) to head (827266c).
⚠️ Report is 25 commits behind head on main.

Additional details and impacted files
@@           Coverage Diff           @@
##             main     #677   +/-   ##
=======================================
  Coverage   94.01%   94.01%           
=======================================
  Files          89       89           
  Lines       13614    13621    +7     
=======================================
+ Hits        12799    12806    +7     
  Misses        815      815           
Flag Coverage Δ
unittests 94.01% <100.00%> (+<0.01%) ⬆️

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:
  • ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

@ghiggi
Copy link
Contributor

ghiggi commented Aug 18, 2025

Hey @djhoese ! It has been a while that I have worked on it.

i wonder if sooner or later we should try to replace the use of spherical class/functions with spherely. I see the development stop 6 months ago, so I don't know what is their plan for the future. But it's a wrapper around s2geometry C++ code, so migrating to spherely should also I guess speed up a lot area intersection stuffs and remapping !

@djhoese
Copy link
Member Author

djhoese commented Aug 18, 2025

Oh awesome. Yeah the people who created it are very active in this type of software.

@spruceboy
Copy link

@djhoese - Thanks for looking at this!

I am likely doing somethign wrong, but I tried those changes on my polar2grid 3.1, but I think some of the other parts of polar2grid must need updating, as it fails for me later in the process with an odd error:

DEBUG    : Reading ['/opt/cspp/polar2grid_v_3_1/libexec/python_runtime/lib/python3.11/site-packages/satpy/etc/writers/awips_tiled.yaml', '/opt/cspp/polar2grid_v_3_1/etc/polar2grid/writers/awips_tiled.yaml']
DEBUG    : Adjusting lettered grid by (377.584617, 336.026470) so it better matches data X/Y
ERROR    : Unexpected error. Enable debug messages (-vvv) or see log file for details.
DEBUG    : Unexpected error exception: 
Traceback (most recent call last):
  File "<frozen runpy>", line 198, in _run_module_as_main
  File "<frozen runpy>", line 88, in _run_code
  File "/opt/cspp/polar2grid_v_3_1/libexec/python_runtime/lib/python3.11/site-packages/polar2grid/glue.py", line 499, in <module>
    sys.exit(main())
             ^^^^^^
  File "/opt/cspp/polar2grid_v_3_1/libexec/python_runtime/lib/python3.11/site-packages/polar2grid/glue.py", line 283, in main
    ret = processor()
          ^^^^^^^^^^^
  File "/opt/cspp/polar2grid_v_3_1/libexec/python_runtime/lib/python3.11/site-packages/polar2grid/glue.py", line 337, in __call__
    return self._run_processing()
           ^^^^^^^^^^^^^^^^^^^^^^
  File "/opt/cspp/polar2grid_v_3_1/libexec/python_runtime/lib/python3.11/site-packages/polar2grid/glue.py", line 384, in _run_processing
    to_save = _save_scenes(scenes_to_save, reader_info, arg_parser._writer_args)
              ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/opt/cspp/polar2grid_v_3_1/libexec/python_runtime/lib/python3.11/site-packages/polar2grid/glue.py", line 231, in _save_scenes
    to_save = _write_scene(
              ^^^^^^^^^^^^^
  File "/opt/cspp/polar2grid_v_3_1/libexec/python_runtime/lib/python3.11/site-packages/polar2grid/glue.py", line 134, in _write_scene
    _write_scene_with_writer(scn, writer_name, data_ids, wargs, to_save)
  File "/opt/cspp/polar2grid_v_3_1/libexec/python_runtime/lib/python3.11/site-packages/polar2grid/glue.py", line 147, in _write_scene_with_writer
    res = scn.save_datasets(writer=writer_name, compute=False, datasets=data_ids, **wargs)
          ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/opt/cspp/polar2grid_v_3_1/libexec/python_runtime/lib/python3.11/site-packages/satpy/scene.py", line 1296, in save_datasets
    return writer.save_datasets(dataarrays, compute=compute, **save_kwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/opt/cspp/polar2grid_v_3_1/libexec/python_runtime/lib/python3.11/site-packages/satpy/writers/awips_tiled.py", line 1634, in save_datasets
    delayeds = self._delay_netcdf_creation(delayed_gen)
               ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/opt/cspp/polar2grid_v_3_1/libexec/python_runtime/lib/python3.11/site-packages/satpy/writers/awips_tiled.py", line 1652, in _delay_netcdf_creation
    delayed_save = dataset_to_save.to_netcdf(output_filename, mode, compute=False)
                   ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/opt/cspp/polar2grid_v_3_1/libexec/python_runtime/lib/python3.11/site-packages/xarray/core/dataset.py", line 2329, in to_netcdf
    return to_netcdf(  # type: ignore  # mypy cannot resolve the overloads:(
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/opt/cspp/polar2grid_v_3_1/libexec/python_runtime/lib/python3.11/site-packages/xarray/backends/api.py", line 1360, in to_netcdf
    dump_to_store(
  File "/opt/cspp/polar2grid_v_3_1/libexec/python_runtime/lib/python3.11/site-packages/xarray/backends/api.py", line 1407, in dump_to_store
    store.store(variables, attrs, check_encoding, writer, unlimited_dims=unlimited_dims)
  File "/opt/cspp/polar2grid_v_3_1/libexec/python_runtime/lib/python3.11/site-packages/xarray/backends/common.py", line 363, in store
    variables, attributes = self.encode(variables, attributes)
                            ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/opt/cspp/polar2grid_v_3_1/libexec/python_runtime/lib/python3.11/site-packages/xarray/backends/common.py", line 452, in encode
    variables, attributes = cf_encoder(variables, attributes)
                            ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/opt/cspp/polar2grid_v_3_1/libexec/python_runtime/lib/python3.11/site-packages/xarray/conventions.py", line 805, in cf_encoder
    new_vars = {k: encode_cf_variable(v, name=k) for k, v in variables.items()}
               ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/opt/cspp/polar2grid_v_3_1/libexec/python_runtime/lib/python3.11/site-packages/xarray/conventions.py", line 805, in <dictcomp>
    new_vars = {k: encode_cf_variable(v, name=k) for k, v in variables.items()}
                   ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/opt/cspp/polar2grid_v_3_1/libexec/python_runtime/lib/python3.11/site-packages/xarray/conventions.py", line 196, in encode_cf_variable
    var = coder.encode(var, name=name)
          ^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/opt/cspp/polar2grid_v_3_1/libexec/python_runtime/lib/python3.11/site-packages/xarray/coding/variables.py", line 472, in encode
    data -= pop_to(encoding, attrs, "add_offset", name=name)
numpy.core._exceptions._UFuncOutputCastingError: Cannot cast ufunc 'subtract' output from dtype('float64') to dtype('int8') with casting rule 'same_kind'

@djhoese
Copy link
Member Author

djhoese commented Aug 18, 2025

@spruceboy I've moved that error/issue back to the P2G issue for now. I think it is unrelated to this fix/PR.

Copy link
Member

@mraspaud mraspaud left a comment

Choose a reason for hiding this comment

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

LGTM!

@mraspaud
Copy link
Member

Indeed, spherely looks really interesting!

@djhoese djhoese merged commit 7ea8338 into pytroll:main Aug 19, 2025
28 checks passed
@djhoese djhoese deleted the bugfix-poly90 branch August 19, 2025 11:55
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

Projects

None yet

Development

Successfully merging this pull request may close these issues.

4 participants