Antimeridian crossing scenes have an incorrect geometry

[emiliodangelo]

stactools is creating weird geometries for scenes crossing the antimeridian.

Using the following query:

( footprint:"Intersects(POLYGON((179.85981205988597 -16.31364078726581,180 -16.31364078726581,180 -16.0913663114235,179.85981205988597 -16.0913663114235,179.85981205988597 -16.31364078726581)))" OR footprint:"Intersects(POLYGON((-180 -16.31364078726581,-179.86673370370877 -16.31364078726581,-179.86673370370877 -16.0913663114235,-180 -16.0913663114235,-180 -16.31364078726581)))" ) AND ( beginPosition:[2023-08-15T00:00:00.000Z TO 2023-08-21T23:59:59.999Z] AND endPosition:[2023-08-15T00:00:00.000Z TO 2023-08-21T23:59:59.999Z] ) AND ( (platformname:Sentinel-2 AND producttype:S2MSI2A))

In Copernicus Open Access Hub, I downloaded several scenes. Two of them, S2A_MSIL2A_20230821T221941_N0509_R029_T01KAB_20230822T021825 and S2A_MSIL2A_20230821T221941_N0509_R029_T01LAC_20230822T021825, crosses the antimeridian.

Let's take S2A_MSIL2A_20230821T221941_N0509_R029_T01KAB_20230822T021825, in its product metadata file we can find the following EXT_POS_LIST:

-16.259071119106963 180 -16.263411193996067 -179.71545 -17.2548501440442 -179.72957 -17.25048161178992 180 -17.238191708196563 179.2392098262128 -16.24776292495267 179.2586253225676 -16.259071119106963 180

Using the following code, I get what I think is the right geometry for the scene:

import antimeridian
from shapely import make_valid, Polygon
from shapely.geometry import mapping

ext_pos_list = [-16.259071119106963, 180, -16.263411193996067, -179.71545, -17.2548501440442, -179.72957, -17.25048161178992, 180, -17.238191708196563, 179.2
392098262128, -16.24776292495267, 179.2586253225676, -16.259071119106963, 180]
footprint = [p[::-1] for p in list(zip(*[iter(round(coord, 6) for coord in ext_pos_list)] * 2))]
# print(footprint) -> [(180, -16.259071), (-179.71545, -16.263411), (-179.72957, -17.25485), (180, -17.250482), (179.23921, -17.238192), (179.258625, -16.247763), (180, -16.259071)]
poly = Polygon(footprint)
# print(poly) -> POLYGON ((180 -16.259071, -179.71545 -16.263411, -179.72957 -17.25485, 180 -17.250482, 179.23921 -17.238192, 179.258625 -16.247763, 180 -16.259071))
fixed_shape = shape(antimeridian.fix_shape(poly))
# print(fixed_shape) -> MULTIPOLYGON (((180 -16.259071, 180 -16.259071, 179.258625 -16.247763, 179.23921 -17.238192, 180 -17.250482, 180 -17.250482, 180 -16.259071)), ((-180 -17.250482, -179.72957 -17.25485, -179.71545 -16.263411, -180 -16.259071, -180 -17.250482)))
valid_shape = make_valid(fixed_shape)
# print(valid_shape) -> MULTIPOLYGON (((180 -16.259071, 180 -16.259071, 179.258625 -16.247763, 179.23921 -17.238192, 180 -17.250482, 180 -17.250482, 180 -16.259071)), ((-180 -17.250482, -179.72957 -17.25485, -179.71545 -16.263411, -180 -16.259071, -180 -17.250482)))

geometry = mapping(valid_shape)
# print(geometry) = {'type': 'MultiPolygon', 'coordinates': [(((180.0, -16.259071), (180.0, -16.259071), (179.258625, -16.247763), (179.23921, -17.238192), (180.0, -17.250482), (180.0, -17.250482), (180.0, -16.259071)),), (((-180.0, -17.250482), (-179.72957, -17.25485), (-179.71545, -16.263411), (-180.0, -16.259071), (-180.0, -17.250482)),)]}
# 

Which result in this geometry:

{
  "type": "MultiPolygon",
  "coordinates": [
    [
      [
        [
          180.0,
          -16.259071
        ],
        [
          180.0,
          -16.259071
        ],
        [
          179.258625,
          -16.247763
        ],
        [
          179.23921,
          -17.238192
        ],
        [
          180.0,
          -17.250482
        ],
        [
          180.0,
          -17.250482
        ],
        [
          180.0,
          -16.259071
        ]
      ]
    ],
    [
      [
        [
          -180.0,
          -17.250482
        ],
        [
          -179.72957,
          -17.25485
        ],
        [
          -179.71545,
          -16.263411
        ],
        [
          -180.0,
          -16.259071
        ],
        [
          -180.0,
          -17.250482
        ]
      ]
    ]
  ]
}

Projecting the geometry on a map:

On the other hand, if I create a STAC item using stac sentinel2 I get the following geometry:

{
  "type": "MultiPolygon",
  "coordinates": [
    [
      [
        [
          180.0,
          -16.2612441
        ],
        [
          179.2584031582179,
          -16.259079947434127
        ],
        [
          179.23921,
          -17.238192
        ],
        [
          180.0,
          -17.250482
        ],
        [
          180.0,
          -17.250482
        ],
        [
          180.0,
          -16.2612441
        ]
      ]
    ],
    [
      [
        [
          -180.0,
          -17.250482
        ],
        [
          -179.72957,
          -17.25485
        ],
        [
          -179.71545,
          -16.263411
        ],
        [
          -180.0,
          -16.2612441
        ],
        [
          -180.0,
          -17.250482
        ]
      ]
    ],
    [
      [
        [
          179.258625,
          -16.247763
        ],
        [
          179.25840315821787,
          -16.259079947434127
        ],
        [
          180.0,
          -16.259071
        ],
        [
          179.258625,
          -16.247763
        ]
      ]
    ]
  ]
}

In a map, it looks like:

Zooming in the map, we can see the geometry is not correct and it doesn't cover the whole area:

I'm not 100% sure if the issue is in stactools-packages/sentinel2or in stac-utils/stactools.

[gadomski]

I can reproduce from current main:

stac sentinel2 create-item S2A_MSIL2A_20230821T221941_N0509_R029_T01KAB_20230822T021825.SAFE .

Produces the following geometry:

I'll look into it, thanks for the report.