Read rasters, indexes and tiling schemas from s3:// as well as disk

pointCollection/geoIndex.py

@@ -190,17 +190,25 @@ def from_list(self, index_list, dir_root=''):
         self.attrs['n_files']=len(fileListTo)
         return self
 
-    def from_file(self, index_file, read_file=False, group='index'):
+    def from_file(self, index_file, read_file=False, group='index', fs=None):
         """
         read geoIndex info from file 'index_file.'
         If read_file is set to False, the file is not read, but the
         h5_file_index attribute of the resulting geoIndex is set to a
         reference to the hdf_file's 'index' attribute.  This seems to be
         faster than reading the whole file.
+
+        index_file may be a local path or a URI (e.g. s3://bucket/key).
+        fs is the filesystem used for a remote index; if None, a session from
+        the default AWS credential chain is used.  Note that an index lives in
+        our own bucket even when the granules it indexes are DAAC holdings, so
+        this is deliberately not the granules' earthaccess session.
         """
 
-        import h5py
-        h5_f = h5py.File(os.path.expanduser(index_file),'r')
+        h5_f = pc.io_utils.open_h5(os.path.expanduser(index_file), 'r',
+                                   fs=fs or (pc.io_utils.get_s3fs(daac=None)
+                                             if pc.io_utils.is_remote_path(index_file)
+                                             else None))
         h5_i = h5_f[group]
         if read_file:
             for bin in h5_i.keys():

pointCollection/grid/data.py

@@ -540,7 +540,10 @@ def from_geotif(self, file, date_format=None, **kwargs):
                 self.time=np.nan
                 return self
         try:
-            ds=gdal.Open(file, gdalconst.GA_ReadOnly)
+            # s3://bucket/key -> /vsis3/bucket/key; local paths pass through.
+            # GDAL's /vsis3 uses the ordinary AWS credential chain, so a raster
+            # in our own bucket needs no earthaccess session.
+            ds=gdal.Open(pc.io_utils.as_gdal_path(file), gdalconst.GA_ReadOnly)
             self.from_gdal(ds, **kwargs)
         except Exception as e:
             if 'verbose' in kwargs and kwargs['verbose']:
@@ -682,14 +685,14 @@ def from_gdal(self, ds, field='z', bands=None, bounds=None, extent=None,
         self.__update_size_and_shape__()
         return self
 
-    def h5_open(self, h5_file, mode='r', compression=None):
+    def h5_open(self, h5_file, mode='r', compression=None, fs=None):
         """
         Open an HDF5 file with or without external compression.
 
         Parameters
         ----------
         h5_file: str
-            HDF5 file
+            HDF5 file.  May be a local path or a URI (e.g. s3://bucket/key).
         mode: str, default 'r'
             Mode of opening the HDF5 file
         compression: str or NoneType, default None
@@ -698,20 +701,34 @@ def h5_open(self, h5_file, mode='r', compression=None):
             - ``None``: file is not externally compressed
             - ``'bzip'``
             - ``'gzip'``
+        fs: s3fs.S3FileSystem or NoneType, default None
+            Filesystem used if h5_file is remote.  If None, a session built
+            from the default AWS credential chain is used -- gridded rasters
+            are ancillary data we own, not DAAC holdings, so an earthaccess
+            session would not grant access to them.  Pass an explicit fs to
+            read a raster out of a DAAC bucket.
         """
         # lazy import of h5py
         import h5py
+        # a remote file is opened as a file object and handed to h5py/bz2/gzip,
+        # all three of which accept one in place of a name
+        if pc.io_utils.is_remote_path(h5_file):
+            if fs is None:
+                fs = pc.io_utils.get_s3fs(daac=None)
+            source = fs.open(h5_file, 'rb')
+        else:
+            source = h5_file
         if (compression is None):
-            return h5py.File(h5_file,mode=mode)
+            return h5py.File(source,mode=mode)
         elif (compression == 'bzip'):
             # read bytes from bzip compressed file
-            with bz2.BZ2File(h5_file) as fd:
+            with bz2.BZ2File(source) as fd:
                 fid = io.BytesIO(fd.read())
                 fid.seek(0)
                 return h5py.File(fid, 'r')
         elif (compression == 'gzip'):
             # read gzip compressed file and extract into in-memory file object
-            with gzip.open(h5_file,'r') as fd:
+            with gzip.open(source,'r') as fd:
                 fid = io.BytesIO(fd.read())
                 fid.seek(0)
                 return h5py.File(fid, 'r')
@@ -863,7 +880,7 @@ def from_h5(self, h5_file, field_mapping=None, group='/', \
         bounds=None,  skip=1, fill_value=None,
         t_axis=None, t_range=None, bands=None,
         compression=None, swap_xy=False, source_fillvalue=None,
-        coord_mapping=None):
+        coord_mapping=None, fs=None):
         """
         Read a raster from an HDF5 file.
 
@@ -903,6 +920,9 @@ def from_h5(self, h5_file, field_mapping=None, group='/', \
             - ``None``: file is not externally compressed
             - ``'bzip'``
             - ``'gzip'``
+        fs: s3fs.S3FileSystem or NoneType, default None
+            Filesystem to use if h5_file is a URI rather than a local path.
+            See h5_open().
         swap_xy: bool, default False
             Swap the orientation of x and y variables in the grid
 
@@ -934,7 +954,7 @@ def from_h5(self, h5_file, field_mapping=None, group='/', \
 
         #default
         yorient=1
-        with self.h5_open(h5_file, mode='r', compression=compression) as h5f:
+        with self.h5_open(h5_file, mode='r', compression=compression, fs=fs) as h5f:
             x = np.array(h5f[group][xname]).ravel()
             y = np.array(h5f[group][yname]).ravel()
             for time_var_name in set(['time','t', timename]):
@@ -1089,14 +1109,14 @@ def from_h5(self, h5_file, field_mapping=None, group='/', \
         self.__update_size_and_shape__()
         return self
 
-    def nc_open(self, nc_file, mode='r', compression=None):
+    def nc_open(self, nc_file, mode='r', compression=None, fs=None):
         """
         Open a netCDF4 file with or without external compression.
 
         Parameters
         ----------
         nc_file: str
-            netCDF4 file
+            netCDF4 file.  May be a local path or a URI (e.g. s3://bucket/key).
         mode: str, default 'r'
             Mode of opening the netCDF4 file
         compression, str or NoneType, default None
@@ -1105,8 +1125,25 @@ def nc_open(self, nc_file, mode='r', compression=None):
             - ``None``: file is not externally compressed
             - ``'bzip'``
             - ``'gzip'``
+        fs: s3fs.S3FileSystem or NoneType, default None
+            Filesystem used if nc_file is remote.  If None, a session built
+            from the default AWS credential chain is used; see h5_open().
         """
         import netCDF4
+        if pc.io_utils.is_remote_path(nc_file):
+            # netCDF4.Dataset takes a name or an in-memory buffer, but not a
+            # file object, so a remote file has to be read whole.  That is the
+            # same thing the compression branches below already do, and it is
+            # why remote reads suit ancillary grids rather than large rasters.
+            if fs is None:
+                fs = pc.io_utils.get_s3fs(daac=None)
+            with fs.open(nc_file, 'rb') as fd:
+                buffer = fd.read()
+            if (compression == 'bzip'):
+                buffer = bz2.decompress(buffer)
+            elif (compression == 'gzip'):
+                buffer = gzip.decompress(buffer)
+            return netCDF4.Dataset(uuid.uuid4().hex, mode=mode, memory=buffer)
         if (compression is None):
             return netCDF4.Dataset(nc_file, mode=mode)
         elif (compression == 'bzip'):
@@ -1125,7 +1162,7 @@ def from_nc(self, nc_file, field_mapping=None, group='',
         bands=None, skip=1,
         fill_value=None,
         meta_only=False,
-        t_axis=None, compression=None):
+        t_axis=None, compression=None, fs=None):
         """
         Read a raster from a netCDF4 file.
 
@@ -1167,6 +1204,9 @@ def from_nc(self, nc_file, field_mapping=None, group='',
             - ``None``: file is not externally compressed
             - ``'bzip'``
             - ``'gzip'``
+        fs: s3fs.S3FileSystem or NoneType, default None
+            Filesystem to use if nc_file is a URI rather than a local path.
+            See nc_open().
 
         Returns
         -------
@@ -1192,7 +1232,7 @@ def from_nc(self, nc_file, field_mapping=None, group='',
         dims=[xname, yname, 't', 'time'] + list(self.coordinates or [])
         t=None
         grid_mapping_name = None
-        with self.nc_open(nc_file,mode='r',compression=compression) as fileID:
+        with self.nc_open(nc_file,mode='r',compression=compression, fs=fs) as fileID:
             # set automasking
             fileID.set_auto_mask(False)
             # check if reading from root group or sub-group

pointCollection/io_utils.py

@@ -60,16 +60,57 @@ def strip_pair_suffix(filename):
 
 def get_s3fs(daac='NSIDC', **kwargs):
     """
-    Return a cached, authenticated s3fs.S3FileSystem for the given DAAC,
-    created via earthaccess.get_s3fs_session(). Sessions are cached by
-    (daac, kwargs) so repeated calls don't re-derive credentials.
+    Return a cached s3fs.S3FileSystem.
+
+    Parameters
+    ----------
+    daac : str or None, default 'NSIDC'
+        If a DAAC name, the session is created with
+        earthaccess.get_s3fs_session(), which supplies the short-lived
+        in-region credentials that DAAC's cloud buckets require.
+        If None, an ordinary s3fs.S3FileSystem() is returned, which picks up
+        whatever the default AWS credential chain provides (environment,
+        ~/.aws, or an instance/task role).  That is the right choice for
+        buckets we own rather than read from a DAAC -- ancillary rasters,
+        masks and tiling schemas on s3://maap-ops-workspace/... -- since
+        earthaccess credentials do not grant access to them.
+
+    Sessions are cached by (daac, kwargs) so repeated calls don't re-derive
+    credentials.
     """
     key = (daac, tuple(sorted(kwargs.items())))
     if key not in _S3FS_CACHE:
-        import earthaccess
-        _S3FS_CACHE[key] = earthaccess.get_s3fs_session(daac=daac, **kwargs)
+        if daac is None:
+            import s3fs
+            _S3FS_CACHE[key] = s3fs.S3FileSystem(**kwargs)
+        else:
+            import earthaccess
+            _S3FS_CACHE[key] = earthaccess.get_s3fs_session(daac=daac, **kwargs)
     return _S3FS_CACHE[key]
 
+def as_gdal_path(filename):
+    """
+    Translate a URI into the /vsi... path GDAL uses for the same object.
+
+    GDAL cannot open an 's3://bucket/key' URI directly, but it reads the same
+    object through its /vsis3/ virtual filesystem, which uses the ordinary AWS
+    credential chain (AWS_* environment variables or ~/.aws).  Local paths and
+    paths that are already /vsi... are returned unchanged, so callers can pass
+    everything through this on the way to gdal.Open().
+    """
+    if not is_remote_path(filename):
+        return filename
+    scheme, _, rest = filename.partition('://')
+    scheme = scheme.lower()
+    if scheme == 's3':
+        return '/vsis3/' + rest
+    if scheme == 'gs':
+        return '/vsigs/' + rest
+    if scheme in ('http', 'https', 'ftp'):
+        return '/vsicurl/' + filename
+    # unknown scheme: hand it to GDAL as-is and let GDAL report the problem
+    return filename
+
 def path_exists(filename, fs=None, assume_remote_exists=True):
     """
     Check whether a local or remote file exists.

pointCollection/scripts/query_ATL11_cloud.py

@@ -69,7 +69,8 @@ def index_path_for_granule(granule_basename, index_root):
     return os.path.join(index_root, subdir, granule_basename)
 
 
-def read_ATL11_granule_cloud_items(s3_url, index_file, xr, yr, fields=None, fs=None, version_mismatch='error'):
+def read_ATL11_granule_cloud_items(s3_url, index_file, xr, yr, fields=None, fs=None,
+                                   index_fs=None, version_mismatch='error'):
     """
     Find the rows of a cloud ATL11 granule falling within [xr, yr], using a
     pre-built per-granule geoIndex to locate them, and return them as the
@@ -91,6 +92,11 @@ def read_ATL11_granule_cloud_items(s3_url, index_file, xr, yr, fields=None, fs=N
     fields : list or dict, optional
     fs : s3fs.S3FileSystem, optional
         reused across calls to avoid re-deriving S3 credentials per granule.
+        This is the session for the GRANULE, i.e. a DAAC (earthaccess) one.
+    index_fs : s3fs.S3FileSystem, optional
+        session for a remote index_file.  The index is ours and lives in our
+        own bucket, so it needs the default AWS credentials rather than the
+        DAAC session in `fs`; if None, one is derived on demand.
     version_mismatch : {'error', 'skip'}, optional
         what to do when the granule indexed by index_file does not match
         the cloud-found granule s3_url (e.g. the index was built from an
@@ -111,12 +117,17 @@ def read_ATL11_granule_cloud_items(s3_url, index_file, xr, yr, fields=None, fs=N
     if version_mismatch not in ('error', 'skip'):
         raise ValueError(f"version_mismatch must be 'error' or 'skip', got {version_mismatch!r}")
 
-    if not os.path.isfile(index_file):
+    if pc.io_utils.is_remote_path(index_file) and index_fs is None:
+        index_fs = pc.io_utils.get_s3fs(daac=None)
+    # os.path.isfile() is False for any URI, so a remote index has to be
+    # checked against the bucket -- otherwise every granule would be reported
+    # missing and skipped, and the tile would come back empty rather than failing
+    if not pc.io_utils.path_exists(index_file, fs=index_fs, assume_remote_exists=False):
         warnings.warn(f'query_ATL11_cloud: missing geoIndex {index_file} for granule {s3_url}, skipping')
         return None
 
     s3_basename = os.path.basename(s3_url)
-    gI = pc.geoIndex().from_file(index_file)
+    gI = pc.geoIndex().from_file(index_file, fs=index_fs)
     indexed_name = pc.io_utils.strip_pair_suffix(gI.attrs.get('file_0'))
     if indexed_name is not None:
         indexed_basename = os.path.basename(indexed_name)
@@ -134,7 +145,8 @@ def read_ATL11_granule_cloud_items(s3_url, index_file, xr, yr, fields=None, fs=N
     return D
 
 
-def read_ATL11_granule_cloud(s3_url, index_file, xr, yr, fields=None, fs=None, version_mismatch='error'):
+def read_ATL11_granule_cloud(s3_url, index_file, xr, yr, fields=None, fs=None,
+                             index_fs=None, version_mismatch='error'):
     """
     Read the rows of a cloud ATL11 granule falling within [xr, yr], using a
     pre-built per-granule geoIndex to locate them, concatenated into a
@@ -147,7 +159,8 @@ def read_ATL11_granule_cloud(s3_url, index_file, xr, yr, fields=None, fs=None, v
     pointCollection.data, or None if the granule was skipped
     """
     D = read_ATL11_granule_cloud_items(s3_url, index_file, xr, yr, fields=fields,
-                                        fs=fs, version_mismatch=version_mismatch)
+                                        fs=fs, index_fs=index_fs,
+                                        version_mismatch=version_mismatch)
     if D is None:
         return None
     return pc.data().from_list(D)

pointCollection/tilingSchema.py

@@ -217,16 +217,33 @@ def to_json(self, json_file):
         with open(json_file,'w') as fh:
             json.dump(scheme_dict, fh, indent=2)
 
-    def from_file(self, scheme_file):
+    def from_file(self, scheme_file, fs=None):
+        """
+        Read a tiling schema from a .json or .h5 file.
+
+        Parameters
+        ----------
+        scheme_file : str
+            Schema file.  May be a local path or a URI (e.g. s3://bucket/key),
+            so a schema can live beside the tiles it describes in a bucket.
+        fs : s3fs.S3FileSystem or NoneType, default None
+            Filesystem to use if scheme_file is remote.  If None, a session
+            built from the default AWS credential chain is used: a schema is
+            ours to write, so it sits in our own bucket rather than a DAAC's,
+            even when the tiles it points at are DAAC holdings.
+        """
+        remote = pc.io_utils.is_remote_path(scheme_file)
+        if remote and fs is None:
+            fs = pc.io_utils.get_s3fs(daac=None)
 
         # choose what kind of file this is:
         if scheme_file.endswith('.json'):
-            with open(scheme_file,'r') as fh:
+            with (fs.open(scheme_file, 'r') if remote
+                  else open(scheme_file, 'r')) as fh:
                 scheme_dict = json.load(fh)
         elif scheme_file.endswith('.h5'):
-            import h5py
             scheme_dict={}
-            with h5py.File(scheme_file,'r') as fh:
+            with pc.io_utils.open_h5(scheme_file, 'r', fs=fs) as fh:
                 if 'tiling_schema' in fh:
                     group='tiling_schema'
                 else: