ENH: use fsspec to allow read remote single-files

[raybellwaves]

Originally asked on gitter: https://gitter.im/rioxarray/community

I see rioxarray can open remote files using a url path (e.g. https://corteva.github.io/rioxarray/latest/examples/read-locks.html?highlight=url#Reading-without-Locks)

It would be good if rioxarray can work with fsspec objects (https://filesystem-spec.readthedocs.io/en/latest/).

There has been recent work on adding this to xarray (pydata/xarray#4823) and geopandas (geopandas/geopandas#1535)

[snowman2]

Sounds like a good addition. @raybellwaves is this something you wouldn't mind working on?

[snowman2]

I suspect the area you would want to change is around here:

rioxarray/rioxarray/_io.py

Line 715 in 7684848

if isinstance(filename, rasterio.io.DatasetReader):

[scottyhq]

Just wanted to link to some issues and explain why this could be helpful, but also why it might be somewhat difficult. We (mostly @jhamman and @martindurant) looked into this a while back b/c of trying to stream data from google requester pays buckets to xarray datasets and rasterio/gdal didn't have that implemented (rasterio/rasterio#1948) but you could do it via fsspec.open().

So this PR was an attempt to make xr.open_rasterio() accept fsspec objects pydata/xarray#4140

But, it was unclear how to actually implement this because it seemed that ultimately GDAL required that file-like objects be read entirely into MemoryFile(), sort of defeating the purpose of efficient header-only initial access. rasterio/rasterio#977 (comment)

Would definitely be nice though to have the common interface (credentials in particular) that fsspec provides for remote data!

[jhamman]

But, it was unclear how to actually implement this because it seemed that ultimately GDAL required that file-like objects be read entirely into MemoryFile(), sort of defeating the purpose of efficient header-only initial access. mapbox/rasterio#977 (comment)

Exactly. This seems to be the crux to making this work right now. I ran out of steam after @sgillies last comment there. I suspect this is something that can be worked around but I don't know how.

[djhoese]

(Disclaimer: coming in with very little knowledge of previous discussions on this or the internals of rasterio)

It seems like the function @sgillies points to returns a VSILFile which is supposed to be anything that implements this interface: https://gdal.org/api/cpl_cpp.html#classVSIVirtualHandle. So theoretically if rasterio could be provide a wrapper class for that (or some other type of base class) that it knew it could just pass to GDAL then wouldn't that work? Granted, this is mixing a lot of python and C++ class stuff which I'm sure gets complicated pretty fast especially when GDAL is involved.

...but its programming, we can do whatever we want, right?

[sgillies]

@djhoese that's a very intriguing idea and I would be interested in looking into this with you. I admire fsspec and its chaining paradigm. GDAL's own chaining is less flexible, offers no control over caching, and due to under-specification of chaining there are some bugs in format driver implementations.

[djhoese]

@sgillies I would be interested in looking this too. I'm not sure I have the time this month (or even next if I'm being honest), but it is interesting enough that I'll probably make time for it.

Just so I'm clear, what @scottyhq and @jhamman wanted is the ability to provide an open file-like object to rasterio and have rasterio intelligently/selectively load data as it is needed. Is that right @jhamman @scottyhq? The solution that we would be attempting is adding a new low-level class in rasterio to pass to GDAL that would allow random/byte-range access to this file-like object. Am I understanding this correctly?

[scottyhq]

@djhoese and @sgillies thanks so much for chiming in! any time you're willing to explore this would be greatly appreciated.

Just so I'm clear, what @scottyhq and @jhamman wanted is the ability to provide an open file-like object to rasterio and have rasterio intelligently/selectively load data as it is needed.

Yes. I definitely feel out of my depth when it comes to the internals, but I can provide a few code examples that I think help guide things. The first order issue (from PRs linked above if you want to follow the paper trail):

# 1st step, improve efficiency here (Wall time: 3.12 s, should be ~100ms to read header)
gcp_settings = dict(token='cloud', project='pangeo-181919', requester_pays=True)
with fsspec.open('gs://pangeo-ncar-soilgrids/TAXNWRB_250m.tif', **gcp_settings) as fobj:
    with rasterio.open(fobj) as src: 
        print(src.profile)
        # also will want to read data subsets

the above snippet will run via a pangeo binder session on GCP if that is helpful:

Ultimately (and please others chime in if there are other ideas here! don't want to highjack @raybellwaves 's issue), I think the goal would be to align with other xarray backends, in such a way that URIs are parsed and handled behind the scenes by fsspec:

# works in xarray>0.17
uri = "s3://era5-pds/zarr/2020/12/data/eastward_wind_at_10_metres.zarr"
aws_settings = dict(anon=True)
fsspec_options = dict(storage_options=aws_options)
ds = xr.open_mfdataset(uri, 
                       engine="zarr", 
                       backend_kwargs=fsspec_options)

And this would be a similar syntax for rioxarray (related to #197):

# Goal (does not currently work)
uri = "gs://pangeo-ncar-soilgrids/TAXNWRB_250m.tif"
gcp_settings = dict(token='cloud', project='pangeo-181919', requester_pays=True)
fsspec_options = dict(storage_options=gcp_settings)
ds = xr.open_dataset(uri, 
                     engine="rioxarray", 
                     backend_kwargs=fsspec_options)

Hope that helps!

[raybellwaves]

May be of interest here https://github.com/gjoseph92/stackstac#roadmap. "Easier access to s3://-style URIs (right now, you'll need to pass in gdal_env=stackstac.DEFAULT_GDAL_ENV.updated(always=dict(session=rio.session.AWSSession(...))))" cc. @gjoseph92

[gjoseph92]

Thanks @raybellwaves. I think @djhoese's idea of basically writing a translation layer from the Python file interface to the GDAL VSIFile interface is very promising, and would have pretty broad applicability. That would be really nice to have. Kind of reminds me of pyarrow.PythonFile, which could be useful to reference.

cc @TomAugspurger

[djhoese]

I was going to ask if anyone knew of a public S3-stored geotiff that I could test with locally, but I'm realizing that I could probably use any local geotiff opened with fsspec, right? For more benchmarking, if anyone has a public S3 geotiff that doesn't charge anyone to read, let me know.

[martindurant]

Not necessarily - some libraries use os.fspath to get an equivalent local filename when given a file-like.

…

On March 27, 2021 11:07:00 AM EDT, David Hoese ***@***.***> wrote: I was going to ask if anyone knew of a public S3-stored geotiff that I could test with locally, but I'm realizing that I could probably use any local geotiff opened with fsspec, right? For more benchmarking, if anyone has a public S3 geotiff that doesn't charge anyone to read, let me know. -- You are receiving this because you were mentioned. Reply to this email directly or view it on GitHub: #246 (comment)

-- Sent from my Android device with K-9 Mail. Please excuse my brevity.

[gjoseph92]

@djhoese here are a few:

• 343x343 uint8 RGB overview https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/13/S/DV/2020/4/S2B_13SDV_20200428_0_L2A/L2A_PVI.tif
• 1830x1830 uint16 1-band https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/13/S/DV/2020/4/S2B_13SDV_20200428_0_L2A/B01.tif
• 10980x10980 uint16 1-band https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/13/S/DV/2020/4/S2B_13SDV_20200428_0_L2A/B02.tif
• 10980x10980 uint16 RGB https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/13/S/DV/2020/4/S2B_13SDV_20200428_0_L2A/TCI.tif

[djhoese]

I got something that compiles! ...now I just need to make it do what it is supposed to. I'll make a PR as soon as it resembles what it should. From then on we can define the exact interfaces and functionality.