diff --git a/examples/save_load_hcs_data.py b/examples/save_load_hcs_data.py
new file mode 100644
index 0000000..4185148
--- /dev/null
+++ b/examples/save_load_hcs_data.py
@@ -0,0 +1,54 @@
+import napari
+import os
+import shutil
+import squidpy as sq
+from ngff_tables_prototype.writer import write_spatial_anndata
+from ngff_tables_prototype.reader import load_to_napari_viewer
+import numpy as np
+import pandas as pd
+from pathlib import Path
+import anndata as ad
+
+# Instructions
+# Download files from here:
+# https://drive.google.com/drive/folders/1TWMcKSQuq_ZEkTPPCLOO4sE6oMlhSGu_?usp=sharing
+# Place it in the ome-ngff-tables-prototype folder and run the script from
+# the example folder
+# Unzip the OME-Zarr file
+output_fpath = Path("../HCS_Example_Data/UZH_2x2_HCS.ome.zarr")
+feature_path = Path("../HCS_Example_Data/FeatureData")
+
+def write_hcs_table() -> None:
+    # Loop over the 4 sites
+    # Check that data was downloaded
+    assert len(os.listdir(feature_path)) == 4, 'Feature data should contain \
+                                            for files. Check instructions above'
+    assert os.path.isdir(output_fpath), 'OME Zarr folder does not exist. Was it unzipped?'
+
+    # TODO: Check if feature data has already been written to the OME-Zarr file
+
+    # Loop through all 4 fields of view,
+    # load the feature data & save to OME-Zarr file
+    for site in range(4):
+        df = pd.read_csv(feature_path / f'HCS_FeatureData_Site{site}.csv')
+        df_x = df.iloc[:, :-2]
+        features_ad = ad.AnnData(df_x.loc[:, df_x.columns != 'Label'])
+        features_ad.obs = pd.concat([df['Label'], df.iloc[:, -2:]], axis = 1)
+        write_spatial_anndata(
+            file_path=output_fpath / 'B' / '03' / str(site),
+            tables_adata=features_ad,
+            tables_region="labels/label_image",
+            tables_instance_key="Label",
+        )
+
+    return
+
+
+if __name__ == "__main__":
+    write_hcs_table()
+    print(output_fpath)
+    viewer = load_to_napari_viewer(
+        file_path=output_fpath / 'B' / '03' / '0',
+        groups=["labels/label_image", "tables/regions_table"],
+    )
+    napari.run()
diff --git a/setup.cfg b/setup.cfg
index 3e74782..296584a 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -38,6 +38,8 @@ install_requires =
     napari-ome-zarr
     ome-zarr
     pandas
+    matplotlib
+    squidpy
 
 
 [options.packages.find]
diff --git a/src/ngff_tables_prototype/writer.py b/src/ngff_tables_prototype/writer.py
index d0437fa..ed96c9b 100644
--- a/src/ngff_tables_prototype/writer.py
+++ b/src/ngff_tables_prototype/writer.py
@@ -127,6 +127,7 @@ def write_spatial_anndata(
     # label group
     label_image: Optional[np.ndarray] = None,
     label_name: str = "label_image",
+    label_axes: list = ['y', 'x'],
     # table group
     tables_adata: Optional[AnnData] = None,
     tables_region: Optional[Union[str, List[str]]] = None,
@@ -158,6 +159,8 @@ def write_spatial_anndata(
         The name of the label image. See ome-zarr-py for details.
     label_image : Union[str, List[str]]
         The label image (i.e. segmentation mask). See ome-zarr-py for details.
+    label_axes: list
+        The labels for the label image axes. See ome-zarr-py for details. Defaults to ['y', 'x'] for 2D images
     tables_adata:
         The :class:`anndata.AnnData` table with gene expression and annotations.
     tables_region
@@ -241,7 +244,7 @@ def write_spatial_anndata(
     if label_image is not None:
         # i.e. segmentation raster masks
         # the function write labels will create the group labels, so we pass the root
-        write_labels(label_image, group=root, name=label_name)
+        write_labels(label_image, group=root, name=label_name, axes=label_axes)
 
     if tables_adata is not None:
         # e.g. expression table