Bbox and Url Map

[rmg55]

@gjoseph92 - really enjoying working with this package FYI!

I had a couple ideas for features that I thought might be worthwhile to discuss.

1. Include an exact: bool argument. This would:
   • Create a list of unique bboxes (snapped) from the STAC object
   • For each unique bbox, open a single cog and read the bounds like:

     with rasterio.open(url) as src:
          bnds = src.bounds

2. For folks that have compute next to the data, I assume there would be an io benefit to remap the raster urls (correct? - I have not tested this yet though). Possibly include a remap_url: func argument. This would be a function that would alter the url/href of the raster allowing to remap the url from /vsicurl to /vsis3 (or other). Alternative would be to build this logic in with (remap_url: str such as s3 or az or gc).

Thoughts? Would be happy to put together a PR if either of these seem like they would fit into stackstac

[gjoseph92]

Thanks @rmg55, glad to hear it's working out for you! Both great suggestions; I'll reply separately for easier threading. If you want to explore either of these in PRs, that would be awesome.

For the exact: bool, what's the use-case/issue you're trying to address here? Sounds to me like you'd use this when you don't fully trust the STAC metadata (or it's missing the proj extension, etc.), but you also don't know the exact output grid parameters you want (so you still want stackstac to pick it for you), so you want to force it to get georeferencing metadata from the data itself instead of STAC?

I think this is one case of a more general issue: wanting/needing to get some metadata out of the files directly, instead of STAC.

Instead of putting this in stack, I'd like to see some sort of patch_from_data or infer_from_data functions, which go look at the data, and then return an updated STAC ItemCollection with the correct metadata filled in.

1. You can save your updated ItemCollection (or even just the JSON patch?) to disk, so you don't have a repeat this slow operation.
2. We keep very explicit boundaries around what will and won't do IO (what's fast & free vs slow & possibly will charge you money).
3. Separate concerns; keep the stack logic focused on one thing.
4. Easier to optimize/parallelize: patch_from_data should probably run itself on Dask, opening datasets in parallel using all the AutoParallelRioReader infrastructure we already have—plus, if your cluster is already near the data, latency (and possibly cost) will be way lower if we can run the IO there. If you want to get the bboxes of 20,000 GeoTIFFs, this will make a huge difference.

So yes, I think this totally fits into stackstac—in fact, I think it's a key part to make it work on real-world data. I'd just like to see it as a separate component, instead of another code path in the stack logic.

[rmg55]

For the exact: bool, what's the use-case/issue you're trying to address here? Sounds to me like you'd use this when you don't fully trust the STAC metadata (or it's missing the proj extension, etc.), but you also don't know the exact output grid parameters you want (so you still want stackstac to pick it for you), so you want to force it to get georeferencing metadata from the data itself instead of STAC?

Yep - mostly for the case when the proj extension is not in the STAC, and using the WGS84 BBOX is not accurate enough to match the raster projected bbox (#5). @matthewhanson probably can speak to this more.

Some use cases are: 1) If you working on ROIs that you want to be sure a pixel is fully encompassed within (hence resampling/warping is a bit problematic). 2) From the scientific perspective, it is awkward to write into the methods (e.g. "we used HLS data, but resampled scenes to a different grid because...) and makes the overall reproducibility a bit more convoluted that just using the original data.

Here an example using Harmonized Landsat Sentinel (HLS) Data from NASA on the LPCloud DAAC that illustrates this.

1. Get STAC (expand for code)

from satsearch import Search
def get_STAC_items(url, collection, dates, bbox):
  results = Search.search(url=url,
                      collections=collection, 
                      datetime=dates,
                      bbox=bbox)

  return(results)

url = 'https://cmr.earthdata.nasa.gov/stac/LPCLOUD/' 
collection = ['HLSS30.v1.5']#'C1711924822-LPCLOUD' #HLS
bbox=[-104.79107047,   40.78311181, -104.67687336,   40.87008987]
dates = '2013-02-01/2021-04-10'

stac_items_hls = get_STAC_items(url,collection,dates,','.join(map(str, bbox))).items()
s_col = stac_items_hls.geojson()['features']

2. Without specifying the bbox

stackstac.stack(s_col,epsg=32613,resolution=30,chunksize=99999,assets=['B01'])

dask.array<fetch_raster_window, shape=(55, 1, 3715, 3717), dtype=float64, chunksize=(1, 1, 3715, 3717), chunktype=numpy.ndarray>
Coordinates:
  * time            (time) object 1594404158627000000 ... 1616004156165000000
    id              (time) <U19 'G1947862363-LPCLOUD' ... 'G2024286982-LPCLOUD'
  * band            (band) <U3 'B01'
  * x               (x) float64 5e+05 5e+05 5e+05 ... 6.114e+05 6.114e+05
  * y               (y) float64 4.601e+06 4.601e+06 ... 4.489e+06 4.489e+06
    start_datetime  (time) <U24 '2020-07-10T18:02:38.627Z' ... '2021-03-17T18...
    end_datetime    (time) <U24 '2020-07-10T18:02:38.627Z' ... '2021-03-17T18...
    eo:cloud_cover  (time) int64 2 0 0 39 1 0 1 1 80 ... 40 0 4 81 0 56 99 62
Attributes:
    spec:        RasterSpec(epsg=32613, bounds=(499950, 4489410, 611460, 4600...
    crs:         epsg:32613
    transform:   | 30.00, 0.00, 499950.00|\n| 0.00,-30.00, 4600860.00|\n| 0.0...
    resolution:  30        

3. With specifying the bbox

with rasterio.open('https://lpdaac.earthdata.nasa.gov/lp-prod-protected/HLSS30.015/HLS.S30.T13TEF.2020274T174141.v1.5.B01.tif') as src:
  b = src.bounds
stackstac.stack(s_col,epsg=32613,resolution=30,chunksize=99999,assets=['B01'],bounds=list(b))

dask.array<fetch_raster_window, shape=(57, 1, 3660, 3660), dtype=float64, chunksize=(1, 1, 3660, 3660), chunktype=numpy.ndarray>
Coordinates:
  * time            (time) object 1594404158627000000 ... 1616608359899000000
    id              (time) <U19 'G1947862363-LPCLOUD' ... 'G2026583234-LPCLOUD'
  * band            (band) <U3 'B01'
  * x               (x) float64 5e+05 5e+05 5.001e+05 ... 6.097e+05 6.098e+05
  * y               (y) float64 4.6e+06 4.6e+06 4.6e+06 ... 4.49e+06 4.49e+06
    start_datetime  (time) <U24 '2020-07-10T18:02:38.627Z' ... '2021-03-24T17...
    end_datetime    (time) <U24 '2020-07-10T18:02:38.627Z' ... '2021-03-24T17...
    eo:cloud_cover  (time) int64 2 0 0 39 1 0 1 1 80 ... 0 4 81 0 56 99 62 35 15
Attributes:
    spec:        RasterSpec(epsg=32613, bounds=(499980, 4490220, 609780, 4600...
    crs:         epsg:32613
    transform:   | 30.00, 0.00, 499980.00|\n| 0.00,-30.00, 4600020.00|\n| 0.0...
    resolution:  30     

Working with this data, I also noticed a strange difference in performance when calling da.mean().compute(). Not sure if this is some sort of request quota limits on the data combined with how GDAL warp_vrt works?

Estimated BBOX from STAC (1945 seconds):

True BBOX (48 seconds):

I like your proposed workflow for this. For many STACs, especially ones with a large temporal component, this should be fairly fast because we would only need to open one file per unique bbox in the catalog (e.g. in the example above). With other non-temporal catalogs (e.g. NAIP), this would certainly be much more onerous, so having the optimized workflow would be good.

[RichardScottOZ]

NASA STAC? Nice.

[RichardScottOZ]

I imagine optimised those = lots of happy people.

[gjoseph92]

For URL rewriting (from https://sentinel-cogs.s3.us-west-2.amazonaws.com/... to s3://sentinel-cogs/...): yes, I think this fits into stackstac... but I'd like to do more research to confirm it's actually beneficial first.

A quick search makes it sound like S3's only API is over HTTP/HTTPS, so I don't know if s3:// URLs have any magic to them that makes them more performant/cheaper than a request to the http:// URL? I have no idea here, I'd love to hear a definitive answer.

Then there's the question of /vsicurl vs /vsis3. Again, I know little here, but I'd want to check that /vsis3 isn't just a wrapper around /vsicurl. Skimming the docs, it sure sounds like it is: between AWS_REGION, AWS_S3_ENDPOINT, and particularly AWS_VIRTUAL_HOSTING, it sounds like it's just generating a URI to pass to /vsicurl.

The main thing /vsis3 does is add auth to the request. So even if there are no performance differences, routing through /vsis3 could still be handy, in case an asset is given as a plain HTTP URL, but actually requires AWS credentials to access. That said, so far, the assets I've seen on AWS given as https:// have all been public, and the private/requester-pays ones have been s3://, but that's obviously not an official standard (I think?).

Stepping back, if the ideas in corteva/rioxarray#246 (comment) around routing fsspec into GDAL come to fruition, maybe we could handle it at the fsspec level instead, and not worry about /vsicurl vs /vsis3.

[rmg55]

I also am unclear on how the fsspec/s3fs/gs versus http impact performance if compute is in the same location as the data. I probably have time in the next week to put together a quick test or perhaps a good question for the Pangeo gitter page?

Stepping back, if the ideas in corteva/rioxarray#246 (comment) around routing fsspec into GDAL come to fruition, maybe we could handle it at the fsspec level instead, and not worry about /vsicurl vs /vsis3.

Agreed!

[RichardScottOZ]

This is what Zac Flamig said: "The S3 API is http/https because it’s REST so I don’t think there’s any advantage to using S3:// vs https from the perspective of performance"

[rmg55]

@RichardScottOZ - thanks for passing along that info. I did some tests with the S2a data on AWS, and /vsicurl had the same io as /vsis3 (https://github.com/rmg55/CloudDAAC_Binders/blob/main/s3_v_http.ipynb). The one edge case would be if NASA throttles egress out of the S3 buckets, then I wonder if we need the extra configs associated with /vsis3 to access the s3 buckets directly when the compute is in the same region - something to ask the NASA DAAC folks about...

[matthewhanson]

If the http URL is actually just the http version of the s3 location then I don't think there's any way to throttle egress so they should always be the same.
What is sometimes done when the bucket is requester pays though is they put http access through Cloudfront where they can throttle egress. This is what is done for Landsat - you can get faster speeds using the requester pays bucket or use their free, but slower http endpoint.

[rmg55]

What is sometimes done when the bucket is requester pays though is they put http access through Cloudfront where they can throttle egress

@matthewhanson thanks for the info. I am under the impression that the NASA Cloud DAAC storage is likely to do the same thing (e.g. throttle the egress). Do you know if that is that correct?

I am attempting to test throughput on the aws landsat data you mentioned. Is this the AWS Open Datasets ( if so, what is the url of the req. pay bucket) or the USGS bucket

I have put together a pangeo binder on us-west-2 (https://github.com/rmg55/CloudDAAC_Binders) to benchmark the throughput from both /vsicurl/ and /vsis3/, but struggling to get the environmental variables correctly specified. @gjoseph92 and @scottyhq any suggestions on this? Should I be using rasterio.env or stackstac.rio_reader.LayeredEnv or system env. variables (os.environ.update(env))?

[gjoseph92]

@rmg55 see rmg55/CloudDAAC_Binders#1, but the short answer is: use stackstac.LayeredEnv here. (It's pretty awkward to use though, I know.)

os.environ.update will work locally, but when using a distributed cluster, your workers workers may be in separate processes (or machines!), so they won't have the same environment variables. The purpose of LayeredEnv is to let you send an environment to all workers. Don't be thrown off by the "layered" part—just think of it as equivalent to rasterio.Env (the actual layering bit is more advanced and mostly used by stackstac internals to optimize some GDAL stuff).

It's also best practice to start from stackstac.DEFAULT_GDAL_ENV and call .updated(...) to add more options, so you can keep the stackstac defaults (unless you have a strong reason to change them).

[rmg55]

@gjoseph92 - thanks for update to stackstac, the pull request to rmg/CloudDAAC_Binders#1, and the useful explanations!

[scottyhq]

thanks for looking into the benchmarking @rmg55 , very interested to see more results from that effort! I agree it would be great to get some more tests against various STAC catalogs. the NASA CMR STAC endpoint is a trove of minimal and often quirky records (but valid). I was just exploring a dataset of COGs which are pre-mosaiced covering all of Greenland (highly compressed b/c mostly nodata), in which case you only need to read the first file to get relevant proj info: https://gist.github.com/scottyhq/c4a4e889b58a0a153dd5fb18bad9f3e8