Bug where boundary coordinates wrap around 180 degrees when geostationary AreaDefinition extends beyond 180 degrees
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)
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
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
- Read the whole issue, then the project's contributing guide.
- Comment on the issue to say you are picking it up — it saves two people doing the same work.
- Fork the repository and make your change on a branch.
- 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