/bounds of objects that cross the Antimeridian
jonaraphael opened this issue · 7 comments
Typically objects are compact single polygons that return reasonable values [minx, miny, maxx, maxy] when you get their bounding box.
Case of interest:
When an object crosses the antimeridian (i.e. ±180º longitude), the typical behavior is splitting it into a multipolygon with disconnected elements on opposite sides of the map.
Current behavior:
The current implementation of the /bounds endpoint does not treat this correctly... It returns a degenerate bbox that looks like [-180, miny, 180, maxy], spanning the whole globe
Preferred behavior:
Instead, according to the GeoJSON spec: https://datatracker.ietf.org/doc/html/rfc7946#section-5.2 , the bbox should have the "wrong" order for the longitude coordinates (i.e. minx > maxx).
Note that it is not simply enough to return[180, miny, -180, maxy]. The endpoint should actually find the leftmost bound of the rightmost polygon as minx, and the rightmost point of the leftmost polygon as maxx.
A frontend application may then choose to handle this case as best they see fit, for instance see: https://github.com/developmentseed/rio-viz/blob/main/rio_viz/templates/index.html#L1202-L1210
Open question:
I'm not sure how this behavior should be defined in more complex situations, where the object is not convex. For instance, a U-shaped polygon that crosses the antimeridian could result in multipolygon with 3 components.
Example:
/bounds?sceneid=S1A_IW_GRDH_1SDV_20230726T183302_20230726T183327_049598_05F6CA_31E7
returns
[-180.0, 61.06949078480844, 180.0, 62.88226850489882]
because it's productInfo.json looks similar to this:
"footprint" : {
"type" : "MultiPolygon",
"coordinates" : [
[ [ [ 180.0, 62.52138872791733 ], [ 179.94979152463551, 62.52658106629492 ], [ 179.46990657679865, 62.57485506167223 ], ..., [ 179.97829453382462, 61.06949078480844 ], [ 180.0, 61.115074335603865 ], [ 180.0, 62.52138872791733 ] ] ],
[ [ [ -180.0, 62.52138872791733 ], [ -180.0, 61.115074335603865 ], [ -179.95528668344312, 61.208976568342585 ], ..., [ -179.34123633603335, 62.451941133348306 ], [ -179.81084543337622, 62.50182719954094 ], [ -180.0, 62.52138872791733 ] ] ] ]
},
😭 there are 2 issues here:
- We create the
boundsin https://github.com/cogeotiff/rio-tiler-pds/blob/main/rio_tiler_pds/sentinel/aws/sentinel1.py#L76 but checking the min/max of the Geometry. In your case, the dataset is crossing the dateline so the geometry is a MultiPolygon (one polygon on each side).
To fix this we can do something like:
def get_bounds(geom: Dict) -> Tuple[float, float, float, float]:
"""Get Bounds from GeoJSON geometry and handle multi polygon crossing the antimeridian line.
ref: https://github.com/cogeotiff/rio-tiler-pds/issues/77
"""
if geom["type"] == "MultiPolygon":
bounds = [
featureBounds({"type": "Polygon", "coordinates": poly})
for poly in geom["coordinates"]
]
minx, miny, maxx, maxy = zip(*bounds)
return [max(minx), min(miny), min(maxx), max(maxy)]
return featureBounds(geom)Which will return [175.38292995222957, 61.06949078480844, -179.34123633603335, 62.88226850489882] for the S1A_IW_GRDH_1SDV_20230726T183302_20230726T183327_049598_05F6CA_31E7 dataset (instead of [-180.0, 61.06949078480844, 180.0, 62.88226850489882])
- There is another issue within rasterio/gdal for the Sentinel-1 dataset. The COG are not directly georeferenced but uses internal GCPS. In rio-tiler we use those GCPS to virtually warp the data to geographic projection and then
read tile. Sadly for GCPS dataset crossing the dateline there is an issue (see: rasterio/rasterio#2892).
Seems like no one has taken any action on your rasterio bug?
Are there any workarounds that can be implemented in rio-tiler-pds while waiting for that to be resolved, or have you already implemented something?
Sadly even if we fix the bounds in rio-tiler-crs (as mentioned in the previous comment) there will still be the issue in rasterio/gdal
@jonaraphael I just publish a new version with a fix for #78 and the multi-polygon bug BUT this won't fix the underlying issue in rasterio: rasterio/rasterio#2892)
Thank you for the progress!
is there anyway you can provide the file in rasterio 2892 without it being on user-pays S3? I'd like to explore this.