Skewed CRS when using Sentinel2Indexer

[KaSchr]

Hey Team - love what you are doing at Clay! I am currently evaluating stacchip to understand if it will serve my needs. Here is what I plan to do as part of my data prep routine:

1. Search S-2 collections given an area of observation (AOI)
2. Chip each retrieved item into 256x256 chips, only retain chips that have <0.1 missing data
3. From the remaining chips, find the chips that intersect the original AOI
4. Get the image data for these chips

I am facing an issue with the Sentinel2Indexer class. Specifically, the chips created by the indexer seem to be slightly shifted in their coordinate system and do not intersect with my original AOI or the underlying S2 geotif. Here are two screenshots to highlight the issues:

1. AOI clearly intersects with the underlying S-2 tile (see screenshot below)

2. Chips as generated using the Sentinel2Indexer seem to be shifted and at a different scale

I have made sure that the CRS for all is the same (EPSG:32720 in this case).

Broadly what I am doing is the following:

I use the following AOI:

coords = [[-61.47502802133589,-10.547903668545445],[-61.47502802133589,-10.749253483741057],[-61.17444254239267,-10.749253483741057],[-61.17444254239267,-10.547903668545445],[-61.47502802133589,-10.547903668545445]]

#convert to gdf
polgyon = Polygon(coords)
aoi_gdf = gpd.GeoDataFrame(index=[0], crs='epsg:4326', geometry=[polgyon])
aoi_geometry = aoi_gdf.geometry

1. I search the S-2 collection for that AOI to get the STAC items

def search_sentinel2_collection(start_date, end_date, aoi_geometry,max_cloud=1,max_items=None):
    """
    Search Sentinel 2 data collection for target_date
    and collect results including meta data in a dictionary.
    This function uses the PySTAC client
    """
    client = pystac_client.Client.open("https://earth-search.aws.element84.com/v1")
    collection = "sentinel-2-l2a"
    
    search_results = client.search(
        collections=[collection],
        query = {"eo:cloud_cover":{"lt":max_cloud}},
        intersects=aoi_geometry.to_crs("EPSG:4326").geometry[0].__geo_interface__, 
        datetime=f"{start_date}/{end_date}",max_items=max_items
    )
        
    return search_results

search_results = search_sentinel2_collection(start_date=analysis_start_date, end_date=analysis_end_date, aoi_geometry=aoi_geometry,max_cloud=5,max_items=1000)

items = search_results.item_collection()
items_list = list(items)

2. I index each item using the Sentinel2Indexer and create the index

indexer = Sentinel2Indexer(items_list[0],chip_size=256,chip_max_nodata=.1)
stac_table = indexer.create_index()

4. For each of the resulting chips I get the bbox and collect them in a geodataframe

import geopandas as gpd
def create_geodf_from_stac(stac_table):
    ls=stac_table.to_pylist()
    #align geometry
    for i in ls:
        geo_raw = i['geometry'][0]
        polygon_coords = [(point['x'], point['y']) for point in geo_raw]
        polygon = Polygon(polygon_coords)
        i['geo_shapely'] = polygon
    
    df = pd.DataFrame(ls)
    df.drop(columns=['geometry'],inplace=True)
    
    geodf = gpd.GeoDataFrame(df,geometry="geo_shapely").set_crs(4326)
    return geodf    


gdf = create_geodf_from_stac(stac_table)

5. I plot the AOI, the gdf and an exemplary geotiff

I have tried this for other AOIs and the the issue remains. I’d appreciate some advice on this. Many thanks!

[isaachcamp]

Hi @KaSchr, probably a bit late at this point but I was having the same problem. I think #10 should have fixed it.