Can't build from static STAC Catalog which has local file paths

[maawoo]

Thanks for reaching out in Jena-Earth-Observation-School/sdc-tools#43, @print-sid8.

The rasteret version I am running is 0.3.12. As the title says, we are using static STAC Catalogs on a fileserver connected to our university HPC system. The data is mirrored from Digital Earth Africa (see here for their metadata explorer of the s2_l2a product and here for an example STAC Item). The Items of our mirrored dataset differ only by using a relative path for each asset's href but seems like this should not be an issue?

First of all, if I run the following:

collection = rasteret.build_from_stac(
    name="sdc-s2-l2a",
    stac_api="file:///geonfs/02_vol3/SaldiDataCube/original_data/S2_L2A/catalog.json",
    collection="S2_L2A_catalog",
    static_catalog=True,
)

It takes 18 minutes to run into this error (not sure what happens in this time, but maybe the error could be raised earlier?):

ValueError: COG header enrichment produced no band metadata. Verify your band_map points to tiled COG assets and that the STAC items include those assets.

When adding a basic band_map dict

band_map = {
    "B01": "B01",
    "B02": "B02",
    "B03": "B03",
    "B04": "B04",
    "B05": "B05",
    "B06": "B06",
    "B07": "B07",
    "B08": "B08",
    "B8A": "B8A",
    "B09": "B09",
    "B11": "B11",
    "B12": "B12",
    "SCL": "SCL",
}

collection = rasteret.build_from_stac(
    name="sdc-s2-l2a",
    stac_api="file:///geonfs/02_vol3/SaldiDataCube/original_data/S2_L2A/catalog.json",
    collection="S2_L2A_catalog",
    static_catalog=True,
    band_map=band_map,
)

I am getting errors that rust panicked:

See error_log1.txt for the full error log.

When changing the stac_api path to "/geonfs/02_vol3/SaldiDataCube/original_data/S2_L2A/catalog.json", the error changes to the following:

OSError: Failed to fetch bytes from /geonfs/02_vol3/SaldiDataCube/original_data/S2_L2A/34JCN/2018/1/S2A_34JCN_20180106_0_L2A/B01.tif: relative URL without a base

See error_log2.txt for the full error log.

Not sure if I should do something differently or if this is a bug in rasteret?

[print-sid8]

Thanks a lot for raising this! and providing such detailed info.

This will require a deeper look.

The rust panics are due to the boundary of rasteret to Obstore (which is a rust aiohttp + multicloud auth handling library).

My initial suspicions by just looking at logs and errors, is either a bug in rasteret, or obstore, or the file itself is slightly corrupt/different than expected.

You can clarify the third guess, have you already read that exact URL in some other tool?

[maawoo]

Each error log complains about a different Item and GeoTIFF asset. I can confirm that both Items and GeoTIFF files are fine / not corrupted. Thanks for having a look into this!

[print-sid8]

@maawoo can you try out the new branch in #28 PR ?

im yet to merge it, im reviewing it, would be good if you can manually build rasteret from this branch at your end.

This time just do the initial build_from_stac command, no need to try band_map unless it gives you some errors.

It shudnt give any errors, nor should it do that long wait.

collection = rasteret.build_from_stac(
    name="sdc-s2-l2a",
    stac_api="file:///geonfs/02_vol3/SaldiDataCube/original_data/S2_L2A/catalog.json",
    collection="S2_L2A_catalog",
    static_catalog=True,
)

[print-sid8]

Just a reminder @maawoo that if you don't provide any filters to build steps, it should try to build/cache every available tif file in the stac catalog. So that will definitely take time. Best to test with some date and geometry filters.

[maawoo]

I tested the new branch on a subset (1 year of a single tile) and the build step worked. However, given that it took around 24 min for ~1/184 of our Sentinel-2 data, I'm not sure if our use case actually fits to rasteret's goals?

I wanted to test loading to Xarray now but it failed with:

ValueError: Record CRS is missing (proj:epsg) so Rasteret cannot transform the WGS84 query geometry into raster CRS. Fix by adding a per-record `proj:epsg` column to your record table, or build with `enrich_cog=True` so Rasteret can infer CRS from GeoTIFF headers when available.

Wanted to rebuild and test again but rasteret.build_from_stac doesn't include an enrich_cog parameter.

Each STAC Item includes proj:code, which should be prioritized. Some older STAC Items might include proj:epsg but this field was deprecated with v1.2.0: https://github.com/stac-extensions/projection

I think rasteret should have already raised an error during the building step when it failed to get CRS information.

[print-sid8]

Thanks for trying it out.

Wanted to rebuild and test again but rasteret.build_from_stac doesn't include an enrich_cog parameter.

build_from_stac doesn't need that param. Rasteret assumes you are sending STAC with COGs in it and will naturally index them.

Each STAC Item includes proj:code

This has been on my roadmap so this will be fixed soon. build_from_table() already captures newer proj keys

However, given that it took around 24 min for ~1/184 of our Sentinel-2 data, I'm not sure if our use case actually fits to rasteret's goals?

Did you mean 1 file out of 184? Can u explain that figure?

Also, the indexing will take time but not 24 minutes for 1 file thats not normal. And rasteret's read speed comes after build phase, when you use the get_xarray/get_numpy and other functions on the Collection you just built.

i will be pushing fixes to same branch in sometime

[print-sid8]

Some patches for this are pushed to same PR branch. Do try again when you have some free time thanks.

This debugging would have been better if I was able to directly access the STAC myself and find all edge cases, but I don't think I can. So thanks for testing again by yourself @maawoo

[maawoo]

build_from_stac doesn't need that param.

Just mentioned that because the error message suggested using it to fix the missing CRS 😉

Did you mean 1 file out of 184? Can u explain that figure?

We have stored 23 tiles (e.g., 34JCN), each containing 8 years of data. So the subset I used was roughly 1/184th of the entire Sentinel-2 data we have stored on our file servers.

I understand that rasteret will result in a performance boost after the build phase, but seems like this will take a while. Correct me if I'm wrong, but it seems like rasteret is made for speeding up workflows that repeatedly access small spatio-temporal subsets of a larger dataset. It is not made for speeding up access to the entire dataset. That's why I came to the conclusion that our use case doesn't really fit rasteret's goals.

[print-sid8]

For large/full data scans it should provide speeds up too.

PyTorch, TorchGeo, Numpy, Geodataframe, Xarray, should all work faster than before.

Xarray risks OOM coz it doesn't stream/iterate, but i guess you are in HPC cluster, you would know your workloads better.

I will be glad to hear about your experience for full scans. Coz it will help me know workloads like yours better, ive got more feedback from folks doing AI ML with gridded/tiled/random access on large country wide rasters, or point sampling of 1000s of points to get timeseries graphs.

For 23 UTM tiles for 8 years of sentinel 2 over africa im sure that's lots of sentinel 2 scenes. So yes build phase can be slow even for subset, maybe a few 1000 scenes are part of your subset? Rasteret does show in logs the progress of build . would love to know the amount of scenes you dealing with in subset.

P.S. I will be travelling for a while now, will be back online in 2 to 3 hours.
Look forward to continuing discussing here

[print-sid8]

Hey we could attribute the long time taken to also the JSON traversing/crawling. Maybe the catalog JSON contains too many keys? I could check that. A full log of build phase from your side will be more useful.