pytroll / pytroll/pyresample

Bug where boundary coordinates wrap around 180 degrees when geostationary AreaDefinition extends beyond 180 degrees

Open
#530 7 comments 0 reactions 0 assignees View on GitHub

Nobody has claimed this yet.

Dominant language
Python
Stars
385
Forks
102
Avg merge
4d 2h
Merged PRs (30d)
9

Description

When the extent of geostationary AreaDefinitions extend beyond 180 degrees boundary coordinates wrap around 180 degrees and derived shapely.Polygons are wrong (this was discovered using the boundary method but due to bug #529 I am using get_geostationary_bounding_box_in_lonlats here instead).

from shapely.geometry import Polygon
from satpy.resample import get_area_def
from pyresample.geometry import get_geostationary_bounding_box_in_lonlats

goes_west = get_area_def("goes_west_abi_p_1km")

x, y = get_geostationary_bounding_box_in_lonlats(geos_west)
coords = list(zip(x, y))
print(coords)
Polygon(coords)

boundary_points_and_polygon
boundary_bokeh_plot

Problem description

In this case the top left corner of the boundary gets wrapped around 180 degrees and therefore the resulting polygon is wrong. It can't be used for geo selecting points and when displayed the polygon get's further simplified because I think polygon lines are not allowed to cross themselves.

The reason for this is that Proj by default wraps coordinates during transformation. This can be changed by adding +over to the proj4 CRS string. In Pyresample WKT CRS is used as should be and as far as I could find out this flag can not be directly set in WKT CRS. Instead if a CRS is initialized with +over the WKT representation gets a REMARK section with the proj4 string which seems to be respected when doing transformations with Proj.

goes_west_copy = goes_west.copy(projection=goes_west.crs.to_proj4() + " +over")

This leads to the expected output below.

Expected Output

boudary_points_and_polygon_expected
boundary_bokeh_expected_output

Versions of Python, package at hand and relevant dependencies

Python: 3.10.8
Pyresample: v1.27.1

Contributor guide

No contributing guide indexed for this repository

First steps

  1. Read the whole issue, then the project's contributing guide.
  2. Comment on the issue to say you are picking it up — it saves two people doing the same work.
  3. Fork the repository and make your change on a branch.
  4. Open a pull request that references the issue number.

Research direction

Start with get_geostationary_bounding_box_in_lonlats and the AreaDefinition/CRS transformation shown in the report, reproducing the GOES-West example. Trace how the WKT CRS is used by Proj and verify the reported +over behavior. Done means boundary coordinates extending beyond 180 degrees remain unwrapped and form the expected Shapely polygon.

Written by the indexing model from the issue text.

Assessment

Tech stack
python
Domain
data
Issue type
Bug
Difficulty
4/5
Estimated time
3-5 days
Activity status
Stale
Clarity
Mostly clear
Newbie friendliness
35/100

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.