diff --git a/examples/atlas_streamlines.ipynb b/examples/atlas_streamlines.ipynb
new file mode 100644
index 0000000..81b3dcf
--- /dev/null
+++ b/examples/atlas_streamlines.ipynb
@@ -0,0 +1,166 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "5823062b-e25f-45aa-96d7-65c83845af61",
+   "metadata": {},
+   "source": [
+    "# Working with cortical streamlines\n",
+    "\n",
+    "In this example we walk through various methods that use the Allen streamlines. Streamlines are the paths that connect the surface of the isocortex to the lower white matter surface while following the curvature of these surfaces."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "46e7e0e7-4295-40a5-a5b3-1bc0fc84f07e",
+   "metadata": {},
+   "source": [
+    "## Depth lookup\n",
+    "We can use the streamlines to find the depth of coordinates in 3D space relative to the surface of the cortex. Note as the streamlines are only defined in the cortex so this lookup will only work for coordinates that lie within the Isocortex of the Allen atlas volume. In this example we obtain the depths in percentage, where 0 % indicates the surface of the isocortex and 100 % the white matter boundary."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "0a3fddd1-9ddd-4dca-b256-3dd7d29bad67",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from iblatlas.streamlines.utils import xyz_to_depth\n",
+    "from iblatlas.atlas import AllenAtlas\n",
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt\n",
+    "\n",
+    "ba = AllenAtlas()\n",
+    "# In this example we are going to get out the coordinates of the voxels that lie within the region MOp\n",
+    "region_id = ba.regions.acronym2index('MOp')[1][0][0]\n",
+    "# MOp is in the beryl mapping so we map our label volume to this mapping\n",
+    "ba.label = ba.regions.mappings['Beryl-lr'][ba.label]\n",
+    "# These are the indices in the volume, order is AP, ML, DV\n",
+    "ixyz = np.where(ba.label == region_id)\n",
+    "# Convert these to xyz coordinates relative to bregma in order ML, AP, DV\n",
+    "xyz = ba.bc.i2xyz(np.c_[ixyz[1], ixyz[0], ixyz[2]])\n",
+    "\n",
+    "depth_per = xyz_to_depth(xyz) * 100\n",
+    "\n",
+    "# We can assign a volume\n",
+    "depth_per_vol = np.full(ba.image.shape, np.nan)\n",
+    "depth_per_vol[ixyz[0], ixyz[1], ixyz[2]] = depth_per\n",
+    "\n",
+    "fig, ax = plt.subplots()\n",
+    "ba.plot_cslice(620/1e6, volume='image', ax=ax)\n",
+    "ba.plot_cslice(620/1e6, volume='volume', region_values=depth_per_vol, ax=ax, cmap='viridis', vmin=0, vmax=100)\n",
+    "ax.set_axis_off()\n",
+    "cbar = fig.colorbar(ax.get_images()[1], ax=ax)\n",
+    "cbar.set_label('Depth from surface of cortex (%)')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "c1da61f0-7281-45d8-8bd7-33be0d373302",
+   "metadata": {},
+   "source": [
+    "We can also do this look up in um relative to the surface of the cortex."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "047509db-fb6b-46f3-8cb2-817b12223bed",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "depth_um = xyz_to_depth(xyz, per=False)\n",
+    "\n",
+    "# We can assign a volume\n",
+    "depth_um_vol = np.full(ba.image.shape, np.nan)\n",
+    "depth_um_vol[ixyz[0], ixyz[1], ixyz[2]] = depth_um\n",
+    "\n",
+    "fig, ax = plt.subplots()\n",
+    "ba.plot_cslice(620/1e6, volume='image', ax=ax)\n",
+    "ba.plot_cslice(620/1e6, volume='volume', region_values=depth_um_vol, ax=ax, cmap='tab20', vmin=0, vmax=2000)\n",
+    "ax.set_axis_off()\n",
+    "cbar = fig.colorbar(ax.get_images()[1], ax=ax)\n",
+    "cbar.set_label('Depth from surface of cortex (um)')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "bdbe5c0b-d7e3-441e-82ab-8bbcbd013bdd",
+   "metadata": {},
+   "source": [
+    "## Flatmap projection\n",
+    "The streamlines can also be used to project a volume or xyz coordinates onto a flattened view of the dorsal cortex. The volume or points can be aggregated along the streamlines in different ways, e.g max, mean, median.\n",
+    "\n",
+    "In this example we will project the Allen Atlas dwi volume onto the dorsal cortex projection. We will take the maximum value along the streamline for each voxel on the surface"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "7b61405d-8589-4cb8-8afc-6f8a87dc4588",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from iblatlas.streamlines.utils import project_volume_onto_flatmap\n",
+    "\n",
+    "ba = AllenAtlas(25)\n",
+    "proj, fig, ax = project_volume_onto_flatmap(ba.image, res_um=25, aggr='max', plot=True, cmap='bone')\n",
+    "ax.set_axis_off()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "702b9abd-f97d-4a4e-967a-e9967e630293",
+   "metadata": {},
+   "source": [
+    "We can also do the same with xyz coordinates (we will use the same xyz coordinates for the region MOp and assign random values to each coordinate)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "19bf31f8-8f22-4f0e-91bb-7998f9379dc6",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from iblatlas.streamlines.utils import project_points_onto_flatmap, get_mask\n",
+    "\n",
+    "# Get a mask to plot our results on, options are 'image', 'annotation' or 'boundary\n",
+    "mask = get_mask('boundary')\n",
+    "\n",
+    "values = np.random.rand(xyz.shape[0])\n",
+    "\n",
+    "proj = project_points_onto_flatmap(xyz, values, aggr='mean', plot=False)\n",
+    "# Mask values with 0 with nan for display purposes\n",
+    "proj[proj == 0] = np.nan\n",
+    "\n",
+    "fig, ax = plt.subplots()\n",
+    "ax.imshow(mask, cmap='Greys')\n",
+    "ax.imshow(proj)\n",
+    "ax.set_axis_off()"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.11.9"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/iblatlas/atlas.py b/iblatlas/atlas.py
index bdbe713..e171b0e 100644
--- a/iblatlas/atlas.py
+++ b/iblatlas/atlas.py
@@ -915,7 +915,8 @@ def _take_remap(vol, ind, axis, mapping):
                 index = bc.xyz2i(np.array([coordinate] * 3))[axis]
             return _take(region_values, index, axis=axis)
 
-    def compute_boundaries(self, values):
+    @staticmethod
+    def compute_boundaries(values):
         """
         Compute the boundaries between regions on slice
         :param values:
@@ -1811,54 +1812,3 @@ def get_bc(res_um=10):
     bc = BrainCoordinates(nxyz=nxyz, xyz0=-bc.i2xyz(iorigin), dxyz=dxyz)
 
     return bc
-
-
-def _download_depth_files(file_name):
-    """
-    These files have been generated using this script
-    https://github.com/int-brain-lab/ibldevtools/blob/master/Mayo/flatmaps/2025-03-20_depths_from_streamlines.py
-    :param file_name:
-    :return:
-    """
-    file_path = BrainAtlas._get_cache_dir().joinpath('depths', file_name)
-    if not file_path.exists():
-        file_path.parent.mkdir(exist_ok=True, parents=True)
-        aws.s3_download_file(f'atlas/depths/{file_path.name}', file_path)
-
-    return np.load(file_path)
-
-
-def xyz_to_depth(xyz, res_um=25):
-    """
-    For a given xyz coordinates find the depth from the surface of the cortex. Note the lookup will only
-    work for xyz cooordinates that are in the Isocortex of the Allen volume. If coordinates outside of this
-    region are given then the depth is returned as nan.
-
-    Parameters
-    ----------
-    xyz : numpy.array
-        An (n, 3) array of Cartesian coordinates. The order is ML, AP, DV and coordinates should be given in meters
-        relative to bregma.
-
-    res_um : float or int
-        The resolution of the brain atlas to do the depth lookup
-
-    Returns
-    -------
-        numpy.array
-        The depths from the surface of the cortex for each cartesian coordinate. If the coordinate does not lie within
-        the Isocortex, depth value returned is nan
-    """
-
-    ind_flat = _download_depth_files(f'depths_ind_{res_um}.npy')
-    depths = _download_depth_files(f'depths_um_{res_um}.npy')
-    bc = get_bc(res_um=res_um)
-
-    ixyz = bc.xyz2i(xyz, mode='clip')
-    iravel = np.ravel_multi_index((ixyz[:, 1], ixyz[:, 0], ixyz[:, 2]), (bc.ny, bc.nx, bc.nz))
-    a, b = ismember(iravel, ind_flat)
-
-    lookup_depths = np.full(iravel.shape, np.nan, dtype=np.float32)
-    lookup_depths[a] = depths[b]
-
-    return lookup_depths
diff --git a/iblatlas/streamlines/__init__.py b/iblatlas/streamlines/__init__.py
new file mode 100644
index 0000000..1a1f548
--- /dev/null
+++ b/iblatlas/streamlines/__init__.py
@@ -0,0 +1,54 @@
+"""A package for working with Allen streamlines data.
+
+This package provides a way to use the Allen streamline data to
+1. Do a lookup of xyz position to depth from the surface of the cortex
+2. Project a volume or xyz positions along the streamlines to display them on a flattened dorsal cortex flatmap
+
+Streamlines are the paths that connect the surface of the isocortex to the white matter surface while following the
+curvature of these surfaces.
+
+For the xyz to depth position lookup the following files have been generated with 10, 25 and 50 um resolution
+
+     depths/
+     ├── depths_ind_{res}.npy
+     ├── depths_per_{res}.npy
+     ├── depths_um_{res}.npy
+
+- The depths_ind_{res}.npy contains the flattened index for all voxels contained within the isocortex of the Allen
+volume at the specified resolution. The ordering of the flattened index is AP, DV, ML
+- The depths_per_{res}.npy contains the depth in percentage from the surface of the cortex of each voxel in the isocortex.
+The percentage values were taken from the laplacian provided by the Allen institute. Voxels at the surface have a percentage
+ depth of 0 while those at the white matter boundary have a percentage depth of 1.
+- The depths_um_{res}.npy contains the depth in um from the surface of the cortex of each voxel in the isocortex. The depth
+values have been computed using the streamlines provided by the Allen institute. For each streamline the length of
+the streamline was computed and voxels assigned a depth based on the steps along the streamline. When multiple streamlines
+pass through a given voxel the median depth has been taken.
+
+The 25 and 50 um resolution datasets have been computed by subsampling the 10 um datasets.
+
+
+To enable the projection onto the dorsal flatmap the following files have been generated with 10, 25 and 50 um resolution
+
+     depths/
+     ├── paths_{res}.pqt
+     ├── dorsal_flatmap.npy
+     ├── dorsal_annotation.npy
+     ├── dorsal_image.npy
+
+
+- The paths_{res}.pqt files contain information about the voxels that correspond to each streamline. They have been
+derived from the dorsal_flatmap_paths_10.h5 paths provided by the Allen institute. We store the values as a pqt file and remove
+values that correspond to no streamline (paths=0) and duplicate voxels per streamline. As we have removed duplicate
+voxels per streamline this means that when taking the average or std along the streamline the behavior will differ
+from the Allen code.
+- The dorsal_flatmap.npy is a 2D numpy array with the dorsal cortex flatmap. The value in each voxel corresponds the
+streamline index that corresponds to that voxel.
+- The dorsal_annotation is 2D numpy array with the dorsal cortex flatmap. The value in each voxel corresponds to the
+region id on the surface of the cortex that corresponds to that voxel
+- The dorsal_image is a 2D numpy array with the dorsal cortex flatmap. The value in each voxel corresponds to the Allen
+image volume projected onto the flatmap display
+
+[1] Wang, Quanxin et al., “The Allen Mouse Brain Common Coordinate Framework: A 3D Reference Atlas”
+ Cell, Volume 181, Issue 4, 936 - 953.e20 May. 2020, https://doi.org/10.1016/j.cell.2020.04.007
+
+"""
\ No newline at end of file
diff --git a/iblatlas/streamlines/streamlines_scraping/prepare_depths.py b/iblatlas/streamlines/streamlines_scraping/prepare_depths.py
new file mode 100644
index 0000000..586cdc0
--- /dev/null
+++ b/iblatlas/streamlines/streamlines_scraping/prepare_depths.py
@@ -0,0 +1,173 @@
+"""
+This script was used to generate the following files
+- depths_um_10.npy
+- depths_um_25.npy
+- depths_um_50.npy
+
+These files are used by the function iblatlas.streamlines.utils.xyz_to_depth with the kwarg per=False to do a lookup of
+xyz position to depth in um from the surface of the cortex.
+
+The original streamlines file that we work off (streamlines_allen.npy) is available to download from here
+https://drive.google.com/drive/u/0/folders/1MbO3n2LAPEocDTDwRKCCJuVyNSvCcaG2
+
+(N.B we also have the IBL version of the streamlines but have decided to use the Allen one for now)
+
+There are cases where multiple streamlines pass through the same voxel. In these cases the voxel has been
+assigned the median depth value
+
+For voxels that don't have any streamlines that pass through them the values are interpolated using nearest neighbour
+
+We also subsample the depth volume to 25 um and 50 um resolution to generate the files for these resolution
+atlases.
+"""
+
+import matplotlib.pyplot as plt
+import numpy as np
+import pandas as pd
+from pathlib import Path
+from scipy.interpolate import NearestNDInterpolator
+from iblatlas.atlas import AllenAtlas
+
+file_path = Path('/Users/admin/int-brain-lab/flatmap')
+
+# -------------------------------------------------------------------------------------------------
+# Depth computation for 10 um
+# -------------------------------------------------------------------------------------------------
+
+ba = AllenAtlas(10)
+# Order is AP, ML, DV
+ba_shape = ba.image.shape
+
+# Convert label volume to cosmos mapping
+idx = ba.regions.acronym2index('Isocortex')[1][0][0]
+ba.label = ba.regions.mappings['Cosmos'][ba.label]
+
+# Find the voxels that lie in Isocortex AP, ML, DV
+xyz = np.where(ba.label == idx)
+# Compute the flattened idx for later
+xyz_flat = np.ravel_multi_index(xyz, ba_shape)
+
+streamlines = np.load(file_path.joinpath('streamlines_allen.npy'))
+# Reorder to match atlas image so AP, ML, DV
+streamlines = streamlines[:, :, [0, 2, 1]]
+
+# Compute cartesian coordinate of top of each streamline
+mlapdv_top = np.c_[ba.bc.i2x(streamlines[:, 0, 1]),
+                   ba.bc.i2y(streamlines[:, 0, 0]),
+                   ba.bc.i2z(streamlines[:, 0, 2])]
+# Compute cartesian coordinate of bottom of each streamline
+mlapdv_bottom = np.c_[ba.bc.i2x(streamlines[:, -1, 1]),
+                      ba.bc.i2y(streamlines[:, -1, 0]),
+                      ba.bc.i2z(streamlines[:, -1, 2])]
+
+# Compute the distance in um of the interval along each streamline
+steps_um = np.sqrt(np.sum((mlapdv_top - mlapdv_bottom) ** 2, axis=1)) * 1e6 / streamlines.shape[1]
+
+del mlapdv_top, mlapdv_bottom, ba
+
+# Get rid of streamlines that have no distance
+kp_strm = np.where(steps_um > 0)[0]
+streamlines = streamlines[kp_strm]
+steps_um = steps_um[kp_strm]
+
+del kp_strm
+
+# Make an array of depths along each streamline
+steps = np.arange(100) + 1
+steps = np.tile(steps.reshape(1, 100), (streamlines.shape[0], 1))
+depths = (steps * steps_um.reshape(steps_um.size, 1)).astype(np.float32)
+
+del steps, steps_um
+
+# Create a dataframe
+df = pd.DataFrame()
+df['inds'] = np.ravel_multi_index(
+    (streamlines[:, :, 0], streamlines[:, :, 1], streamlines[:, :, 2]), ba_shape).flatten()
+df['depths'] = depths.flatten()
+
+del depths, streamlines
+
+# Group by index in volume
+df_grp = df.groupby('inds')
+
+# Find the median value per voxel
+agg_val = df_grp.median()
+
+# Create a volume and set the values to the median values
+vol = np.zeros(np.prod(ba_shape), dtype=np.float32)
+vol[xyz_flat] = np.nan
+vol[agg_val.index.values] = agg_val.depths.values
+vol = vol.reshape(ba_shape)
+
+# Find the indices where we have data (only consider Isocortex)
+coords = np.where(np.bitwise_and(~np.isnan(vol), vol > 0))
+vals = vol[coords[0], coords[1], coords[2]]
+coords = np.array([coords[0], coords[1], coords[2]]).T
+
+# These are the indices in the Isocortex where we don't have values for therefore need to interpolate
+fill_coords = np.where(np.isnan(vol))
+fill_coords = np.array([fill_coords[0], fill_coords[1], fill_coords[2]]).T
+
+# Define a nearest neighbour interpolator
+interpolator = NearestNDInterpolator(coords, vals)
+fill_vals = interpolator(fill_coords)
+
+# Fill the volume with the interpolated values
+vol[fill_coords[:, 0], fill_coords[:, 1], fill_coords[:, 2]] = fill_vals
+vals = vol.flat[xyz_flat].astype(np.float32)
+
+np.save(file_path.joinpath('depths_um_10.npy'), vals)
+
+# -------------------------------------------------------------------------------------------------
+# Downsample to 25 and 50 um atlas
+# -------------------------------------------------------------------------------------------------
+
+# Load in 10um results
+ba = AllenAtlas(10)
+
+depths = np.load(file_path.joinpath('depths_um_10.npy'))
+
+# unravel so that they are in ixyz
+ixyz = np.unravel_index(xyz_flat, ba.image.shape)
+# Find the cartesian coordinates of all indices that we have depths for
+xyz = ba.bc.i2xyz(np.c_[ixyz[1], ixyz[0], ixyz[2]])
+
+del ba, ixyz
+
+for res in [25, 50]:
+
+    # Transform to 25 um atlas
+    ba = AllenAtlas(res)
+    df = pd.DataFrame()
+    df['depths'] = depths
+
+    # Find the index in the downsampled volume for the cartesian coordinates
+    ixyz = ba.bc.xyz2i(xyz)
+    df['ind'] = np.ravel_multi_index((ixyz[:, 1], ixyz[:, 0], ixyz[:, 2]), ba.image.shape)
+
+    # Group by index and find the median value
+    df_grp = df.groupby('ind')
+    df_med = df_grp.depths.median()
+    depths_res = df_med.values
+    xyz_flat_res = df_med.index.values
+
+    np.save(file_path.joinpath(f'depths_um_{res}.npy'), depths_res.astype(np.float32))
+
+    ibl_vol = np.zeros(np.prod(ba.image.shape))
+    ibl_vol[xyz_flat_res] = depths_res
+    ibl_vol = ibl_vol.reshape(ba.image.shape)
+    ibl_vol[ibl_vol == 0] = np.nan
+
+    ap = -2000 / 1e6
+    cmap = 'tab20'
+    vmin = 0
+    vmax = 2000
+    fig, ax = plt.subplots(figsize=(8, 5))
+    ax.imshow(np.moveaxis(ba.image[ba.bc.y2i(ap), :, :], 0, 1), extent=ba.extent(axis=1), cmap='gray')
+    im = ax.imshow(np.moveaxis(ibl_vol[ba.bc.y2i(ap), :, :], 0, 1),
+                   extent=ba.extent(axis=1), cmap=cmap, vmin=vmin, vmax=vmax)
+    ba.plot_cslice(ap, volume='boundary', ax=ax)
+    cbar = fig.colorbar(im, ax=ax, ticks=np.arange(vmin, vmax, 100))
+    cbar.set_ticklabels(np.arange(vmin, vmax, 100))
+    cbar.set_label('Depth from surface (um)')
+    plt.show()
diff --git a/iblatlas/streamlines/streamlines_scraping/prepare_laplacian.py b/iblatlas/streamlines/streamlines_scraping/prepare_laplacian.py
new file mode 100644
index 0000000..ba02147
--- /dev/null
+++ b/iblatlas/streamlines/streamlines_scraping/prepare_laplacian.py
@@ -0,0 +1,149 @@
+"""
+This script was used to generate the following files
+- depths_ind_10.npy
+- depths_ind_25.npy
+- depths_ind_50.npy
+- depths_per_10.npy
+- depths_per_25.npy
+- depths_per_50.npy
+
+These files are used by the function iblatlas.streamlines.utils.xyz_to_depth with the kwarg per=True to do a lookup of
+xyz position to percentage depth from the surface of the cortex.
+
+The original laplacian file that we work off (laplacian_10.nrrd) is provided by Allen and available to download
+from here  https://download.alleninstitute.org/informatics-archive/current-release/mouse_ccf/cortical_coordinates/ccf_2017/
+
+(N.B we also have the IBL version of the Laplacian but have decided to use the Allen one for now)
+
+Compared to the full laplacian volume provided by Allen we store
+- The flattened index of voxels contained in the Isocortex (flattened index uses IBL order i.e AP, ML, DV)
+    e.g depths_ind_10.npy
+- The percentage depth from surface of cortex for each flattened index
+    e.g depths_per_10.npy
+
+We also subsample the laplacian volume to 25 um and 50 um resolution to generate the files for these resolution
+atlases.
+"""
+
+import matplotlib.pyplot as plt
+import numpy as np
+import nrrd
+import pandas as pd
+from pathlib import Path
+
+from iblatlas.atlas import AllenAtlas
+
+file_path = Path('/Users/admin/Downloads/flatmap')
+
+# -------------------------------------------------------------------------------------------------
+# Prepare data for 10 um resolution
+# -------------------------------------------------------------------------------------------------
+
+# Load in the 10 um Allen Atlas
+ba = AllenAtlas(10)
+ba_shape = ba.image.shape
+# Translate the label volume to cosmos mapping
+idx = ba.regions.acronym2index('Isocortex')[1][0][0]
+ba.label = ba.regions.mappings['Cosmos'][ba.label]
+# Find the voxels that are contained within the isocortex, order AP, ML, DV
+xyz = np.where(ba.label == idx)
+# Compute the flattened idx that we will save
+xyz_flat = np.ravel_multi_index(xyz, ba_shape).astype(np.int32)
+
+del ba
+
+# Load in the 10 um resolution laplacian, orientation AP, DV, ML
+laplacian, _ = nrrd.read(file_path.joinpath('laplacian_10.nrrd'))
+vals = laplacian[xyz[0], xyz[2], xyz[1]]
+vals = vals.flatten().astype(np.float32)
+
+del laplacian, xyz
+
+# We just save the values for the voxels contained in isocortex
+np.save(file_path.joinpath('depths_ind_10.npy'), xyz_flat)
+np.save(file_path.joinpath('depths_per_10.npy'), vals)
+
+# Check it makes sense
+vol = np.full(np.prod(ba_shape), np.nan)
+vol[xyz_flat] = vals
+vol = vol.reshape(ba_shape)
+
+ba = AllenAtlas(10)
+ap = -2000 / 1e6
+fig, ax = plt.subplots(figsize=(8, 5))
+ax.imshow(np.moveaxis(ba.image[ba.bc.y2i(ap), :, :], 0, 1),
+          extent=ba.extent(axis=1), cmap='gray')
+im = ax.imshow(np.moveaxis(vol[ba.bc.y2i(ap), :, :], 0, 1),
+               extent=ba.extent(axis=1))
+plt.show()
+
+del vol, xyz_flat, vals, ba
+
+# -------------------------------------------------------------------------------------------------
+# Downsample to 25 um and 50 um
+# -------------------------------------------------------------------------------------------------
+
+ba = AllenAtlas(10)
+xyz_flat = np.load(file_path.joinpath('depths_ind_10.npy'))
+vals = np.load(file_path.joinpath('depths_per_10.npy'))
+# unravel so that they are in ixyz
+ixyz = np.unravel_index(xyz_flat, ba.image.shape)
+# Find the cartesian coordinates of all indices that we have depths for
+xyz = ba.bc.i2xyz(np.c_[ixyz[1], ixyz[0], ixyz[2]])
+
+del ba, ixyz
+
+for res in [25, 50]:
+
+    # Transform to 25 um atlas
+    ba = AllenAtlas(res)
+    df = pd.DataFrame()
+    df['vals'] = vals
+
+    # Find the index in the down-sampled volume for the cartesian coordinates
+    ixyz = ba.bc.xyz2i(xyz)
+    df['ind'] = np.ravel_multi_index((ixyz[:, 1], ixyz[:, 0], ixyz[:, 2]), ba.image.shape)
+
+    # Group by index and find the median value
+    df_grp = df.groupby('ind')
+    df_med = df_grp.vals.median()
+    vals_res = df_med.values
+    xyz_flat_res = df_med.index.values
+    np.save(file_path.joinpath(f'depths_ind_{res}.npy'), xyz_flat_res.astype(np.int32))
+    np.save(file_path.joinpath(f'depths_per_{res}.npy'), vals_res.astype(np.float32))
+
+    vol = np.full(np.prod(ba.image.shape), np.nan)
+    vol[xyz_flat_res] = vals_res
+    vol = vol.reshape(ba.image.shape)
+
+    ap = -2000 / 1e6
+    fig, ax = plt.subplots(figsize=(8, 5))
+    ax.imshow(np.moveaxis(ba.image[ba.bc.y2i(ap), :, :], 0, 1), extent=ba.extent(axis=1), cmap='gray')
+    im = ax.imshow(np.moveaxis(vol[ba.bc.y2i(ap), :, :], 0, 1),
+                   extent=ba.extent(axis=1))
+    plt.show()
+
+
+# -------------------------------------------------------------------------------------------------
+# Check that the laplacian looks okay across the different resolutions
+# -------------------------------------------------------------------------------------------------
+
+fig, axs = plt.subplots(1, 3, sharex=True, sharey=True)
+ap = -2000 / 1e6
+
+for i, res in enumerate([10, 25, 50]):
+    xyz_flat = np.load(file_path.joinpath(f'depths_ind_{res}.npy'))
+    vals = np.load(file_path.joinpath(f'depths_per_{res}.npy'))
+    ba = AllenAtlas(res)
+    ba_shape = ba.image.shape
+
+    vol = np.full(np.prod(ba_shape), np.nan)
+    vol[xyz_flat] = vals
+    vol = vol.reshape(ba_shape)
+
+    ax = axs[i]
+    ax.imshow(np.moveaxis(ba.image[ba.bc.y2i(ap), :, :], 0, 1), extent=ba.extent(axis=1), cmap='gray')
+    im = ax.imshow(np.moveaxis(vol[ba.bc.y2i(ap), :, :], 0, 1),
+                   extent=ba.extent(axis=1))
+
+plt.show()
diff --git a/iblatlas/streamlines/streamlines_scraping/prepare_streamlines.py b/iblatlas/streamlines/streamlines_scraping/prepare_streamlines.py
new file mode 100644
index 0000000..8f0978b
--- /dev/null
+++ b/iblatlas/streamlines/streamlines_scraping/prepare_streamlines.py
@@ -0,0 +1,140 @@
+"""
+This script was used to generate the following files
+- paths_10.pqt
+- paths_25.pqt
+- paths_50.pqt
+- dorsal_flatmap.npy
+
+These files are used by the module iblatlas.streamlines to project 3D volumes or points onto the dorsal cortex
+flatmap using the streamline paths
+
+The original paths file that we work from (dorsal_flatmap_paths_10.h5) is provided by Allen and available to download
+from here  https://download.alleninstitute.org/informatics-archive/current-release/mouse_ccf/cortical_coordinates/ccf_2017/
+
+Compared to Allen the paths are stored in a pandas Dataframes that has columns
+- paths - indicates the streamline index
+- lookup - indicates the voxel in the 3D Allen volume for each path
+
+Compared to the streamlines provided by Allen we remove duplicate voxels per streamline, e.g compare
+paths[677974] from Allen to paths[677974] from IBL
+
+We have also subsampled the streamlines to the 25 um and 50 um volumes
+"""
+
+# -------------------------------------------------------------------------------------------------
+# Process for 10 um
+# -------------------------------------------------------------------------------------------------
+
+import h5py
+import pandas as pd
+import numpy as np
+from iblatlas.atlas import AllenAtlas, get_bc
+from iblatlas.streamlines.utils import project_volume_onto_flatmap
+from pathlib import Path
+file_path = Path('/Users/admin/Downloads/flatmap')
+
+ba = AllenAtlas(10)
+ba_shape = ba.image.shape
+
+del ba
+
+f1 = h5py.File(file_path.joinpath('dorsal_flatmap_paths_10.h5'), 'r+')
+paths = f1['paths'][:]
+f1.close()
+
+# Paths from Allen are in order AP, DV, ML
+paths = np.unravel_index(paths, (1320, 800, 1140))
+# Reorder so that the indices are in IBL order AP, ML, DV
+paths = np.ravel_multi_index((paths[0], paths[2], paths[1]), ba_shape)
+path_shape = paths.shape
+
+# When the value = 0 in a path it is locations where the streamline has ended
+# Find the longest streamline
+
+path_df = pd.DataFrame()
+paths_idx = np.repeat(np.arange(paths.shape[0]), paths.shape[1])
+path_df['lookup'] = paths.flatten().astype(np.int32)
+path_df['paths'] = paths_idx.astype(np.int32)
+
+# Remove the rows where the lookup = 0. These indicate locations where there is no voxel for that depth
+# e.g shorter streamlines
+kp_idx = paths.flat != 0
+path_df = path_df[kp_idx]
+# Remove duplicate voxels within the same streamline
+path_df = path_df.drop_duplicates().reset_index(drop=True)
+
+path_df.to_parquet(file_path.joinpath('paths_10.pqt'))
+
+
+# -------------------------------------------------------------------------------------------------
+# Downsample for 25 um and 50 um
+# -------------------------------------------------------------------------------------------------
+
+res = 50  # 25 or 50
+
+f1 = h5py.File(file_path.joinpath('dorsal_flatmap_paths_10.h5'), 'r+')
+paths = f1['paths'][:]
+f1.close()
+
+bc = get_bc(10)
+bc_res = get_bc(res)
+
+# Paths are in order AP, DV, ML
+path_shape = paths.shape
+paths = paths.ravel()
+paths = np.unravel_index(paths, (1320, 800, 1140))
+
+xyz = bc.i2xyz(np.array([paths[2], paths[0], paths[1]]).T)
+del paths
+ixyz = bc_res.xyz2i(xyz)
+del xyz
+ixyz = np.ravel_multi_index((ixyz[:, 1], ixyz[:, 0], ixyz[:, 2]), (bc_res.ny, bc_res.nx, bc_res.nz))
+
+paths = ixyz.reshape(path_shape)
+
+del ixyz
+
+# When the value = 0 in a path it is locations where the streamline has ended
+# Find the longest streamline
+path_df = pd.DataFrame()
+paths_idx = np.repeat(np.arange(paths.shape[0]), paths.shape[1])
+path_df['lookup'] = paths.flatten().astype(np.int32)
+path_df['paths'] = paths_idx.astype(np.int32)
+
+# Remove the rows where the lookup = 0. These indicate locations where there is no voxel for that depth
+# e.g shorter streamlines
+kp_idx = paths.flat != 0
+path_df = path_df[kp_idx]
+# Remove duplicate voxels within the same streamline
+path_df = path_df.drop_duplicates()
+
+path_df = path_df.drop_duplicates().reset_index(drop=True)
+
+path_df.to_parquet(file_path.joinpath(f'paths_{res}.pqt'))
+
+
+# -------------------------------------------------------------------------------------------------
+# Extract the flatmap file
+# -------------------------------------------------------------------------------------------------
+
+f1 = h5py.File(file_path.joinpath('dorsal_flatmap_paths_10.h5'), 'r+')
+flatmap = f1['flat'][:]
+f1.close()
+
+np.save(file_path.joinpath('dorsal_flatmap.npy'), flatmap)
+
+
+# -------------------------------------------------------------------------------------------------
+# Get the image and annotation files
+# -------------------------------------------------------------------------------------------------
+ba = AllenAtlas(10)
+
+# Annotation
+img = project_volume_onto_flatmap(ba.label, res_um=10, aggr='first', plot=False)
+img[0, 0] = 0
+np.save(file_path.joinpath('dorsal_annotation.npy'), img.astype(np.int16))
+
+# Image
+img = project_volume_onto_flatmap(ba.image, res_um=10, aggr='first', plot=False)
+img[0, 0] = 0
+np.save(file_path.joinpath('dorsal_image.npy'), img.astype(np.float16))
diff --git a/iblatlas/streamlines/utils.py b/iblatlas/streamlines/utils.py
new file mode 100644
index 0000000..98c649c
--- /dev/null
+++ b/iblatlas/streamlines/utils.py
@@ -0,0 +1,312 @@
+import pandas as pd
+import numpy as np
+from iblutil.numerical import ismember
+import matplotlib.pyplot as plt
+from iblatlas.atlas import get_bc, BrainAtlas, aws, AllenAtlas
+from iblatlas.regions import BrainRegions
+
+
+def _download_depth_files(file_name):
+    """
+    Download and return path to relevant file
+    :param file_name:
+    :return:
+    """
+    file_path = BrainAtlas._get_cache_dir().joinpath('depths', file_name)
+    if not file_path.exists():
+        file_path.parent.mkdir(exist_ok=True, parents=True)
+        aws.s3_download_file(f'atlas/depths/{file_path.name}', file_path)
+
+    return file_path
+
+
+def xyz_to_depth(xyz, per=True, res_um=25):
+    """
+    For a given xyz coordinates return the depth from the surface of the cortex. The depth is returned
+    as a percentage if per=True and in um if per=False. Note the lookup will only work for xyz cooordinates
+    that are in the Isocortex of the Allen volume. If coordinates outside of this region are given then
+    the depth is returned as nan.
+
+    Parameters
+    ----------
+    xyz : numpy.array
+        An (n, 3) array of Cartesian coordinates. The order is ML, AP, DV and coordinates should be given in meters
+        relative to bregma.
+
+    per : bool
+        Whether to do the lookup in percentage from the surface of the cortex or depth in um from the surface of the cortex.
+
+    res_um : float or int
+        The resolution of the brain atlas to do the depth lookup
+
+    Returns
+    -------
+        numpy.array
+        The depths from the surface of the cortex for each cartesian coordinate. If the coordinate does not lie within
+        the Isocortex, depth value returned is nan
+    """
+
+    ind_flat = np.load(_download_depth_files(f'depths_ind_{res_um}.npy'))
+    depth_file = f'depths_per_{res_um}.npy' if per else f'depths_um_{res_um}.npy'
+    depths = np.load(_download_depth_files(depth_file))
+    bc = get_bc(res_um=res_um)
+
+    ixyz = bc.xyz2i(xyz, mode='clip')
+    iravel = np.ravel_multi_index((ixyz[:, 1], ixyz[:, 0], ixyz[:, 2]), (bc.ny, bc.nx, bc.nz))
+    a, b = ismember(iravel, ind_flat)
+
+    lookup_depths = np.full(iravel.shape, np.nan, dtype=np.float32)
+    lookup_depths[a] = depths[b]
+
+    return lookup_depths
+
+
+def get_mask(volume='annotation', br=None):
+    """
+    Generate a mask to plot results onto
+
+    Parameters:
+    -----------
+    volume : str, optional
+        The type of volume to project. Options are:
+        - 'image': Projects the anatomical image using max intensity.
+        - 'annotation': Projects the labeled regions and maps them to RGB colors.
+        - 'boundary': Projects labeled regions and extracts anatomical boundaries.
+        Default is 'annotation'.
+
+    br : BrainRegions, optional
+        An instance of the BrainRegions If None, a default BrainRegions is initialized.
+
+    Returns:
+    --------
+    img : np.ndarray
+        The resulting 2D flatmap projection image, either grayscale, RGB, or binary mask,
+        depending on the selected volume type.
+    """
+
+    br = br or BrainRegions()
+    if volume == 'image':
+        img = np.load(_download_depth_files('dorsal_image.npy'))
+    elif volume == 'annotation':
+        img = np.load(_download_depth_files('dorsal_annotation.npy'))
+        img = br.rgb[img]
+    elif volume == 'boundary':
+        img = np.load(_download_depth_files('dorsal_annotation.npy'))
+        img = AllenAtlas.compute_boundaries(img)
+
+    return img
+
+
+def validate_aggr(aggr: str) -> None:
+    """
+    Validates if the provided aggregation type is valid.
+
+    Parameters:
+    ----------
+    aggr : str
+        The aggregation method to validate (e.g., 'mean', 'sum').
+
+    Raises:
+    ------
+    AssertionError
+        If the aggregation type is not one of the allowed values.
+    """
+    poss_aggrs = ['sum', 'count', 'mean', 'std', 'median', 'min', 'max', 'first', 'last']
+    assert aggr in poss_aggrs, f"Aggregation must be one of {poss_aggrs}."
+
+
+def project_volume_onto_flatmap(vol: np.ndarray, res_um: int = 25, aggr: str = 'mean', plot: bool = True,
+                                cmap: str = 'viridis', clevels: tuple = None, ax: plt.Axes = None) -> np.ndarray:
+    """
+    Projects a 3D volume onto a 2D flatmap by aggregating values along streamline paths.
+
+    Parameters:
+    ----------
+    vol : np.ndarray
+        A 3D numpy array representing the volume data to be projected.
+
+    res_um : int
+        The resolution of the volume. Must be one of 10, 25 or 50.
+
+    aggr : str, optional
+        The aggregation method ('sum', 'count', 'mean', etc.), default is 'mean'.
+
+    plot : bool, optional
+        Whether to plot the resulting projection, default is True.
+
+    cmap : str, optional
+        The colormap to use for the plot, default is 'viridis'.
+
+    clevels : tuple, optional
+        The color limits to use for the plot, default is None.
+
+    ax : matplotlib.axes.Axes, optional
+        The axes on which to plot, default is None.
+
+    Returns:
+    -------
+    np.ndarray
+        The projected 2D array onto the flatmap.
+    matplotlib.figure.Figure
+        Matplotlib figure object if plot=True, otherwise None.
+    matplotlib.axes.Axes
+        Matplotlib axes object if plot=True, otherwise None.
+    """
+    bc = get_bc(res_um)
+    assert vol.shape == (bc.ny, bc.nx, bc.nz), f"Volume does not have expected shape of {(bc.ny, bc.nx, bc.nz)}"
+
+    # Validate the aggregation type
+    validate_aggr(aggr)
+
+    # Load the streamline paths
+    path_df = pd.read_parquet(_download_depth_files(f'paths_{res_um}.pqt'))
+
+    # Extract values from the volume using the path lookup
+    path_df['vals'] = vol.flat[path_df['lookup'].values]
+
+    # Aggregate the values along each path
+    flat_df = path_df.groupby('paths').vals.agg(aggr)
+
+    # Project the aggregated data onto the flatmap
+    return _project_onto_flatmap(flat_df, plot=plot, cmap=cmap, clevels=clevels, ax=ax)
+
+
+def project_points_onto_flatmap(xyz: np.ndarray, values: np.ndarray, res_um: int = 25, aggr: str = 'mean', plot: bool = True,
+                                cmap: str = 'viridis', clevels: tuple = None, ax: plt.Axes = None) -> np.ndarray:
+    """
+    Projects 3D points with associated values onto a 2D flatmap.
+
+    Parameters:
+    ----------
+    xyz : np.ndarray
+        An array containing xyz coordinates of the points to be projected. xyz values should be given in metres
+
+    values : np.ndarray
+        A 1D array of values to associate with the points.
+
+    res : int
+        The resolution to load the corresponding streamline paths.
+
+    aggr : str, optional
+        The aggregation method ('sum', 'count', 'mean', etc.), default is 'mean'.
+
+    plot : bool, optional
+        Whether to plot the resulting projection, default is True.
+
+    cmap : str, optional
+        The colormap to use for the plot, default is 'viridis'.
+
+    clevels : tuple, optional
+        The color limits to use for the plot, default is None.
+
+    ax : matplotlib.axes.Axes, optional
+        The axes on which to plot, default is None.
+
+    Returns:
+    -------
+    np.ndarray
+        The projected 2D array onto the flatmap.
+    matplotlib.figure.Figure
+        Matplotlib figure object if plot=True, otherwise None.
+    matplotlib.axes.Axes
+        Matplotlib axes object if plot=True, otherwise None.
+    """
+    # Ensure that xyz and values have matching dimensions
+    assert xyz.shape[0] == values.size, "xyz must have the same number of rows as values."
+
+    # Validate the aggregation type
+    validate_aggr(aggr)
+
+    # Get the boundary coordinates for the given resolution
+    bc = get_bc(res_um)
+
+    # Convert coordinates xyz to indices in volume
+    ixyz = bc.xyz2i(xyz, mode='clip')
+
+    # Create DataFrame of values and their corresponding flattened indices
+    val_df = pd.DataFrame()
+    val_df['vals'] = values
+    val_df['lookup'] = np.ravel_multi_index((ixyz[:, 1], ixyz[:, 0], ixyz[:, 2]), (bc.ny, bc.nx, bc.nz))
+
+    # Remove entries with invalid lookups
+    val_df = val_df[val_df['lookup'] != 0]
+
+    # Load streamline paths
+    path_df = pd.read_parquet(_download_depth_files(f'paths_{res_um}.pqt'))
+
+    # Restrict dataframes to overlapping locations
+    vals_in_paths, _ = ismember(val_df['lookup'].values, path_df['lookup'].values)
+    val_df = val_df[vals_in_paths]
+
+    # Keep only paths that match the val lookups
+    paths_in_vals, _ = ismember(path_df['lookup'].values, val_df['lookup'].values)
+    path_df = path_df[paths_in_vals]
+
+    # Merge path data with values and aggregate by path
+    flat_df = path_df.merge(val_df, on='lookup', how='left').groupby('paths').vals.agg(aggr)
+
+    # Project the aggregated data onto the flatmap
+    return _project_onto_flatmap(flat_df, plot=plot, cmap=cmap, clevels=clevels, ax=ax)
+
+
+def _project_onto_flatmap(flat_df: pd.Series, plot: bool = True, cmap: str = 'viridis',
+                          clevels: tuple = None, ax: plt.Axes = None) -> np.ndarray:
+    """
+    Function to project aggregated data onto a 2D flatmap.
+
+    Parameters:
+    ----------
+    flat_df : pd.Series
+        The data to project onto the flatmap.
+
+    plot : bool, optional
+        Whether to plot the resulting projection, default is True.
+
+    cmap : str, optional
+        The colormap to use for the plot, default is 'viridis'.
+
+    clevels : tuple, optional
+        The color limits to use for the plot, default is None.
+
+    ax : matplotlib.axes.Axes, optional
+        The axes on which to plot, default is None.
+
+    Returns:
+    -------
+    np.ndarray
+        The projected 2D array onto the flatmap.
+    plt.figure
+        Matplotlib figure object if plot=True, otherwise None.
+    plt.axes
+        Matplotlib axes object if plot=True, otherwise None.
+    """
+    # Load the flatmap
+    flatmap = np.load(_download_depth_files('dorsal_flatmap.npy'))
+
+    # Find the indices in the flatmap corresponding to the projected data
+    _, b = ismember(flat_df.index, flatmap.flat)
+
+    # Initialize the projection array with zeros
+    proj = np.zeros(np.prod(flatmap.shape))
+
+    # Assign the values from flat_df to the projection array
+    proj[b] = flat_df.values
+
+    # Reshape the 1D projection into the 2D flatmap shape
+    proj = proj.reshape(flatmap.shape)
+
+    # Plot the result if requested
+    if plot:
+        if ax:
+            fig = ax.get_figure()
+        else:
+            fig, ax = plt.subplots()
+
+        if clevels is None:
+            clevels = (np.nanmin(proj), np.nanmax(proj))
+
+        ax.imshow(proj, cmap=cmap, vmin=clevels[0], vmax=clevels[1])
+
+        return proj, fig, ax
+    else:
+        return proj
diff --git a/iblatlas/tests/test_atlas.py b/iblatlas/tests/test_atlas.py
index 88c4c7b..c906471 100644
--- a/iblatlas/tests/test_atlas.py
+++ b/iblatlas/tests/test_atlas.py
@@ -6,7 +6,8 @@
 
 import iblatlas.atlas
 from iblatlas.atlas import (BrainCoordinates, Trajectory, BrainAtlas,
-                            Insertion, ALLEN_CCF_LANDMARKS_MLAPDV_UM, AllenAtlas, get_bc, xyz_to_depth)
+                            Insertion, ALLEN_CCF_LANDMARKS_MLAPDV_UM, AllenAtlas, get_bc)
+from iblatlas.streamlines.utils import xyz_to_depth
 from iblatlas.regions import BrainRegions, FranklinPaxinosRegions
 from iblatlas.plots import prepare_lr_data, reorder_data
 from iblutil.numerical import ismember
@@ -753,7 +754,7 @@ def test_conversions(self):
 
 class TestDepths(unittest.TestCase):
 
-    @mock.patch('iblatlas.atlas._download_depth_files',
+    @mock.patch('iblatlas.streamlines.utils.np.load',
                 side_effect=[np.array([35838117, 41837813]), np.array([400, 600])])
     def test_depth_lookup(self, mock_file):
         bc = get_bc(25)
@@ -761,7 +762,7 @@ def test_depth_lookup(self, mock_file):
         idx = np.array([[245, 274, 37], [286, 327, 53], [225, 287, 170]])  # ap, ml, dv
         xyz = np.c_[bc.i2x(idx[:, 1]), bc.i2y(idx[:, 0]), bc.i2z(idx[:, 2])]
 
-        depths = xyz_to_depth(xyz)
+        depths = xyz_to_depth(xyz, per=False)
 
         assert np.array_equal(depths, np.array([400, 600, np.nan]), equal_nan=True)
 
diff --git a/pyproject.toml b/pyproject.toml
index f675a5f..8e2213d 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -45,4 +45,5 @@ include-package-data = true
 
 [tool.setuptools.packages.find]
 # exclude packages matching these glob patterns
-exclude = ["iblatlas.tests*", "examples*", "iblatlas.genomics.gene_expression_scrapping"]
+exclude = ["iblatlas.tests*", "examples*", "iblatlas.genomics.gene_expression_scrapping",
+    "iblatlas.streamlines.streamlines_scraping"]