Limit to number of files that can be merged with merge_arrays?

[JamesSample]

Is there a limit to the number of data arrays that can be merged by merge_arrays?

I'm trying to merge 146 DEM tiles that all have the same cell size and CRS:

import glob
import rioxarray as rio
from rioxarray.merge import merge_arrays

search_path = f"/home/jovyan/shared/01_datasets/spatial/dtm10_utm32/raw/*.tif"
flist = sorted(glob.glob(search_path))
print(f"{len(flist)} files to process.")

print("Opening files...")
srcs = [rio.open_rasterio(fpath, mask_and_scale=True) for fpath in flist]

print("Merging tiles...")
rds = merge_arrays(srcs)

This produces the following output:

146 files to process.
Opening files...
Merging tiles...

And then gives this error:

---------------------------------------------------------------------------
RasterioIOError                           Traceback (most recent call last)
/tmp/ipykernel_5294/3413469784.py in <module>
      7 
      8 print("Merging tiles...")
----> 9 rds = merge_arrays(srcs)

/opt/conda/lib/python3.9/site-packages/rioxarray/merge.py in merge_arrays(dataarrays, bounds, res, nodata, precision, method, crs, parse_coordinates)
    160 
    161     # use rasterio to merge
--> 162     merged_data, merged_transform = _rio_merge(
    163         rioduckarrays,
    164         **{key: val for key, val in input_kwargs.items() if val is not None},

/opt/conda/lib/python3.9/site-packages/rasterio/merge.py in merge(datasets, bounds, res, nodata, dtype, precision, indexes, output_count, resampling, method, target_aligned_pixels, dst_path, dst_kwds)
    202 
    203         try:
--> 204             first_colormap = first.colormap(1)
    205         except ValueError:
    206             first_colormap = None

rasterio/_base.pyx in rasterio._base.DatasetBase.colormap()

rasterio/_base.pyx in rasterio._base.DatasetBase.band()

rasterio/_base.pyx in rasterio._base.DatasetBase.handle()

RasterioIOError: Dataset is closed: /home/jovyan/shared/01_datasets/spatial/dtm10_utm32/raw/6403_1_10m_z32.tif

However, if I modify the last line of the code to only consider the first 50 datasets in the list, everything works OK:

rds = merge_arrays(srcs[:50])

6403_1_10m_z32.tif is the first raster in the list, so it shouldn't be closed (and it isn't when I just process the first 50 items).

My system has plenty of resources available (total memory consumption peaks at around 4% based on the output from top). The following rasterio code also works successfully:

    # Merge all datsets for this zone, keeping original CRS
    print("Merging tiles...")
    srcs = [rasterio.open(fpath) for fpath in flist]
    out_meta = srcs[0].meta
    out_image, out_trans = merge(srcs)
    out_meta.update(
        {
            "driver": "GTiff",
            "height": out_image.shape[1],
            "width": out_image.shape[2],
            "transform": out_trans,
            "compress": "lzw",
        }
    )

    with rasterio.open(mosaic_gtiff, "w", **out_meta) as dest:
        dest.write(out_image)
    [src.close() for src in srcs]

but I'm producing a notebook for relative beginners and the rioxarray code is much nicer/more accessible.

I'm on Linux with rioxarray=0.10.1 installed from conda-forge.

Thanks in advance for any help or advice! :)

[snowman2]

This is what currently happens:

rioxarray/rioxarray/merge.py

Line 38 in 033668f

self.colormap = rio_file.colormap

What probably needs to occur is one of these 2 options:

1. Load the contents of colormap into a _colormap attribute and add a colormap method to return the contents.
2. Add colormap method that calls self.xds.rio._manager.acquire() before calling colormap

[snowman2]

Workaround: https://xarray.pydata.org/en/stable/generated/xarray.set_options.html

file_cache_maxsize is defaulted to 128, if you bump it up to more than the number of files, it should work.

[snowman2]

PR for reference.

[JamesSample]

Installed from your colormap branch and it works like a charm - thank you!