From 72ec08b4008aa8d74cdc9c485ac42f93455dbfa3 Mon Sep 17 00:00:00 2001 From: Ben Smith Date: Wed, 30 Sep 2026 22:22:55 +0000 Subject: [PATCH 1/2] mosaic: read tiles in place from S3, in a process pool; fix add_to_band For mosaicking tiles that live on S3 (ATL1415 on MAAP DPS): - io_utils.glob_remote(pattern): the remote glob.glob, returning sorted URIs. - make_mosaic.py: --directory may be a URI, listed with glob_remote and read in place; --output must then be a local absolute path (h5py writes it). New --block_size (default for remote tiles: DEFAULT_REMOTE_BLOCK_SIZE) and -j/--workers. - grid.mosaic.from_list(block_size=None, workers=1): block_size reaches from_h5/from_nc for remote files -- without it s3fs's 50 MiB default reads most of each tile for a few fields. workers > 1 reads the files in a process pool (forkserver: fork after s3fs has started fails, "This class is not fork-safe"; threads give nothing, h5py's lock serializes the reads). Tiles are still added in list order, at most 2 x workers ahead, so the mosaic is identical to the serial one. workers=1 with no block_size is the original code path. Measured on 557 ATL1415 GL tiles on S3 (MAAP ADE): 40 km average 357 s serial -> 44.5 s with 8 workers; z0 (weighted) 802 s -> 195 s with 4 workers; both bit-identical to the serial mosaics. Each worker costs ~0.35 GiB resident. Fix: add_to_band sliced an in-memory mosaic with item[:,:,in_band], which returns a plain grid.data (grid.data.__copy__), then failed on update_spacing; it now converts back with mosaic().from_grid(), as the grid.data branch beside it already did. Hit by from_list(by_band=True) with any 3-D in-memory mosaic input. tests/test_mosaic_list.py: every read option against the original path (weighted by band / all bands / 2-D, replace, bands, all fields, unreadable tile, in-memory inputs), remote reads through fsspec's memory filesystem, block_size reaching open_remote, make_mosaic.py on a remote directory and its output check, and the add_to_band regression. Co-Authored-By: Claude Opus 5.5 --- pointCollection/grid/mosaic.py | 275 +++++++++++++++++++++---- pointCollection/io_utils.py | 26 +++ pointCollection/scripts/make_mosaic.py | 33 ++- tests/test_mosaic_list.py | 179 ++++++++++++++++ 4 files changed, 468 insertions(+), 45 deletions(-) create mode 100644 tests/test_mosaic_list.py diff --git a/pointCollection/grid/mosaic.py b/pointCollection/grid/mosaic.py index f77c1ba..0c275db 100644 --- a/pointCollection/grid/mosaic.py +++ b/pointCollection/grid/mosaic.py @@ -5,16 +5,145 @@ Routines for creating a weighted mosaic from a series of tiles UPDATE HISTORY: + Updated 09/2026: from_list reads remote files with a block_size, and in a + process pool (workers); tiles are still added in list order Updated 06/2023: calculate x and y arrays using np.arange and spacing updated 03/2021: change scheme for calculating weights, raised cosine as default Updated 03/2020: check number of dimensions of z if only a single band Written 03/2020 """ +import collections +import itertools +import os +import pickle import numpy as np from .data import data import pointCollection as pc +# extensions whose readers (from_h5, from_nc) take a block_size for a remote file +_BLOCK_SIZE_FORMATS = ('h5', 'hdf', 'hdf5', 'nc', 'netcdf') + +def _read_kwargs(item, block_size, **kwargs): + """ + from_file keyword arguments for one item: block_size is added only for a + remote HDF5 or netCDF file, the readers that take it (a local file has no + blocks, and from_geotif has no such argument). + """ + if block_size is not None and pc.io_utils.is_remote_path(item) and \ + os.path.splitext(item)[1][1:].lower() in _BLOCK_SIZE_FORMATS: + kwargs['block_size'] = block_size + return kwargs + +def _read_item(job): + """ + Read one file for a mosaic: (item, meta_only, kwargs) -> (grid, exception). + + Module-level so a process pool can pickle it. An exception is returned, + not raised, so the caller can raise it where the serial read would have + -- inside the same try/except, with the same consequence for the tile. + """ + item, meta_only, kwargs = job + try: + if meta_only: + return pc.grid.data().from_file(item, meta_only=True, **kwargs), None + return pc.grid.mosaic().from_file(item, **kwargs), None + except Exception as exc: + try: + pickle.dumps(exc) + except Exception: + # an exception that cannot cross the process boundary would fail + # the whole pool; carry its text instead + exc = RuntimeError(f'{type(exc).__name__}: {exc}') + return None, exc + +# How a reading pool starts its workers. NOT plain fork: by the time a mosaic +# reads remote tiles the parent has an s3fs event-loop thread (listing the +# tiles starts it), and a child forked from a multi-threaded process can +# deadlock -- s3fs itself refuses, "This class is not fork-safe". forkserver +# forks each worker from a clean single-threaded server that has imported +# pointCollection once, so a worker still starts fast (40 GL tiles, 8 +# workers: 9 s, after a one-off ~10 s server start; serial 32 s). +_START_METHOD = 'forkserver' + +def _init_worker(): + """ + Drop any s3fs session a worker inherited (possible only under a fork + start method): its event loop and connections belong to the parent. The + worker builds its own on first use. + """ + pc.io_utils._S3FS_CACHE.clear() + +def _ordered_reads(pool, jobs, window): + """ + Yield _read_item(job) for each job, IN ORDER, with at most `window` reads + in flight -- so results are consumed in list order (the summation order + of the serial loop) and at most `window` tiles wait in memory. + """ + jobs = iter(jobs) + pending = collections.deque(pool.submit(_read_item, job) + for job in itertools.islice(jobs, window)) + while pending: + result = pending.popleft().result() + for job in itertools.islice(jobs, 1): + pending.append(pool.submit(_read_item, job)) + yield result + +class _TileReader: + """ + Reads the string items of a mosaic's input list, serially or in a process + pool, always yielding in list order. Non-string items (grids already in + memory) pass through untouched. + + Threads would not help: h5py holds its global lock for the whole of a + read from a Python file object, so remote reads in threads run one at a + time (measured: 8 threads = serial). Processes do not share that lock. + """ + def __init__(self, workers=1, block_size=None): + self.workers = max(1, int(workers or 1)) + self.block_size = block_size + self.pool = None + + def __enter__(self): + if self.workers > 1: + import concurrent.futures + import multiprocessing + method = _START_METHOD + if method not in multiprocessing.get_all_start_methods(): + method = 'spawn' + context = multiprocessing.get_context(method) + if method == 'forkserver': + # effective only before the server starts; it then lasts for + # the life of this process, so later pools start at once + context.set_forkserver_preload(['pointCollection']) + self.pool = concurrent.futures.ProcessPoolExecutor( + self.workers, mp_context=context, initializer=_init_worker) + return self + + def __exit__(self, *exc_info): + if self.pool is not None: + self.pool.shutdown() + self.pool = None + + def read(self, items, meta_only=False, **kwargs): + """ + yield (item, grid, exception) for each item of `items`, in order; + grid is the item itself (exception None) for a non-string item + """ + items = list(items) + jobs = [(item, meta_only, _read_kwargs(item, self.block_size, **kwargs)) + for item in items if isinstance(item, str)] + if self.pool is None: + results = map(_read_item, jobs) + else: + results = _ordered_reads(self.pool, jobs, 2*self.workers) + for item in items: + if isinstance(item, str): + grid, exc = next(results) + yield item, grid, exc + else: + yield item, item, None + class mosaic(data): def __init__(self, spacing=None, **kwargs): #self.x=None @@ -124,13 +253,24 @@ def setup_bounds_from_list(self, in_list, group=None, fields=None, bounds=None, - bands=None): - - for item in in_list.copy(): + bands=None, + reader=None): + """ + Set the mosaic's extent and spacing from its inputs, removing from + in_list any that cannot be read or fall outside bounds. reader, a + _TileReader, reads the files' metadata (in a pool, if it has one); + without one they are read here, one at a time. + """ + if reader is None: + reader = _TileReader() + metadata = reader.read(in_list.copy(), meta_only=True, group=group, bands=bands) + for item, meta, read_error in metadata: if isinstance(item, str): # read tile grid from file try: - temp=pc.grid.data().from_file(item, group=group, meta_only=True, bands=bands) + if read_error is not None: + raise read_error + temp=meta if bounds is not None: temp=temp.cropped(*bounds) if temp is not None and (len(temp.x)>0) and (len(temp.y) > 0): @@ -161,7 +301,7 @@ def setup_bounds_from_list(self, in_list, self.__update_size_and_shape__() - def setup_fields(self, item, group=None, fields=None, bands=None): + def setup_fields(self, item, group=None, fields=None, bands=None, reader=None): ''' Set up fields based on an input data item ''' @@ -170,7 +310,11 @@ def setup_fields(self, item, group=None, fields=None, bands=None): # read data grid from the first tile HDF5, use it to set the field dimensions if isinstance(item, str): - prototype=pc.grid.mosaic().from_file(item, group=group, fields=fields, bands=bands) + if reader is None: + reader = _TileReader() + _, prototype, read_error = next(reader.read([item], group=group, fields=fields, bands=bands)) + if read_error is not None: + raise read_error else: if bands is not None: item=item[:,:,bands] @@ -372,7 +516,8 @@ def add_to_band(self, item, fields, group=None, if in_band is None: temp=item else: - temp=item[:,:,in_band] + # slicing returns a plain grid.data (grid.data.__copy__) + temp=pc.grid.mosaic().from_grid(item[:,:,in_band]) else: if in_band is None: temp=pc.grid.mosaic().from_grid(item) @@ -548,12 +693,15 @@ def from_list(self, in_list, verbose=False, spacing=[None, None], bands=None, + block_size=None, + workers=1, ): """ Generate a mosaic from a list of inputs. Inputs can be strings (indicating files) or pointCollection.grid or - pointCollection.mosaic objects. + pointCollection.mosaic objects. A string may be a URI + (e.g. s3://bucket/key.h5). Parameters ---------- @@ -576,6 +724,21 @@ def from_list(self, in_list, bands : iterable, optional Bands (e.g. time slices) to read from each input, in order. If not specified, all bands in each input are read. The default is None. + block_size : int, optional + Bytes per range request for a remote HDF5 or netCDF input. None + leaves the filesystem's default, which for s3fs is 50 MiB: a + mosaic reads a few fields from each tile, so with the default a + remote tile is read almost whole. io_utils.DEFAULT_REMOTE_BLOCK_SIZE + suits this read. Ignored for local files. The default is None. + workers : int, optional + Read the files in this many processes (a pool). Tiles are still + added in list order, so the result is the same as a serial read; + at most 2 x workers read tiles wait in memory. Remote reads are + latency-bound, so this is where the speed-up is. Each worker + costs its own interpreter and imports, ~0.35 GiB resident + (measured, 2026-09): budget workers x 0.35 GiB on top of the + mosaic. The default is 1: files are read one at a time, in this + process. Returns ------- @@ -585,43 +748,73 @@ def from_list(self, in_list, """ weight = (pad is not None and pad > 0) or (feather is not None and feather>0) - self.setup_bounds_from_list(in_list, group=group, fields=fields, bounds=bounds, bands=bands) - message = self.setup_fields(in_list[0], group=group, fields=fields, bands=bands) - if message is not None: - return message - # check if using a weighted summation scheme for calculating mosaic - if weight: - if by_band: - if len(self.shape)>2: - band_list=range(self.shape[2]) + with _TileReader(workers=workers, block_size=block_size) as reader: + # with neither option the loops below hand the file names to add, + # add_to_band and replace, which read them: the original code path + prefetch = reader.pool is not None or block_size is not None + + def items(**read_kwargs): + """(item, what to add, bands for the add, read error) in order""" + if not prefetch: + for item in in_list: + yield item, item, read_kwargs.get('bands'), None + return + for item, grid, read_error in reader.read(in_list, **read_kwargs): + # a file read here is already band-selected; an in-memory + # grid is band-selected by the method it goes to + yield item, grid, (None if isinstance(item, str) else read_kwargs.get('bands')), read_error + + self.setup_bounds_from_list(in_list, group=group, fields=fields, bounds=bounds, + bands=bands, reader=reader) + message = self.setup_fields(in_list[0], group=group, fields=fields, bands=bands, + reader=reader) + if message is not None: + return message + # add, add_to_band and replace read self.fields when fields is None + read_fields = fields if fields is not None else self.fields.copy() + # check if using a weighted summation scheme for calculating mosaic + if weight: + if by_band: + if len(self.shape)>2: + band_list=range(self.shape[2]) + else: + band_list=[None] + for band in band_list: + self.invalid = np.ones(self.dimensions[0:2],dtype=bool) + self.weight = np.zeros((self.dimensions[0],self.dimensions[1])) + # if specific input bands were requested, map the output + # band index back to the corresponding input band + in_band = bands[band] if (bands is not None and band is not None) else band + in_bands = None if in_band is None else [in_band] + for item, grid, grid_bands, read_error in items(group=group, fields=read_fields, bands=in_bands): + if read_error is not None: + raise read_error + if prefetch and isinstance(item, str): + self.add_to_band(grid, group=group, fields=fields, pad=pad, feather=feather, in_band=None, out_band=band) + else: + self.add_to_band(item, group=group, fields=fields, pad=pad, feather=feather, in_band=in_band, out_band=band) + self.normalize(band=band) else: - band_list=[None] - for band in band_list: self.invalid = np.ones(self.dimensions[0:2],dtype=bool) self.weight = np.zeros((self.dimensions[0],self.dimensions[1])) - # if specific input bands were requested, map the output - # band index back to the corresponding input band - in_band = bands[band] if (bands is not None and band is not None) else band - for item in in_list: - self.add_to_band(item, group=group, fields=fields, pad=pad, feather=feather, in_band=in_band, out_band=band) - self.normalize(band=band) + # for each file in the list + for item, grid, grid_bands, read_error in items(group=group, fields=read_fields, bands=bands): + try: + if read_error is not None: + raise read_error + self.add(grid, group=group, fields=fields, pad=pad, feather=feather, bands=grid_bands) + except Exception as e: + print(f"mosaic.from_list : problem with {item} for group={group} and fields={fields}") + print(e) + self.normalize() else: - self.invalid = np.ones(self.dimensions[0:2],dtype=bool) - self.weight = np.zeros((self.dimensions[0],self.dimensions[1])) + # overwrite the mosaic with each subsequent tile # for each file in the list - for item in in_list: - try: - self.add(item, group=group, fields=fields, pad=pad, feather=feather, bands=bands) - except Exception as e: - print(f"mosaic.from_list : problem with {item} for group={group} and fields={fields}") - print(e) - self.normalize() - else: - # overwrite the mosaic with each subsequent tile - # for each file in the list - self.invalid = np.ones(self.dimensions[0:2],dtype=bool) - for item in in_list: - self.replace(item, group=group, fields=fields, bands=bands) - self.normalize(by_weight=False) + self.invalid = np.ones(self.dimensions[0:2],dtype=bool) + for item, grid, grid_bands, read_error in items(group=group, fields=read_fields, bands=bands): + if read_error is not None: + raise read_error + self.replace(grid, group=group, fields=fields, bands=grid_bands) + self.normalize(by_weight=False) return self diff --git a/pointCollection/io_utils.py b/pointCollection/io_utils.py index 3b13267..274014a 100644 --- a/pointCollection/io_utils.py +++ b/pointCollection/io_utils.py @@ -367,6 +367,32 @@ def as_gdal_path(filename): # unknown scheme: hand it to GDAL as-is and let GDAL report the problem return filename +def glob_remote(pattern, fs=None): + """ + The remote counterpart of glob.glob: list the objects matching a URI + pattern such as 's3://bucket/tiles/matched/E*.h5'. + + Parameters + ---------- + pattern : str + URI with glob wildcards in its key. + fs : fsspec filesystem, optional + filesystem to list with. If None, a cached session on the default AWS + credential chain is used (get_s3fs(daac=None)): a remote glob is for + buckets we own, as a gridded read of an s3:// file is. + + Returns + ------- + list of str + matching URIs, with the pattern's scheme, SORTED -- unlike glob.glob, + whose order is whatever the directory gives -- so a mosaic built from + the list sums its tiles in the same order every time. + """ + scheme = pattern.partition('://')[0] + if fs is None: + fs = get_s3fs(daac=None) + return sorted(f'{scheme}://' + fs._strip_protocol(path) for path in fs.glob(pattern)) + def path_exists(filename, fs=None, assume_remote_exists=True): """ Check whether a local or remote file exists. diff --git a/pointCollection/scripts/make_mosaic.py b/pointCollection/scripts/make_mosaic.py index beba1f6..0f44658 100755 --- a/pointCollection/scripts/make_mosaic.py +++ b/pointCollection/scripts/make_mosaic.py @@ -8,7 +8,8 @@ COMMAND LINE OPTIONS: --help: list the command line options - -d X, --directory X: directory to run + -d X, --directory X: directory to run; may be a URI (e.g. s3://bucket/dir), + whose tiles are then read in place (-O must be a local absolute path) -g X, --glob_string X: quoted string to pass to glob to find the files --proj4: projection string for the tiles and output mosaic -r X, --range X: valid range of tiles to read [xmin,xmax,ymin,ymax] @@ -26,8 +27,11 @@ -s, --show: create plot of output mosaic -m X, --mode X: Local permissions mode of the output mosaic -N --ignore_Nodata: ignore nodata values (except Nan) in inputs + --block_size X: bytes per range request for remote tiles + -j X, --workers X: read the tiles in X processes UPDATE HISTORY: + Updated 09/2026: remote (URI) directories; --block_size and --workers Updated 05/2024: allow cropping in time for 3D fields Updated 01/2021: added option for setting projection attributes Updated 10/2021: added option for using a non-weighted summation @@ -120,6 +124,11 @@ def main(): parser.add_argument('--mode','-m', type=lambda x: int(x,base=8), default=0o775, help='permissions mode of output mosaic') + parser.add_argument('--block_size', type=int, default=None, + help='bytes per range request when the tiles are remote (default: ' + 'pointCollection.io_utils.DEFAULT_REMOTE_BLOCK_SIZE)') + parser.add_argument('--workers','-j', type=int, default=1, + help='read the tiles in this many processes (default 1)') try: assert(len(sys.argv)>1) args=parser.parse_args() @@ -140,13 +149,27 @@ def main(): if args.verbose: print("searching in directory "+args.directory+" with glob string:"+"["+str(args.glob_string)+"]") # find list of valid files + # A URI directory (e.g. s3://bucket/region) is listed with a remote glob, + # and the tiles are read in place. The output is written with h5py, which + # needs a local file, so it must then be given as a local absolute path: + # a relative one would be joined onto the URI below. + remote = pc.io_utils.is_remote_path(args.directory) + if remote: + if pc.io_utils.is_remote_path(args.output) or not os.path.isabs(args.output): + parser.error(f'--directory {args.directory} is remote, so --output must be a ' + f'local absolute path, not {args.output}') + if args.block_size is None: + args.block_size = pc.io_utils.DEFAULT_REMOTE_BLOCK_SIZE + find_files = pc.io_utils.glob_remote + else: + find_files = glob.glob if isinstance(args.glob_string, str): - initial_file_list = glob.glob(args.directory +'/'+args.glob_string) + initial_file_list = find_files(args.directory +'/'+args.glob_string) else: initial_file_list = [] for glob_string in args.glob_string: - initial_file_list += glob.glob(args.directory +'/'+glob_string) + initial_file_list += find_files(args.directory +'/'+glob_string) if args.verbose: print(f"initial file list contains {len(initial_file_list)} files") @@ -176,7 +199,9 @@ def main(): group=args.in_group, pad=args.pad, feather=args.feather, - by_band=args.by_band) + by_band=args.by_band, + block_size=args.block_size, + workers=args.workers) if isinstance(mosaic, str): if args.verbose: print(f"pc.grid.mosaic failed for group {args.in_group} and fields {args.fields} with message:") diff --git a/tests/test_mosaic_list.py b/tests/test_mosaic_list.py new file mode 100644 index 0000000..ca3eb58 --- /dev/null +++ b/tests/test_mosaic_list.py @@ -0,0 +1,179 @@ +""" +grid.mosaic.from_list and make_mosaic.py: reading the tiles with a block_size, +in a process pool (workers), and from a remote directory must give exactly the +mosaic the original serial, local read gives. + +No network: remote files live in fsspec's in-memory filesystem. Only a +FORKED worker sees it (and the patched session getter), so the remote tests +with workers switch the pool to fork; the local ones run the default, +forkserver. +""" +import importlib +import os +import sys +import fsspec +import numpy as np +import pytest +import pointCollection as pc + +CENTERS = [0., 40., 80.] +T = np.arange(4.) + + +def write_tile(filename, xc, yc, rng): + """a 60 x 60 tile at 2-unit spacing, a 2-D and a 3-D group, some NaNs""" + x = xc + np.arange(-30., 30.1, 2.) + y = yc + np.arange(-30., 30.1, 2.) + z0 = rng.normal(size=(y.size, x.size)) + dz = rng.normal(size=(y.size, x.size, T.size)) + dz[rng.random(dz.shape) < 0.05] = np.nan + pc.grid.data().from_dict({'x': x, 'y': y, 'z0': z0, + 'cell_area': np.ones_like(z0)}).to_h5(filename, group='z0', replace=True) + pc.grid.data().from_dict({'x': x, 'y': y, 't': T, 'dz': dz}).to_h5(filename, group='dz', replace=False) + + +@pytest.fixture(scope='module') +def tiles(tmp_path_factory): + d = tmp_path_factory.mktemp('tiles') + rng = np.random.default_rng(3) + files = [] + for xc in CENTERS: + for yc in CENTERS: + files.append(str(d / f'E{int(xc)}_N{int(yc)}.h5')) + write_tile(files[-1], xc, yc, rng) + return sorted(files) + + +@pytest.fixture +def memory_tiles(tiles, monkeypatch): + """the same tiles in a memory filesystem, returned as memory:// URIs""" + fs = fsspec.filesystem('memory') + uris = [] + for f in tiles: + uri = 'memory://bucket/region/matched/' + os.path.basename(f) + with open(f, 'rb') as fh: + fs.pipe(uri, fh.read()) + uris.append(uri) + monkeypatch.setattr(pc.io_utils, 'get_s3fs', lambda daac=None, **kw: fs) + # pc.grid.mosaic is the class; the module is shadowed by it + monkeypatch.setattr(importlib.import_module('pointCollection.grid.mosaic'), '_START_METHOD', 'fork') + yield uris + fs.rm('memory://bucket', recursive=True) + + +CASES = { + 'weighted by band': dict(group='dz', fields=['dz'], pad=4, feather=8, by_band=True), + 'weighted all bands': dict(group='dz', fields=['dz'], pad=4, feather=8, by_band=False), + 'weighted 2-D': dict(group='z0', fields=['z0', 'cell_area'], pad=4, feather=8, by_band=False), + 'replace': dict(group='z0', fields=['z0'], pad=None, feather=None), + 'selected bands': dict(group='dz', fields=['dz'], pad=4, feather=8, by_band=True, bands=[1, 3]), + 'all fields': dict(group='z0', fields=None, pad=4, feather=8, by_band=False), +} + + +def mosaic(files, **kwargs): + return pc.grid.mosaic().from_list(list(files), **kwargs) + + +def assert_same(a, b): + assert a.fields == b.fields + assert np.array_equal(a.x, b.x) and np.array_equal(a.y, b.y) + for field in a.fields: + assert np.array_equal(getattr(a, field), getattr(b, field), equal_nan=True), field + + +@pytest.mark.parametrize('case', CASES) +@pytest.mark.parametrize('options', [dict(workers=3), dict(block_size=4096), dict(workers=2, block_size=4096)], + ids=['workers', 'block_size', 'both']) +def test_local_read_options_change_nothing(tiles, case, options): + assert_same(mosaic(tiles, **CASES[case]), mosaic(tiles, **CASES[case], **options)) + + +# fork warns once the test process has threads (fsspec's loop); that is the +# reason the default is forkserver, and harmless for an in-memory filesystem +FORK_WARNING = pytest.mark.filterwarnings('ignore:This process .* is multi-threaded:DeprecationWarning') + + +@FORK_WARNING +@pytest.mark.parametrize('case', CASES) +def test_remote_tiles_give_the_local_mosaic(tiles, memory_tiles, case): + reference = mosaic(tiles, **CASES[case]) + assert_same(reference, mosaic(memory_tiles, **CASES[case], block_size=4096)) + assert_same(reference, mosaic(memory_tiles, **CASES[case], block_size=4096, workers=3)) + + +def test_default_start_method_is_forkserver(): + # plain fork deadlocks or fails ("not fork-safe") after s3fs has started + assert importlib.import_module('pointCollection.grid.mosaic')._START_METHOD == 'forkserver' + + +def test_block_size_reaches_the_remote_open(memory_tiles, monkeypatch): + seen = [] + open_remote = pc.io_utils.open_remote + def spy(filename, *args, block_size=None, **kwargs): + seen.append(block_size) + return open_remote(filename, *args, block_size=block_size, **kwargs) + monkeypatch.setattr(pc.io_utils, 'open_remote', spy) + mosaic(memory_tiles, **CASES['weighted 2-D'], block_size=12345) + assert seen and set(seen) == {12345} + + +def test_unreadable_tile_is_dropped_the_same_way(tiles, tmp_path): + bad = str(tmp_path / 'E999_N999.h5') + with open(bad, 'w') as fh: + fh.write('not an hdf5 file') + files = tiles[:4] + [bad] + tiles[4:] + reference = mosaic(files, **CASES['weighted by band']) + assert_same(reference, mosaic(files, **CASES['weighted by band'], workers=3)) + assert_same(reference, mosaic(tiles, **CASES['weighted by band'])) + + +@pytest.mark.parametrize('kwargs', [dict(by_band=True), dict(by_band=True, bands=[0, 2]), + dict(by_band=False), dict(by_band=False, bands=[0, 2])], + ids=['by band', 'by band, selected bands', 'all bands', 'selected bands']) +def test_grids_in_memory_pass_through(tiles, kwargs): + grids = [pc.grid.mosaic().from_file(f, group='dz', fields=['dz']) for f in tiles[:3]] + mixed = grids + tiles[3:] + kwargs = dict(group='dz', fields=['dz'], pad=4, feather=8, **kwargs) + assert_same(mosaic(mixed, **kwargs), mosaic(mixed, **kwargs, workers=2, block_size=4096)) + + +@pytest.mark.parametrize('bands', [None, [1, 3]]) +def test_in_memory_mosaic_by_band_matches_its_file(tiles, bands): + # add_to_band used to slice an in-memory mosaic to a plain grid.data and + # fail on update_spacing + grids = [pc.grid.mosaic().from_file(f, group='dz', fields=['dz']) for f in tiles] + kwargs = dict(group='dz', fields=['dz'], pad=4, feather=8, by_band=True, bands=bands) + assert_same(mosaic(tiles, **kwargs), mosaic(grids, **kwargs)) + + +def test_glob_remote_is_sorted_and_keeps_the_scheme(memory_tiles): + found = pc.io_utils.glob_remote('memory://bucket/region/matched/E*_N0.h5') + assert found == sorted(found) + assert [os.path.basename(f) for f in found] == ['E0_N0.h5', 'E40_N0.h5', 'E80_N0.h5'] + assert all(f.startswith('memory://') for f in found) + + +def run_make_mosaic(monkeypatch, argv): + from pointCollection.scripts import make_mosaic + monkeypatch.setattr(sys, 'argv', ['make_mosaic.py'] + argv) + make_mosaic.main() + + +@FORK_WARNING +def test_make_mosaic_remote_directory(tiles, memory_tiles, tmp_path, monkeypatch): + common = ['-g', 'E*.h5', '-w', '-p', '4', '-f', '8', '--in_group', 'dz/', '-F', 'dz', '-R'] + run_make_mosaic(monkeypatch, ['-d', os.path.dirname(tiles[0]), '-O', str(tmp_path / 'local.h5')] + common) + run_make_mosaic(monkeypatch, ['-d', 'memory://bucket/region/matched', '-O', str(tmp_path / 'remote.h5'), + '-j', '2'] + common) + local = pc.grid.data().from_h5(str(tmp_path / 'local.h5'), group='dz') + remote = pc.grid.data().from_h5(str(tmp_path / 'remote.h5'), group='dz') + assert np.array_equal(local.dz, remote.dz, equal_nan=True) + + +@pytest.mark.parametrize('output', ['mosaic.h5', 's3://bucket/mosaic.h5']) +def test_make_mosaic_remote_directory_needs_a_local_absolute_output(memory_tiles, monkeypatch, output): + with pytest.raises(SystemExit) as exit_info: + run_make_mosaic(monkeypatch, ['-d', 'memory://bucket/region/matched', '-g', 'E*.h5', + '--in_group', 'z0/', '-F', 'z0', '-O', output]) + assert exit_info.value.code == 2 From 2b69f1ab8ece2803d6da055c0f2053d065fa583e Mon Sep 17 00:00:00 2001 From: Ben Smith Date: Wed, 30 Sep 2026 22:41:31 +0000 Subject: [PATCH 2/2] make_mosaic: sort a local glob, as glob_remote does CI failed test_make_mosaic_remote_directory: the local run took its tiles in the directory's order (glob.glob), the remote run in sorted order, and a weighted mosaic depends on summation order in the last bit (reversing 9 test tiles: max |diff| 2.2e-16, same NaNs). The test passed on the ADE only because that directory happened to list in sorted order. Sorting makes a local mosaic the same on every filesystem and the same as the remote one. It can differ in the last bit from a mosaic made before, from an unsorted directory listing. test_make_mosaic_sorts_a_local_glob reproduces CI's condition anywhere (a glob that returns the files reversed); it fails without the sort. Co-Authored-By: Claude Opus 5.5 --- pointCollection/scripts/make_mosaic.py | 6 +++++- tests/test_mosaic_list.py | 17 +++++++++++++++++ 2 files changed, 22 insertions(+), 1 deletion(-) diff --git a/pointCollection/scripts/make_mosaic.py b/pointCollection/scripts/make_mosaic.py index 0f44658..12d3b90 100755 --- a/pointCollection/scripts/make_mosaic.py +++ b/pointCollection/scripts/make_mosaic.py @@ -162,7 +162,11 @@ def main(): args.block_size = pc.io_utils.DEFAULT_REMOTE_BLOCK_SIZE find_files = pc.io_utils.glob_remote else: - find_files = glob.glob + # SORTED, like glob_remote: a weighted mosaic sums its tiles in list + # order, and a different order changes the last bit, so an unsorted + # glob (whatever order the directory gives) made the same mosaic + # differ between filesystems, and between local and remote tiles + find_files = lambda pattern: sorted(glob.glob(pattern)) if isinstance(args.glob_string, str): initial_file_list = find_files(args.directory +'/'+args.glob_string) diff --git a/tests/test_mosaic_list.py b/tests/test_mosaic_list.py index ca3eb58..15c18a1 100644 --- a/tests/test_mosaic_list.py +++ b/tests/test_mosaic_list.py @@ -171,6 +171,23 @@ def test_make_mosaic_remote_directory(tiles, memory_tiles, tmp_path, monkeypatch assert np.array_equal(local.dz, remote.dz, equal_nan=True) +def test_make_mosaic_sorts_a_local_glob(tiles, tmp_path, monkeypatch): + # a weighted mosaic depends on summation order in the last bit; a + # directory that lists its files in another order (as CI's did) must not + # change the result + import glob + from pointCollection.scripts import make_mosaic + common = ['-d', os.path.dirname(tiles[0]), '-g', 'E*.h5', '-w', '-p', '4', '-f', '8', + '--in_group', 'dz/', '-F', 'dz', '-R'] + run_make_mosaic(monkeypatch, common + ['-O', str(tmp_path / 'listed.h5')]) + listed_glob = glob.glob # make_mosaic.glob IS the glob module: keep the real one + monkeypatch.setattr(make_mosaic.glob, 'glob', lambda pattern: listed_glob(pattern)[::-1]) + run_make_mosaic(monkeypatch, common + ['-O', str(tmp_path / 'reversed.h5')]) + listed = pc.grid.data().from_h5(str(tmp_path / 'listed.h5'), group='dz') + reversed_ = pc.grid.data().from_h5(str(tmp_path / 'reversed.h5'), group='dz') + assert np.array_equal(listed.dz, reversed_.dz, equal_nan=True) + + @pytest.mark.parametrize('output', ['mosaic.h5', 's3://bucket/mosaic.h5']) def test_make_mosaic_remote_directory_needs_a_local_absolute_output(memory_tiles, monkeypatch, output): with pytest.raises(SystemExit) as exit_info: