<div dir="ltr"><div>From what I could tell, you are attempting to re-project a bounding box across the antimeridian.</div><div><br></div><div>Here is one way to handle that (<a href="https://pyproj4.github.io/pyproj/stable/api/transformer.html#pyproj.transformer.Transformer.transform_bounds">https://pyproj4.github.io/pyproj/stable/api/transformer.html#pyproj.transformer.Transformer.transform_bounds</a>):</div><div>>>> import shapely.geometry<br>>>> from pyproj import Transformer<br>>>> transformer = Transformer.from_crs(32601, 4326, always_xy=True)<br>>>> def bounding_polygon(left, bottom, right, top):<br>... if right < left:<br>... return shapely.geometry.MultiPolygon(<br>... [<br>... shapely.geometry.box(left, bottom, 180, top),<br>... shapely.geometry.box(-180, bottom, right, top),<br>... ]<br>... )<br>... return shapely.geometry.box(left, bottom, right, top)<br>... <br>>>> shapely.geometry.mapping(bounding_polygon(*transformer.transform_bounds(377120, 7577600, 418080, 7618560)))<br>{'type': 'MultiPolygon', 'coordinates': [(((180.0, 68.28494231541386), (180.0, 68.66685731991039), (179.97430985556775, 68.66685731991039), (179.97430985556775, 68.28494231541386), (180.0, 68.28494231541386)),), (((-178.98557986959838, 68.28494231541386), (-178.98557986959838, 68.66685731991039), (-180.0, 68.66685731991039), (-180.0, 68.28494231541386), (-178.98557986959838, 68.28494231541386)),)]}</div><div><br></div><div>Alternatively, you could use rasterio.warp.transform_bounds (rasterio 1.3+ and GDAL 3.4+).<br><a href="https://github.com/rasterio/rasterio/issues/2313">https://github.com/rasterio/rasterio/issues/2313</a><br><br></div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Tue, Mar 29, 2022 at 3:30 PM Robin Wilson <<a href="mailto:robin@rtwilson.com">robin@rtwilson.com</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div><div style="font-family:Helvetica,Arial;font-size:13px">Hi,</div><div style="font-family:Helvetica,Arial;font-size:13px"><br></div><div style="font-family:Helvetica,Arial;font-size:13px">I’ve been running into problems creating STAC items for images that cross the antimeridian. The geometries for the images are coming out as spanning the width of the whole world (from -180 to +180 degrees). With the help of the author of the tool I was using to create the STAC items (rio-stac), I’ve narrowed it down to a simple test case of reprojecting geometries (which is what the tool is doing internally).</div><div style="font-family:Helvetica,Arial;font-size:13px"><br></div><div style="font-family:Helvetica,Arial;font-size:13px">If you put the following GeoJSON for an antimeridian-crossing polygon in WGS 84 / UTM zone 1N in a file called test_32601.geojson:</div><div style="font-family:Helvetica,Arial;font-size:13px"><br></div><div style="font-family:Helvetica,Arial;font-size:13px">```</div><div style="font-family:Helvetica,Arial;font-size:13px"><div style="margin:0px">{"type":"Polygon","coordinates":[[[377120,7577600],[418080,7577600],[418080,7618560],[377120,7618560],[377120,7577600]]]}</div><div>```</div><div><br></div><div>and then run the following:</div><div><br></div><div>```</div><div><div>ogr2ogr -f GeoJSON -s_srs epsg:32601 -t_srs epsg:4326 out_4326.geojson test_32601.geojson</div></div><div>```</div><div><br></div><div>The contents of the output file will be:</div><div><br></div><div>```</div><div>{ "type": "Polygon", "coordinates": [ [ [ -179.018099960087909, 68.666857319910392 ], [ -178.985579869598382, 68.299719485120988 ], [ -179.976982644969041, 68.284942315413858 ], [ 179.974309855567753, 68.651801088224246 ], [ 180.0, 68.652259945574556 ], [ 180.0, 68.652259945541758 ], [ 180.0, 68.652259944541754 ], [ -180.0, 68.652259944541754 ], [ -180.0, 68.652259945541758 ], [ -180.0, 68.652259945610126 ], [ -179.018099960087909, 68.666857319910392 ] ] ] }</div><div>```</div><div><br></div><div>In that output geometry there are longitudes of 180 and -180, so when plotted on a map it spans the whole width of the world.</div><div><br></div><div>The author of the rio-stac library did a PR to try and fix this problem (<a href="https://github.com/developmentseed/rio-stac/pull/30" target="_blank">https://github.com/developmentseed/rio-stac/pull/30</a>), but it doesn’t seem to have fixed it for this example, as the underlying GDAL code used by his library (through rasterio, I believe) is giving this result.</div><div><br></div><div>I find this result unexpected, but I have no knowledge of how GDAL deals with the antimeridian. Is this to be expected, or is this potentially a bug? If it isn’t a bug, is there any way around this?</div><div><br></div><div><div>Best regards,</div><div><br></div><div>Robin</div></div><div><br></div><div>Dr Robin Wilson</div><div><a href="http://www.rtwilson.com" target="_blank">www.rtwilson.com</a></div><div><br></div></div><br><div></div></div>_______________________________________________<br>
gdal-dev mailing list<br>
<a href="mailto:gdal-dev@lists.osgeo.org" target="_blank">gdal-dev@lists.osgeo.org</a><br>
<a href="https://lists.osgeo.org/mailman/listinfo/gdal-dev" rel="noreferrer" target="_blank">https://lists.osgeo.org/mailman/listinfo/gdal-dev</a><br>
</blockquote></div><br clear="all"><br>-- <br><div dir="ltr" class="gmail_signature"><div dir="ltr"><div><div dir="ltr"><div><div dir="ltr"><div><div dir="ltr"><font style="font-family:arial,helvetica,sans-serif">Alan Snow</font><br></div></div></div></div></div></div></div></div>