Internal Server Error when querying for non-square polygons in Liberia

[jamescw19]

I'm trying to access some data from the io-lulc-annual-v02 collection through the Python SDK. Whenever I use non-square polygons as the intersects input, I get a 500 Internal Server Error response from MSPC. This seems to be the case when running queries over one of my test areas, Liberia (I have not tested in other regions around this), but the search works as expected if I use the same shape polygons in areas like New Zealand or the UK. A minimal example is shown below:

import pystac_client
import planetary_computer

catalog = pystac_client.Client.open(
    "https://planetarycomputer.microsoft.com/api/stac/v1", 
    modifier=planetary_computer.sign_inplace
)

# A triangular Polygon in Liberia
triangle_liberia = {
  "coordinates": [
    [
      [
        -10.638401890350423,
        6.390713138106577
      ],
      [
        -10.638577763341829,
        6.39056202611188
      ],
      [
        -10.63829380174181,
        6.390563846739283
      ],
      [
        -10.638401890350423,
        6.390713138106577
      ]
    ]
  ],
  "type": "Polygon"
}

# A square Polygon in Liberia
square_liberia = {
  "coordinates": [
    [
      [
        -10.638039152306249,
        6.390847864424117
      ],
      [
        -10.638039152306249,
        6.389995810894121
      ],
      [
        -10.637337492351634,
        6.389995810894121
      ],
      [
        -10.637337492351634,
        6.390847864424117
      ],
      [
        -10.638039152306249,
        6.390847864424117
      ]
    ]
  ],
  "type": "Polygon"
}

# A triangular Polygon in New Zealand
triangle_nz = {
  "coordinates": [
    [
      [
        175.6264093174089,
        -38.940290230303006
      ],
      [
        175.62522338888976,
        -38.941362220141755
      ],
      [
        175.627627298052,
        -38.9413372903292
      ],
      [
        175.6264093174089,
        -38.940290230303006
      ]
    ]
  ],
  "type": "Polygon"
}

# Fails
catalog.search(collections=["io-lulc-annual-v02"], intersects=triangle_liberia, max_items=1).item_collection()
# Succeeds
catalog.search(collections=["io-lulc-annual-v02"], intersects=square_liberia, max_items=1).item_collection()
# Succeeds
catalog.search(collections=["io-lulc-annual-v02"], intersects=triangle_nz, max_items=1).item_collection()
# Succeeds
catalog.search(collections=["sentinel-2-l2a"], intersects=triangle_liberia, max_items=1).item_collection()

[777arc]

Very strange, I'm able to reproduce it, nothing in the logs pops out but clearly something is going on, will keep looking

[ghidalgo3]

@jamescw19 I think Marc is out but our internal investigation revealed that the library that handles the geometry processing is truncating decimal numbers from your polygon coordinates and creating polygons with zero area. Can you make your search polygon larger and reduce the number of decimal points to no more than 0.0001 degrees?

Unrelated: whenever I see a WSG84 coordinate with that many decimal points producing creating an error, my first action is to reduce the number of decimal points to see if it still happens with a lower resolution. I found this answer a good guide to know how many decimal points I need for a WSG84 coordinate https://gis.stackexchange.com/questions/8650/measuring-accuracy-of-latitude-and-longitude/8674#8674

[jamescw19]

Thanks for the link, very interesting!
I gave this a try again with the following polygon:

{
    "coordinates": [[
        [-10.7, 6.4], 
        [-10.7, 6.3], 
        [-10.6, 6.3], 
        [-10.7, 6.4]
    ]],
    "type": "Polygon"
}

Same error.

I tried tweaking the latitude values (longitude showed no issues) and found that the failures occurred between [-8.0, 8.0]. Anything above 8.1 and I got items back. Same with setting longitude to 101.x. So there's something up with that band 8 degrees either side of the equator, with this data source

[777arc]

It might be related to this issue https://trac.osgeo.org/postgis/ticket/5583 which was fixed in a newer version of GEOS, so we'll attempt to bump the version of PostGIS we're using and see if it fixes it, so standby.