diff --git a/.gitignore b/.gitignore
index 2cc96b5..8dfa32c 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,5 +1,6 @@
 __pycache__/
 .pytest_cache/
 examples/large_geotiff/
+examples/multi_netcdf/
 mapshader.egg-info/
 mapshader/.version
diff --git a/example_multiband/create_multiband_netcdf.py b/example_multiband/create_multiband_netcdf.py
deleted file mode 100644
index dabf519..0000000
--- a/example_multiband/create_multiband_netcdf.py
+++ /dev/null
@@ -1,71 +0,0 @@
-# Create multiple netcdf files that are adjacent tiles containing the same
-# multiple variables.
-
-import numpy as np
-import os
-import rioxarray  # noqa: F401
-import xarray as xr
-
-
-lat_limits = (0, 3)  # Integer min and max latitude
-lon_limits = (0, 4)  # Integer min and max longitude
-ny, nx = 4, 5        # Grid points per tile.
-
-scale = 10
-offset = 1
-
-
-dx = scale / nx
-dy = scale / ny
-
-rng = np.random.default_rng(92741)
-
-red_min = np.nan
-red_max = np.nan
-green_min = np.nan
-green_max = np.nan
-blue_min = np.nan
-blue_max = np.nan
-
-for lat in range(lat_limits[0], lat_limits[1]):
-    for lon in range(lon_limits[0], lon_limits[1]):
-        # Randomly remove some tiles.
-        if rng.random(1)[0] > 0.8:
-            continue
-
-        x = np.linspace(lon*scale + offset + dx/2, (lon + 1)*scale + offset - dx/2, nx)
-        y = np.linspace(lat*scale + offset + dy/2, (lat + 1)*scale + offset - dy/2, ny)
-        y = y[::-1]
-        filename = os.path.join("example_multiband", f"dummy_{lon}_{lat}.nc")
-
-        red = rng.random(size=(ny, nx), dtype=np.float32)
-        green = rng.random(size=(ny, nx), dtype=np.float32) + lon
-        blue = rng.random(size=(ny, nx), dtype=np.float32) + lat
-
-        dims = ["y", "x"]
-        ds = xr.Dataset(
-            data_vars=dict(
-                red=(dims, red),
-                green=(dims, green),
-                blue=(dims, blue),
-            ),
-            coords=dict(y=y, x=x),
-            attrs=dict(
-                description="dummy data",
-                crs="epsg:4326",
-            ),
-        )
-
-        ds.to_netcdf(filename)
-        print(f"Written file {filename}", ds.rio.bounds(), x.min(), y.min(), x.max(), y.max())
-
-        red_min = np.nanmin((red_min, red.min()))
-        red_max = np.nanmax((red_max, red.max()))
-        green_min = np.nanmin((green_min, green.min()))
-        green_max = np.nanmax((green_max, green.max()))
-        blue_min = np.nanmin((blue_min, blue.min()))
-        blue_max = np.nanmax((blue_max, blue.max()))
-
-print(f"red limits: {red_min} {red_max}")
-print(f"green limits: {green_min} {green_max}")
-print(f"blue limits: {blue_min} {blue_max}")
diff --git a/example_multiband/multiband_netcdf.yaml b/example_multiband/multiband_netcdf.yaml
deleted file mode 100644
index d7bd520..0000000
--- a/example_multiband/multiband_netcdf.yaml
+++ /dev/null
@@ -1,98 +0,0 @@
----
-
-metadata:
-  version: 1
-
-sources:
-  - name: Dummy Multiband Netcdf Red
-    key: multiband-netcdf-red
-    text: Dummy multiband netcdf red
-    description: Dummy multiband netcdf red
-    geometry_type: raster
-    shade_how: linear
-    cmap:
-      - white
-      - red
-    span:
-      - 0
-      - 1
-    raster_padding: 0
-    raster_interpolate: linear
-    xfield: geometry
-    yfield: geometry
-    filepath: example_multiband/dummy_*.nc
-    band: red
-    transforms:
-      - name: reproject_raster
-        args:
-          epsg: 3857
-      - name: build_raster_overviews
-        args:
-          levels:
-            0: 256
-            1: 512
-            2: 1024
-    service_types:
-      - tile
-
-  - name: Dummy Multiband Netcdf Green
-    key: multiband-netcdf-green
-    text: Dummy multiband netcdf green
-    description: Dummy multiband netcdf green
-    geometry_type: raster
-    shade_how: linear
-    cmap:
-      - yellow
-      - green
-    span:
-      - 0
-      - 4
-    raster_padding: 0
-    raster_interpolate: linear
-    xfield: geometry
-    yfield: geometry
-    filepath: example_multiband/dummy_*.nc
-    band: green
-    transforms:
-      - name: reproject_raster
-        args:
-          epsg: 3857
-      - name: build_raster_overviews
-        args:
-          levels:
-            0: 256
-            1: 512
-            2: 1024
-    service_types:
-      - tile
-
-  - name: Dummy Multiband Netcdf Blue
-    key: multiband-netcdf-blue
-    text: Dummy multiband netcdf blue
-    description: Dummy multiband netcdf blue
-    geometry_type: raster
-    shade_how: linear
-    cmap:
-      - yellow
-      - blue
-    span:
-      - 0
-      - 3
-    raster_padding: 0
-    raster_interpolate: linear
-    xfield: geometry
-    yfield: geometry
-    filepath: example_multiband/dummy_*.nc
-    band: blue
-    transforms:
-      - name: reproject_raster
-        args:
-          epsg: 3857
-      - name: build_raster_overviews
-        args:
-          levels:
-            0: 256
-            1: 512
-            2: 1024
-    service_types:
-      - tile
diff --git a/example_multiband/run_multiband_netcdf.py b/example_multiband/run_multiband_netcdf.py
deleted file mode 100644
index e1536f8..0000000
--- a/example_multiband/run_multiband_netcdf.py
+++ /dev/null
@@ -1,46 +0,0 @@
-from mapshader.sources import MapSource
-from mapshader.core import render_map
-import PIL
-import yaml
-
-
-yaml_file = "example_multiband/multiband_netcdf.yaml"
-
-with open(yaml_file, 'r') as f:
-    content = f.read()
-    config_obj = yaml.safe_load(content)
-    source_objs = config_obj['sources']
-
-sources = []
-for source_obj in source_objs:
-    sources.append(MapSource.from_obj(source_obj))
-
-print(sources)
-
-
-# Want tile containing slightly positive lat and lon
-def run(x, y, z, index, j):
-    print("XYZ", x, y, z)
-
-    img = render_map(source, x=x, y=y, z=z, height=256, width=256)  # xarray Image 8x8
-
-    bytes = img.to_bytesio()
-    pillow_img = PIL.Image.open(bytes)
-    pillow_img.save(f"out_{j}_{index}.png", format="png")
-
-
-def xyz_contains_data(z):
-    ntiles = 2**z  # In both x and y directions.
-    x = ntiles // 2
-    y = max(x-1, 0)
-    return x, y, z
-
-
-for j in range(3):
-    source = sources[j].load()
-
-    i = 0
-    for z in range(0, 4):
-        x, y, z = xyz_contains_data(z)
-        run(x, y, z, i, j)
-        i += 1
diff --git a/examples/multi_netcdf.yaml b/examples/multi_netcdf.yaml
new file mode 100644
index 0000000..f45eec2
--- /dev/null
+++ b/examples/multi_netcdf.yaml
@@ -0,0 +1,69 @@
+---
+sources:
+  - name: Multi Netcdf constant
+    key: multi-netcdf-constant
+    text: multi netcdf constant
+    description: multi netcdf constant
+    geometry_type: raster
+    shade_how: linear
+    cmap:
+      - white
+      - red
+    span:
+      - 0
+      - 1
+    raster_interpolate: linear
+    xfield: geometry
+    yfield: geometry
+    filepath: examples/multi_netcdf/multi_*.nc
+    band: constant
+    transforms:
+      - name: reproject_raster
+        args:
+          epsg: 3857
+      - name: build_raster_overviews
+        args:
+          levels:
+            0: 256
+            1: 512
+            2: 1024
+            3: 2048
+            4: 4096
+            5: 8192
+            6: 16384
+    service_types:
+      - tile
+
+  - name: Multi Netcdf increasing
+    key: multi-netcdf-increasing
+    text: multi netcdf increasing
+    description: multi netcdf increasing
+    geometry_type: raster
+    shade_how: linear
+    cmap:
+      - black
+      - yellow
+    span:
+      - 0
+      - 10
+    raster_interpolate: linear
+    xfield: geometry
+    yfield: geometry
+    filepath: examples/multi_netcdf/multi_*.nc
+    band: increasing
+    transforms:
+      - name: reproject_raster
+        args:
+          epsg: 3857
+      - name: build_raster_overviews
+        args:
+          levels:
+            0: 256
+            1: 512
+            2: 1024
+            3: 2048
+            4: 4096
+            5: 8192
+            6: 16384
+    service_types:
+      - tile
diff --git a/mapshader/multifile.py b/mapshader/multifile.py
index 67069a7..d383b2d 100644
--- a/mapshader/multifile.py
+++ b/mapshader/multifile.py
@@ -280,14 +280,13 @@ def load_bounds(self, xmin, ymin, xmax, ymax, band, transforms):
                     da.rio.set_crs(crs, inplace=True)
                     arrays.append(da)
 
-        if len(arrays) == 1:
-            merged = arrays[0]
-        else:
-            merged = merge_arrays(arrays)
+            if len(arrays) == 1:
+                merged = arrays[0]
+            else:
+                merged = merge_arrays(arrays)
 
-        merged = merged.squeeze()
+            merged = merged.squeeze()
 
-        with self._lock:
             merged = self._apply_transforms(merged, transforms)
 
         return merged
diff --git a/mapshader/tests/test_multi_netcdf.py b/mapshader/tests/test_multi_netcdf.py
new file mode 100644
index 0000000..2b906cd
--- /dev/null
+++ b/mapshader/tests/test_multi_netcdf.py
@@ -0,0 +1,93 @@
+from glob import glob
+import numpy as np
+import os
+import pytest
+import rioxarray  # noqa: F401
+import subprocess
+import xarray as xr
+
+
+def check_and_create_multi_netcdf():
+    directory = os.path.join("examples", "multi_netcdf")
+    if not os.path.isdir(directory):
+        print(f"Creating directory {directory}")
+        os.makedirs(directory)
+
+    # Number of files.
+    nx_tile = 50
+    ny_tile = 20
+
+    filenames = glob(os.path.join(directory, "multi_*.nc"))
+    if len(filenames) == nx_tile*ny_tile:
+        # If the target directory contains the correct number of netcdf files,
+        # assume it is correct and do not regenerate it.
+        return
+
+    # Grid points per tile.
+    ny, nx = 50, 100
+    lon_start = -90.0
+    lat_start = -40.0
+    scale = 1.5  # Degrees per tile.
+    dx = scale / nx
+    dy = scale / ny
+
+    rng = np.random.default_rng(92741)
+    dims = ["y", "x"]
+
+    for j in range(ny_tile):
+        for i in range(nx_tile):
+            x = np.linspace(lon_start + i*scale + dx/2, lon_start + (i+1)*scale - dx/2, nx)
+            y = np.linspace(lat_start + j*scale + dy/2, lat_start + (j+1)*scale - dy/2, ny)
+            y = y[::-1]
+            filename = f"multi_{j:03}_{i:03}.nc"
+            full_filename = os.path.join(directory, filename)
+
+            # 2 bands.
+            constant = rng.random(size=(ny, nx), dtype=np.float32)
+            increasing = rng.random(size=(ny, nx), dtype=np.float32) + j
+
+            ds = xr.Dataset(
+                data_vars=dict(
+                    constant=(dims, constant),
+                    increasing=(dims, increasing),
+                ),
+                coords=dict(y=y, x=x),
+                attrs=dict(
+                    description="Test multi netcdf data",
+                    crs="epsg:4326",
+                ),
+            )
+
+            ds.to_netcdf(full_filename)
+            print(f"Written file {full_filename}")
+
+
+@pytest.mark.large
+def test_multi_netcdf_create_overviews():
+    check_and_create_multi_netcdf()
+
+    yaml_file = os.path.join("examples", "multi_netcdf.yaml")
+    cmd = ["mapshader", "build-raster-overviews", yaml_file]
+    subprocess.run(cmd)
+
+
+@pytest.mark.large
+@pytest.mark.parametrize("z, shape", [
+    [0, (28, 57)],
+    [1, (80, 145)],
+    [2, (158, 288)],
+    [3, (312, 572)],
+    [4, (619, 1141)],
+    [5, (1235, 2278)],
+    [6, (2467, 4554)],
+])
+@pytest.mark.parametrize("band", ("constant", "increasing"))
+def test_multi_netcdf_overview_files(band, z, shape):
+    directory = os.path.join("examples", "multi_netcdf", "overviews")
+    filename = os.path.join(directory, f"{z}_{band}.nc")
+    assert os.path.isfile(filename)
+
+    ds = xr.open_dataset(filename, chunks=dict(y=512, x=512))
+    da = ds[band]
+    assert da.dtype == np.float32
+    assert da.shape == shape