Problematic stackstac requirement of deprecated proj:epsg attribute

[MathewNWSH]

Hi, While using stackstac, I noticed that it is not compatible with the proj extension 2.0. Specifically, based on the knowledge from the extension's repository: the proj:epsg attribute has been deprecated and replaced with proj:code.

This is problematic because the deprecated proj:epsg is still being used by Stackstac.

Attached is an example of a catalog utilizing the proj extension 2.0: https://stac.dataspace.copernicus.eu/v1/collections/sentinel-2-l2a/items/S2A_MSIL2A_20241107T113311_N0511_R080_T31VDH_20241107T123948

ss you can see, the projection is provided via proj:code but remains inaccessible to StackSTAC, which still looks for the hardcoded proj:epsg.

example stackstac output:

[MarcYin]

Extracting and specifying the epsg code will cause inconsistency in the output array shape, if latest pystac (pystac==1.13) is used.

Bellow is an example taken from: https://stackstac.readthedocs.io/en/latest/basic.html

def extract_epsg(item):
    """
    Extract EPSG code from STAC item properties.
    Returns an integer EPSG code.
    Raises NotImplementedError if proj:code uses a non-EPSG authority.
    Raises ValueError if no code found or invalid format.
    Related to github issue:
    # https://github.com/gjoseph92/stackstac/issues/262
    
    !!! Do not use this function, as specifing epsg can cause in consistencies in the output
    array shape to the one in the STAC item!!!
    

    """
    proj_code = item.properties.get('proj:code')
    if proj_code:
        # Example proj_code: "EPSG:32630"
        try:
            authority, code = proj_code.split(':', 1)
            if authority.upper() != 'EPSG':
                raise NotImplementedError(
                    f"Unsupported projection authority '{authority}' in proj:code. "
                    "Only 'EPSG' is implemented."
                )
            epsg_code = int(code)
            return epsg_code
        except (ValueError, AttributeError):
            raise ValueError(f"Invalid proj:code format: {proj_code}")
    
    epsg_code = item.properties.get('proj:epsg')
    if epsg_code:
        if isinstance(epsg_code, int):
            return epsg_code
        else:
            raise ValueError(f"Invalid proj:epsg value: {epsg_code}")
    
    raise ValueError("EPSG code not found in item properties (checked proj:code and proj:epsg).")


import stackstac
import pystac_client

lon, lat = -105.78, 35.79

URL = "https://earth-search.aws.element84.com/v1"
catalog = pystac_client.Client.open(URL)
items = catalog.search(
    intersects=dict(type="Point", coordinates=[lon, lat]),
    collections=["sentinel-2-l2a"],
    datetime="2020-03-01/2020-06-01"
).item_collection()

# to filter only use N0500 processing baseline product
items = [i for i in items if 'N0500' in i.properties.get('s2:product_uri')]
print(len(items))

epsg = extract_epsg(items[0])
print(f"EPSG code: {epsg}")

stack = stackstac.stack(items, epsg=epsg, resolution=10)
stack

Expected to be 10980x10980 for array in x and y direction.

If maintaining the exact DataArray shape is critical, it’s recommended to pin pystac to version 1.11.0, as done in this commit: https://github.com/stac-utils/xpystac/commit/466152fe9911b04c08bb648b6dda21943a763280.

Interestingly, the shape inconsistency occurs even when the correct EPSG code is specified. This suggests that the issue may come from how the output bounds are calculated, potentially related to the code here: https://github.com/gjoseph92/stackstac/blob/38571909cd83d1f3da0c20e44450560e174ff845/stackstac/prepare.py#L208