Antimeridian Edge Cases #207
Description
Hi,
we have observed that some antimeridian issues are left for some Sentinel-2 tiles.
Here's a quick example: S2A_MSIL2A_20230409T002611_N0509_R102_T01WCT_20230409T040056.SAFE
Also on CDSE in the STAC Browser it is displayed incorrectly: https://browser.stac.dataspace.copernicus.eu/collections/sentinel-2-l2a/items/[S2A_MSIL2A_20230409T002611_N0509_R102_T01WCT_20230409T040056](https://browser.stac.dataspace.copernicus.eu/collections/sentinel-2-l2a/items/S2A_MSIL2A_20230409T002611_N0509_R102_T01WCT_20230409T040056?.language=de)?.language=de
We discovered it here:
scene_path = 'S2A_MSIL2A_20230409T002611_N0509_R102_T01WCT_20230409T040056.SAFE'
from stactools.sentinel2.stac import create_item
stac = create_item(scene_path)
gives
lib/python3.13/site-packages/antimeridian/_implementation.py:441: FixWindingWarning: The exterior ring of this shape is wound clockwise. Since this is a common error in real-world geometries, this package is reversing the exterior coordinates of the input shape before running its algorithm. If you know that your input shape is correct (i.e. if your data encompasses both poles), pass `fix_winding=False`.
FixWindingWarning.warn()
---------------------------------------------------------------------------
Exception Traceback (most recent call last)
Cell In[5], line 1
----> 1 stac = create_item(scene_path)
File [~/micromamba/envs/s2_dateline/lib/python3.13/site-packages/stactools/sentinel2/stac.py:145](https://portal.terrabyte.lrz.de/node/hpdar03c02s03.cos.lrz.de/36749/~/micromamba/envs/s2_dateline/lib/python3.13/site-packages/stactools/sentinel2/stac.py#line=144), in create_item(granule_href, tolerance, additional_providers, read_href_modifier, asset_href_prefix)
137 # sometimes, antimeridian and/or polar crossing scenes on some platforms end up
138 # with geometries that cover the inverse area that they should, so nearly the
139 # entire globe. This has been seen to have different behaviors on different
(...) 142 # is unreasonably large. Typical areas will no greater than 3, whereas an
143 # incorrect globe-covering geometry will have an area for 61110.
144 if (ga := geometry.area) > 100:
--> 145 raise Exception(f"Area of geometry is {ga}, which is too large to be correct.")
147 bbox = [round(v, COORD_ROUNDING) for v in antimeridian.bbox(geometry)]
149 item = pystac.Item(
150 id=metadata.scene_id,
151 geometry=shapely_mapping(geometry),
(...) 154 properties={"created": now_to_rfc3339_str()},
155 )
Exception: Area of geometry is 64800.0, which is too large to be correct.
We've also tried increasing the tolerance of geometry simplification but it didn't change: stac = create_item(scene_path, tolerance = 10)
The error comes from: make_valid_geometry()
We've checked on which tiles we get this specific error and there seems to be kind of a system:
| N | R | T | Correction Status | Count |
|---|---|---|---|---|
| N0500 | R116 | T01WCM | Corrected | 28 |
| N0509 | R116 | T01WCM | Corrected | 4 |
| N0510 | R102 | T01WCT | Corrected | 3 |
| N0509 | R102 | T01WCT | Corrected | 2 |
| N0510 | R116 | T01WCM | Corrected | 1 |
| N0511 | R036 | T37XDC | Not Corrected | 1 |
| N0511 | R061 | T48RUN | Not Corrected | 1 |
| N0511 | R093 | T35UPB | Not Corrected | 1 |
| N0511 | R102 | T01WCT | Corrected | 1 |
| N0511 | R122 | T32SQD | Not Corrected | 1 |
| N0511 | R124 | T22LFJ | Not Corrected | 1 |
| N0511 | R132 | T51SUA | Not Corrected | 1 |
It usually looks like this. One point is across the dateline. The bbox will span the globe.
I've defined a heuristic to solve this - but not sure if it is the best way or if you have more experience/better ideas:
Here's what the fix will look like:
def crosses_dateline(geometry):
"""
Check if a Polygon crosses the dateline (+/-180 degrees) by having both positive and negative longitudes,
and those longitudes are near 180 degrees. This function excludes crossings near 0 degrees.
Parameters:
geometry (shapely.geometry.Polygon): A Polygon geometry to check.
Returns:
Boolean: True if the polygon crosses the dateline (±180 longitude), otherwise False.
"""
if isinstance(geometry, Polygon):
# Get the longitudes of the exterior ring
longitudes = [coord[0] for coord in geometry.exterior.coords]
# Check if there are both positive and negative longitudes
has_positive = any(lon > 0 for lon in longitudes)
has_negative = any(lon < 0 for lon in longitudes)
# Check if the longitudes are near the dateline (close to ±180 degrees)
near_dateline = any(abs(lon) > 170 for lon in longitudes)
# The polygon crosses the dateline if it has both positive and negative longitudes,
# and at least one longitude is close to 180 degrees (±170 or more).
return has_positive and has_negative and near_dateline
else:
return False # Return False if the geometry is not a Polygon
# Example use:
# gdf_crossing = gdf[gdf['geometry'].apply(crosses_dateline)]
def correct_dateline_crossing(geometry):
"""Corrects polygons that cross the dateline. There are 3 possible cases:
Case 1: If all longitudes have the same sign, skip the geometry. Return original.
Case 2: If there's exactly one positive or one negative longitude. Replace the sign.
Case 3: If there are two positive and two negative longitudes. Calculate distance to 180. Then whichever distance is smaller replace those longitudes with 180 and flip their sign.
Parameters:
geometry (shapely.geometry.Polygon): A Polygon that needs to be corrected identified by crosses_dateline().
Returns:
geometry (shapely.geometry.Polygon): The corrected Polygon, not crossing the dateline anymore.
"""
# Ensure the geometry is a Polygon
if not isinstance(geometry, Polygon):
return geometry
# Extract the longitudes and latitudes from the geometry
coords = list(geometry.exterior.coords)
lons = [lon for lon, lat in coords]
lats = [lat for lon, lat in coords]
# Ensure the polygon's first and last coordinates are the same
if coords[0] != coords[-1]:
raise ValueError("Polygon is not closed properly; first and last coordinates must be identical.")
# Remove the last coordinate (duplicate of the first)
lons.pop()
lats.pop()
# Check the signs of all longitudes
pos_lons = [lon for lon in lons if lon > 0]
neg_lons = [lon for lon in lons if lon < 0]
# Case 1: If all longitudes have the same sign, skip the geometry
if len(pos_lons) == 0 or len(neg_lons) == 0:
return geometry
# Case 2: If there's exactly one positive or one negative longitude
elif len(pos_lons) == 1:
# Only one positive longitude, replace it with 180
new_lons = [-180 if lon > 0 else lon for lon in lons]
elif len(neg_lons) == 1:
# Only one negative longitude, replace it with -180
new_lons = [180 if lon < 0 else lon for lon in lons]
# Case 3: If there are two positive and two negative longitudes
elif len(pos_lons) == 2 and len(neg_lons) == 2:
# Calculate distances to 180 for positive and negative longitudes
pos_distances = [180 - lon for lon in pos_lons]
neg_distances = [abs(-180 + lon) for lon in neg_lons]
if sum(pos_distances) < sum(neg_distances):
# Replace positive longitudes with 180 and flip their sign
new_lons = [180 if lon > 0 else lon for lon in lons]
new_lons = [-180 if lon == 180 else lon for lon in new_lons]
else:
# Replace negative longitudes with -180 and flip their sign
new_lons = [-180 if lon < 0 else lon for lon in lons]
new_lons = [180 if lon == -180 else lon for lon in new_lons]
else:
# Other cases, return the geometry unchanged
return geometry
# Close the polygon loop by adding the first coordinate at the end
new_coords = [(lon, lat) for lon, lat in zip(new_lons, lats)]
new_coords.append(new_coords[0]) # Ensure the polygon is closed
return Polygon(new_coords)
# Example usage:
# Assuming you have a GeoDataFrame 'gdf':
# gdf['geometry'] = gdf['geometry'].apply(correct_dateline_crossing)
Could be fixed in make_valid_geometry():
sentinel2/src/stactools/sentinel2/stac.py
Line 807 in b48dbb7
if (ga := geometry.area) > 100 & (crosses_dateline(geometry)):
raise Warning(f"Geometry is crossing the dateline and area of geometry is {ga}, which is too large to be correct. It will be corrected.")
geometry = correct_dateline_crossing(geometry)
if (ga := geometry.area) > 100:
raise Execption(f"Area of geometry is {ga}, which is too large to be correct.")
if (crosses_dateline(geometry)) :
raise Exception(f"Geometry is still crossing the dateline."
# i guess antimeridian.bbox(geometry) makes a bbox and returns it as a list [xmin, ymin, xmax), ymax] (https://github.com/gadomski/antimeridian/blob/main/src/antimeridian/_implementation.py#L748)
# bbox = [round(v, COORD_ROUNDING) for v in antimeridian.bbox(geometry)]
# We also still have to do that
bbox = [round(v, COORD_ROUNDING) for v in geometry.total_bounds]