Repository: czbiohub-sf/iohub (language: python)

Issue:
zarr store written by iohub.writer doesn't seem to be compliant with OME-HCS format
Using napari-ome-zarr plugin to view the dataset written by the writer module causes an error that suggests missing metadata:
` napari --plugin napari-ome-zarr /hpc/projects/CompMicro/rawdata/hummingbird/Janie/2022_03_15_orgs_nuc_mem_63x_04NA/all_21_3.zarr`

causes the following error:
```
ERROR Failed to load Row_0/Col_271/0/0
Traceback (most recent call last):
  File "/hpc/user_apps/comp_micro/conda_envs/recorder/lib/python3.9/site-packages/fsspec/mapping.py", line 137, in __getitem__
    result = self.fs.cat(k)
  File "/hpc/user_apps/comp_micro/conda_envs/recorder/lib/python3.9/site-packages/fsspec/spec.py", line 811, in cat
    return self.cat_file(paths[0], **kwargs)
  File "/hpc/user_apps/comp_micro/conda_envs/recorder/lib/python3.9/site-packages/fsspec/spec.py", line 710, in cat_file
    with self.open(path, "rb", **kwargs) as f:
  File "/hpc/user_apps/comp_micro/conda_envs/recorder/lib/python3.9/site-packages/fsspec/spec.py", line 1094, in open
    f = self._open(
  File "/hpc/user_apps/comp_micro/conda_envs/recorder/lib/python3.9/site-packages/fsspec/implementations/local.py", line 175, in _open
    return LocalFileOpener(path, mode, fs=self, **kwargs)
  File "/hpc/user_apps/comp_micro/conda_envs/recorder/lib/python3.9/site-packages/fsspec/implementations/local.py", line 273, in __init__
    self._open()
  File "/hpc/user_apps/comp_micro/conda_envs/recorder/lib/python3.9/site-packages/fsspec/implementations/local.py", line 278, in _open
    self.f = open(self.path, mode=self.mode)
FileNotFoundError: [Errno 2] No such file or directory: '/hpc/projects/CompMicro/rawdata/hummingbird/Janie/2022_03_15_orgs_nuc_mem_63x_04NA/all_21_3.zarr/Row_0/Col_271/0/0/.zarray'

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/hpc/user_apps/comp_micro/conda_envs/recorder/lib/python3.9/site-packages/zarr/storage.py", line 1376, in __getitem__
    return self.map[key]
  File "/hpc/user_apps/comp_micro/conda_envs/recorder/lib/python3.9/site-packages/fsspec/mapping.py", line 141, in __getitem__
    raise KeyError(key)
KeyError: 'Row_0/Col_271/0/0/.zarray'
```

This .zarray json file is found in the reference data. It can be viewed using:
`napari --plugin napari-ome-zarr /local/scratch/groups/cmanalysis.grp/zarr_ref_datasets/2551.zarr/` 

I copied the reference dataset with the following aws cli command:
`aws --no-sign-request --endpoint-url=https://uk1s3.embassy.ebi.ac.uk/ s3 cp --recursive --exclude=".git/*" s3://idr/zarr/v0.1/plates/2551.zarr /local/scratch/groups/cmanalysis.grp/zarr_ref_datasets/2551.zarr`


Discussion / hints:
`iohub.writer` offers optional support to write HCS-compliant stores that needs to be explicitly specified:
https://github.com/czbiohub/iohub/blob/2564737c400d38c0856b63f51680e5c222e083ff/iohub/writer.py#L22-L26

@mattersoflight Do you know how is this dataset generated? If it's a `recOrder` acquisition then it does not use the HCS option.
In addition to this, it does seem that the HCS writer writes to an old spec which may not be supported by newer versions of `napari-ome-zarr`. (See also https://github.com/mehta-lab/waveorder/issues/94).

While `waveorder` may be constrained by backwards-compatibility issues, an upgrade to the newest spec (v0.4 or v0.5) should be implemented in `iohub`.

Fix patch:
diff --git a/examples/Writer_Example.ipynb b/examples/Writer_Example.ipynb
deleted file mode 100644
index 25c848ff..00000000
--- a/examples/Writer_Example.ipynb
+++ /dev/null
@@ -1,883 +0,0 @@
-{
- "cells": [
-  {
-   "cell_type": "code",
-   "execution_count": 1,
-   "id": "292ecf7c",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "import os\n",
-    "import sys\n",
-    "p = os.path.abspath('../..')\n",
-    "if p not in sys.path:\n",
-    "    sys.path.append(p)\n",
-    "\n",
-    "from waveorder.io.writer import WaveorderWriter\n",
-    "import numpy as np\n",
-    "import zarr"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "id": "2d4a3894",
-   "metadata": {},
-   "source": [
-    "# Writing data with no HCS Specification"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "id": "2bceef7c",
-   "metadata": {},
-   "source": [
-    "### Initialize the writer\n",
-    "\n",
-    "This will instantiate the writer class, no zarr stores, subfolders, etc have been created yet.  This step allows you to define which type of writer you want and let it know your directory to save the data"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 4,
-   "id": "e73b0ce8",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "writer = WaveorderWriter('/Users/cameron.foltz/Desktop/Test', hcs=False, hcs_meta=None, verbose=True)"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "id": "af804973",
-   "metadata": {},
-   "source": [
-    "### Set Position and Create Zarr\n",
-    "\n",
-    "'create_zarr_root()' will create an empty zarr store within the save directory.  No arrays have been added here.  Name of the store must be input by the user."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 5,
-   "id": "5196b4a0",
-   "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "Creating new zarr store at /Users/cameron.foltz/Desktop/Test/waveOrder_Writer_Example.zarr\n"
-     ]
-    }
-   ],
-   "source": [
-    "writer.create_zarr_root('waveOrder_Writer_Example.zarr')"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "id": "d35ea316",
-   "metadata": {},
-   "source": [
-    "### Create some fake data to save later"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 6,
-   "id": "5cf8cb26",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "P = 10\n",
-    "T = 3\n",
-    "C = 2\n",
-    "Z = 65\n",
-    "Y = 256\n",
-    "X = 256\n",
-    "\n",
-    "data = np.random.rand(P, T, C, Z, Y, X)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 7,
-   "id": "e70549cb",
-   "metadata": {},
-   "outputs": [
-    {
-     "data": {
-      "text/plain": [
-       "(10, 3, 2, 65, 256, 256)"
-      ]
-     },
-     "execution_count": 7,
-     "metadata": {},
-     "output_type": "execute_result"
-    }
-   ],
-   "source": [
-    "np.shape(data)"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "id": "99679028",
-   "metadata": {},
-   "source": [
-    "### initialize the array which sets the data saving parameters\n",
-    "\n",
-    "Here is where you will initialize the array in your current writer.store .  data_shape should match the (T, C, Z, Y, X) dimensions of your data.  Positions are separated by zarr subgroups under the root zarr store.\n",
-    "\n",
-    "__chunk_size__ determines how zarr will chunk your data.  This means that when you later try to load the data, it will load one chunk at a time with this specified size.  To have the chunk be one z-slice, you would set chunk_size = (1,1,1,Y,X)\n",
-    "\n",
-    "__chan_names__ describe the names of the channels of your data in the order in which they will be written.\n",
-    "\n",
-    "__clims__ corresponds to the the display contrast limits in the metadata for every channel, if none, default values will be used"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 8,
-   "id": "cccfe609",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "data_shape = (3,2,65,256,256)\n",
-    "chunk_size = (1,1,1,256,256)\n",
-    "chan_names = ['DAPI', 'Phase3D']\n",
-    "clims = [(0,1000),(-0.3, 0.3)]\n",
-    "dtype = 'float32'"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 9,
-   "id": "41d7cc52",
-   "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "Creating and opening subgroup Row_0/Col_0/Pos_000\n"
-     ]
-    }
-   ],
-   "source": [
-    "position = 0\n",
-    "writer.init_array(position, data_shape, chunk_size, chan_names, dtype, clims, position_name=None, overwrite=False)"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "id": "435bf54c",
-   "metadata": {},
-   "source": [
-    "### Write the data along specific dimensions"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 10,
-   "id": "201cc9ac",
-   "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "Opening subgroup Row_0/Col_0/Pos_000\n"
-     ]
-    }
-   ],
-   "source": [
-    "writer.write(data[0], p=position, t=slice(0, 3), c=slice(0, 2), z=slice(0, 65))"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "id": "57f06b96",
-   "metadata": {},
-   "source": [
-    "# Complex, 'Position First' Acquisition"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "id": "fc1995e1",
-   "metadata": {},
-   "source": [
-    "Let's say you have a dataset with 10 position with each position containing a data of size (T, C, Z, Y, X) = (3, 2, 65, 256, 256) and you want to save the first time point of every position first, and then move on to the next timepoint.  This will show you how to properly loop through this example case."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 11,
-   "id": "10a6a00a",
-   "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "Creating new zarr store at /Users/cameron.foltz/Desktop/Test/waveOrder_Writer_Example_MultiPos.zarr\n"
-     ]
-    }
-   ],
-   "source": [
-    "writer = WaveorderWriter('/Users/cameron.foltz/Desktop/Test')\n",
-    "writer.create_zarr_root('waveOrder_Writer_Example_MultiPos.zarr')"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 12,
-   "id": "60fc225c",
-   "metadata": {
-    "scrolled": false
-   },
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "Writing Position 0 at Time 0\n",
-      "Writing Position 1 at Time 0\n",
-      "Writing Position 2 at Time 0\n",
-      "Writing Position 3 at Time 0\n",
-      "Writing Position 4 at Time 0\n",
-      "Writing Position 5 at Time 0\n",
-      "Writing Position 6 at Time 0\n",
-      "Writing Position 7 at Time 0\n",
-      "Writing Position 8 at Time 0\n",
-      "Writing Position 9 at Time 0\n",
-      "Writing Position 0 at Time 1\n",
-      "Writing Position 1 at Time 1\n",
-      "Writing Position 2 at Time 1\n",
-      "Writing Position 3 at Time 1\n",
-      "Writing Position 4 at Time 1\n",
-      "Writing Position 5 at Time 1\n",
-      "Writing Position 6 at Time 1\n",
-      "Writing Position 7 at Time 1\n",
-      "Writing Position 8 at Time 1\n",
-      "Writing Position 9 at Time 1\n",
-      "Writing Position 0 at Time 2\n",
-      "Writing Position 1 at Time 2\n",
-      "Writing Position 2 at Time 2\n",
-      "Writing Position 3 at Time 2\n",
-      "Writing Position 4 at Time 2\n",
-      "Writing Position 5 at Time 2\n",
-      "Writing Position 6 at Time 2\n",
-      "Writing Position 7 at Time 2\n",
-      "Writing Position 8 at Time 2\n",
-      "Writing Position 9 at Time 2\n"
-     ]
-    }
-   ],
-   "source": [
-    "# Must first initialize all of the arrays, since you will be returning to them \n",
-    "# at different times in the data-saving process.\n",
-    "num_pos = 10\n",
-    "\n",
-    "for i in range(num_pos):\n",
-    "    writer.init_array(i, data_shape, chunk_size, chan_names, dtype, clims, position_name = None, overwrite=False)    \n",
-    "# Now begin the rime loop \n",
-    "for t in range(3):\n",
-    "    \n",
-    "    # At each time point, write data for all positions\n",
-    "    for pos in range(num_pos):\n",
-    "        \n",
-    "        print(f'Writing Position {pos} at Time {t}')\n",
-    "        \n",
-    "        #write this position data\n",
-    "        writer.write(data[pos, t], p=pos, t=t, c=slice(0, 2), z=slice(0, 65))"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "id": "33b316e6",
-   "metadata": {},
-   "source": [
-    "# Writing Data with HCS Specification"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "id": "914401b6",
-   "metadata": {},
-   "source": [
-    "##### create an HCS 4x4 Plate with rows A, B, C, D and columns 1, 2, 3, 4.  Each well then contains 2 FOV's."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 13,
-   "id": "b0a1b8ef",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "hcs_meta = {'plate':{\n",
-    "  'acquisitions': [{'id': 1,\n",
-    "                     'maximumfieldcount': 2,\n",
-    "                     'name': 'Dataset',\n",
-    "                     'starttime': 0}],\n",
-    "  'columns': [{'name': '1'},\n",
-    "           {'name': '2'},\n",
-    "           {'name': '3'},\n",
-    "           {'name': '4'}],\n",
-    "  'field_count': 2,\n",
-    "  'name': 'MultiWell_Plate_Example',\n",
-    "  'rows': [{'name': 'A'},\n",
-    "          {'name': 'B'},\n",
-    "          {'name': 'C'},\n",
-    "          {'name': 'D'}],\n",
-    "  'version': '0.1',\n",
-    "  'wells': [{'path': 'A/1'},\n",
-    "            {'path': 'A/2'},\n",
-    "            {'path': 'A/3'},\n",
-    "            {'path': 'A/4'},\n",
-    "            {'path': 'B/1'},\n",
-    "            {'path': 'B/2'},\n",
-    "            {'path': 'B/3'},\n",
-    "            {'path': 'B/4'},\n",
-    "            {'path': 'C/1'},\n",
-    "            {'path': 'C/2'},\n",
-    "            {'path': 'C/3'},\n",
-    "            {'path': 'C/4'},\n",
-    "            {'path': 'D/1'},\n",
-    "            {'path': 'D/2'},\n",
-    "            {'path': 'D/3'},\n",
-    "            {'path': 'D/4'}]},\n",
-    "           \n",
-    "       \n",
-    "'well': [{'images': [{'path': 'FOV1'}, {'path': 'FOV2'}]},\n",
-    "         {'images': [{'path': 'FOV1'}, {'path': 'FOV2'}]},\n",
-    "         {'images': [{'path': 'FOV1'}, {'path': 'FOV2'}]},\n",
-    "         {'images': [{'path': 'FOV1'}, {'path': 'FOV2'}]},\n",
-    "         {'images': [{'path': 'FOV1'}, {'path': 'FOV2'}]},\n",
-    "         {'images': [{'path': 'FOV1'}, {'path': 'FOV2'}]},\n",
-    "         {'images': [{'path': 'FOV1'}, {'path': 'FOV2'}]},\n",
-    "         {'images': [{'path': 'FOV1'}, {'path': 'FOV2'}]},\n",
-    "         {'images': [{'path': 'FOV1'}, {'path': 'FOV2'}]},\n",
-    "         {'images': [{'path': 'FOV1'}, {'path': 'FOV2'}]},\n",
-    "         {'images': [{'path': 'FOV1'}, {'path': 'FOV2'}]},\n",
-    "         {'images': [{'path': 'FOV1'}, {'path': 'FOV2'}]},\n",
-    "         {'images': [{'path': 'FOV1'}, {'path': 'FOV2'}]},\n",
-    "         {'images': [{'path': 'FOV1'}, {'path': 'FOV2'}]},\n",
-    "         {'images': [{'path': 'FOV1'}, {'path': 'FOV2'}]},\n",
-    "         {'images': [{'path': 'FOV1'}, {'path': 'FOV2'}]}]\n",
-    "    }\n",
-    "    "
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 14,
-   "id": "56758d8b",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "writer = WaveorderWriter('/Users/cameron.foltz/Desktop/Test', hcs=True, hcs_meta=hcs_meta, verbose=True)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 15,
-   "id": "7e17b39a",
-   "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "Creating new zarr store at /Users/cameron.foltz/Desktop/Test/waveOrder_Writer_ExampleHCS.zarr\n"
-     ]
-    }
-   ],
-   "source": [
-    "writer.create_zarr_root('waveOrder_Writer_ExampleHCS.zarr')"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 16,
-   "id": "5ece924b",
-   "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "/\n",
-      " ├── A\n",
-      " │   ├── 1\n",
-      " │   │   ├── FOV1\n",
-      " │   │   └── FOV2\n",
-      " │   ├── 2\n",
-      " │   │   ├── FOV1\n",
-      " │   │   └── FOV2\n",
-      " │   ├── 3\n",
-      " │   │   ├── FOV1\n",
-      " │   │   └── FOV2\n",
-      " │   └── 4\n",
-      " │       ├── FOV1\n",
-      " │       └── FOV2\n",
-      " ├── B\n",
-      " │   ├── 1\n",
-      " │   │   ├── FOV1\n",
-      " │   │   └── FOV2\n",
-      " │   ├── 2\n",
-      " │   │   ├── FOV1\n",
-      " │   │   └── FOV2\n",
-      " │   ├── 3\n",
-      " │   │   ├── FOV1\n",
-      " │   │   └── FOV2\n",
-      " │   └── 4\n",
-      " │       ├── FOV1\n",
-      " │       └── FOV2\n",
-      " ├── C\n",
-      " │   ├── 1\n",
-      " │   │   ├── FOV1\n",
-      " │   │   └── FOV2\n",
-      " │   ├── 2\n",
-      " │   │   ├── FOV1\n",
-      " │   │   └── FOV2\n",
-      " │   ├── 3\n",
-      " │   │   ├── FOV1\n",
-      " │   │   └── FOV2\n",
-      " │   └── 4\n",
-      " │       ├── FOV1\n",
-      " │       └── FOV2\n",
-      " └── D\n",
-      "     ├── 1\n",
-      "     │   ├── FOV1\n",
-      "     │   └── FOV2\n",
-      "     ├── 2\n",
-      "     │   ├── FOV1\n",
-      "     │   └── FOV2\n",
-      "     ├── 3\n",
-      "     │   ├── FOV1\n",
-      "     │   └── FOV2\n",
-      "     └── 4\n",
-      "         ├── FOV1\n",
-      "         └── FOV2\n"
-     ]
-    }
-   ],
-   "source": [
-    "print(writer.store.tree())"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 17,
-   "id": "1878a5da",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "P = 32\n",
-    "T = 3\n",
-    "C = 2\n",
-    "Z = 11\n",
-    "Y = 256\n",
-    "X = 256\n",
-    "\n",
-    "data = np.random.rand(P, T, C, Z, Y, X)\n",
-    "\n",
-    "data_shape = (3,2,11,256,256)\n",
-    "chunk_size = (1,1,1,256,256)\n",
-    "chan_names = ['DAPI', 'Phase3D']\n",
-    "clims = [(0,1000),(-0.3, 0.3)]\n",
-    "dtype = 'float32'"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 18,
-   "id": "0aed4289",
-   "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "Opening subgroup A/1/FOV1\n",
-      "Opening subgroup A/1/FOV2\n",
-      "Opening subgroup A/2/FOV1\n",
-      "Opening subgroup A/2/FOV2\n",
-      "Opening subgroup A/3/FOV1\n",
-      "Opening subgroup A/3/FOV2\n",
-      "Opening subgroup A/4/FOV1\n",
-      "Opening subgroup A/4/FOV2\n",
-      "Opening subgroup B/1/FOV1\n",
-      "Opening subgroup B/1/FOV2\n",
-      "Opening subgroup B/2/FOV1\n",
-      "Opening subgroup B/2/FOV2\n",
-      "Opening subgroup B/3/FOV1\n",
-      "Opening subgroup B/3/FOV2\n",
-      "Opening subgroup B/4/FOV1\n",
-      "Opening subgroup B/4/FOV2\n",
-      "Opening subgroup C/1/FOV1\n",
-      "Opening subgroup C/1/FOV2\n",
-      "Opening subgroup C/2/FOV1\n",
-      "Opening subgroup C/2/FOV2\n",
-      "Opening subgroup C/3/FOV1\n",
-      "Opening subgroup C/3/FOV2\n",
-      "Opening subgroup C/4/FOV1\n",
-      "Opening subgroup C/4/FOV2\n",
-      "Opening subgroup D/1/FOV1\n",
-      "Opening subgroup D/1/FOV2\n",
-      "Opening subgroup D/2/FOV1\n",
-      "Opening subgroup D/2/FOV2\n",
-      "Opening subgroup D/3/FOV1\n",
-      "Opening subgroup D/3/FOV2\n",
-      "Opening subgroup D/4/FOV1\n",
-      "Opening subgroup D/4/FOV2\n",
-      "Writing Position 0 at Time 0\n",
-      "Opening subgroup A/1/FOV1\n",
-      "Writing Position 1 at Time 0\n",
-      "Opening subgroup A/1/FOV2\n",
-      "Writing Position 2 at Time 0\n",
-      "Opening subgroup A/2/FOV1\n",
-      "Writing Position 3 at Time 0\n",
-      "Opening subgroup A/2/FOV2\n",
-      "Writing Position 4 at Time 0\n",
-      "Opening subgroup A/3/FOV1\n",
-      "Writing Position 5 at Time 0\n",
-      "Opening subgroup A/3/FOV2\n",
-      "Writing Position 6 at Time 0\n",
-      "Opening subgroup A/4/FOV1\n",
-      "Writing Position 7 at Time 0\n",
-      "Opening subgroup A/4/FOV2\n",
-      "Writing Position 8 at Time 0\n",
-      "Opening subgroup B/1/FOV1\n",
-      "Writing Position 9 at Time 0\n",
-      "Opening subgroup B/1/FOV2\n",
-      "Writing Position 10 at Time 0\n",
-      "Opening subgroup B/2/FOV1\n",
-      "Writing Position 11 at Time 0\n",
-      "Opening subgroup B/2/FOV2\n",
-      "Writing Position 12 at Time 0\n",
-      "Opening subgroup B/3/FOV1\n",
-      "Writing Position 13 at Time 0\n",
-      "Opening subgroup B/3/FOV2\n",
-      "Writing Position 14 at Time 0\n",
-      "Opening subgroup B/4/FOV1\n",
-      "Writing Position 15 at Time 0\n",
-      "Opening subgroup B/4/FOV2\n",
-      "Writing Position 16 at Time 0\n",
-      "Opening subgroup C/1/FOV1\n",
-      "Writing Position 17 at Time 0\n",
-      "Opening subgroup C/1/FOV2\n",
-      "Writing Position 18 at Time 0\n",
-      "Opening subgroup C/2/FOV1\n",
-      "Writing Position 19 at Time 0\n",
-      "Opening subgroup C/2/FOV2\n",
-      "Writing Position 20 at Time 0\n",
-      "Opening subgroup C/3/FOV1\n",
-      "Writing Position 21 at Time 0\n",
-      "Opening subgroup C/3/FOV2\n",
-      "Writing Position 22 at Time 0\n",
-      "Opening subgroup C/4/FOV1\n",
-      "Writing Position 23 at Time 0\n",
-      "Opening subgroup C/4/FOV2\n",
-      "Writing Position 24 at Time 0\n",
-      "Opening subgroup D/1/FOV1\n",
-      "Writing Position 25 at Time 0\n",
-      "Opening subgroup D/1/FOV2\n",
-      "Writing Position 26 at Time 0\n",
-      "Opening subgroup D/2/FOV1\n",
-      "Writing Position 27 at Time 0\n",
-      "Opening subgroup D/2/FOV2\n",
-      "Writing Position 28 at Time 0\n",
-      "Opening subgroup D/3/FOV1\n",
-      "Writing Position 29 at Time 0\n",
-      "Opening subgroup D/3/FOV2\n",
-      "Writing Position 30 at Time 0\n",
-      "Opening subgroup D/4/FOV1\n",
-      "Writing Position 31 at Time 0\n",
-      "Opening subgroup D/4/FOV2\n",
-      "Writing Position 0 at Time 1\n",
-      "Opening subgroup A/1/FOV1\n",
-      "Writing Position 1 at Time 1\n",
-      "Opening subgroup A/1/FOV2\n",
-      "Writing Position 2 at Time 1\n",
-      "Opening subgroup A/2/FOV1\n",
-      "Writing Position 3 at Time 1\n",
-      "Opening subgroup A/2/FOV2\n",
-      "Writing Position 4 at Time 1\n",
-      "Opening subgroup A/3/FOV1\n",
-      "Writing Position 5 at Time 1\n",
-      "Opening subgroup A/3/FOV2\n",
-      "Writing Position 6 at Time 1\n",
-      "Opening subgroup A/4/FOV1\n",
-      "Writing Position 7 at Time 1\n",
-      "Opening subgroup A/4/FOV2\n",
-      "Writing Position 8 at Time 1\n",
-      "Opening subgroup B/1/FOV1\n",
-      "Writing Position 9 at Time 1\n",
-      "Opening subgroup B/1/FOV2\n",
-      "Writing Position 10 at Time 1\n",
-      "Opening subgroup B/2/FOV1\n",
-      "Writing Position 11 at Time 1\n",
-      "Opening subgroup B/2/FOV2\n",
-      "Writing Position 12 at Time 1\n",
-      "Opening subgroup B/3/FOV1\n",
-      "Writing Position 13 at Time 1\n",
-      "Opening subgroup B/3/FOV2\n",
-      "Writing Position 14 at Time 1\n",
-      "Opening subgroup B/4/FOV1\n",
-      "Writing Position 15 at Time 1\n",
-      "Opening subgroup B/4/FOV2\n",
-      "Writing Position 16 at Time 1\n",
-      "Opening subgroup C/1/FOV1\n",
-      "Writing Position 17 at Time 1\n",
-      "Opening subgroup C/1/FOV2\n",
-      "Writing Position 18 at Time 1\n",
-      "Opening subgroup C/2/FOV1\n",
-      "Writing Position 19 at Time 1\n",
-      "Opening subgroup C/2/FOV2\n",
-      "Writing Position 20 at Time 1\n",
-      "Opening subgroup C/3/FOV1\n",
-      "Writing Position 21 at Time 1\n",
-      "Opening subgroup C/3/FOV2\n",
-      "Writing Position 22 at Time 1\n",
-      "Opening subgroup C/4/FOV1\n",
-      "Writing Position 23 at Time 1\n",
-      "Opening subgroup C/4/FOV2\n",
-      "Writing Position 24 at Time 1\n",
-      "Opening subgroup D/1/FOV1\n",
-      "Writing Position 25 at Time 1\n",
-      "Opening subgroup D/1/FOV2\n",
-      "Writing Position 26 at Time 1\n",
-      "Opening subgroup D/2/FOV1\n",
-      "Writing Position 27 at Time 1\n",
-      "Opening subgroup D/2/FOV2\n",
-      "Writing Position 28 at Time 1\n",
-      "Opening subgroup D/3/FOV1\n",
-      "Writing Position 29 at Time 1\n",
-      "Opening subgroup D/3/FOV2\n",
-      "Writing Position 30 at Time 1\n",
-      "Opening subgroup D/4/FOV1\n",
-      "Writing Position 31 at Time 1\n",
-      "Opening subgroup D/4/FOV2\n",
-      "Writing Position 0 at Time 2\n",
-      "Opening subgroup A/1/FOV1\n",
-      "Writing Position 1 at Time 2\n",
-      "Opening subgroup A/1/FOV2\n",
-      "Writing Position 2 at Time 2\n",
-      "Opening subgroup A/2/FOV1\n",
-      "Writing Position 3 at Time 2\n",
-      "Opening subgroup A/2/FOV2\n",
-      "Writing Position 4 at Time 2\n",
-      "Opening subgroup A/3/FOV1\n",
-      "Writing Position 5 at Time 2\n",
-      "Opening subgroup A/3/FOV2\n",
-      "Writing Position 6 at Time 2\n",
-      "Opening subgroup A/4/FOV1\n",
-      "Writing Position 7 at Time 2\n",
-      "Opening subgroup A/4/FOV2\n",
-      "Writing Position 8 at Time 2\n",
-      "Opening subgroup B/1/FOV1\n",
-      "Writing Position 9 at Time 2\n",
-      "Opening subgroup B/1/FOV2\n",
-      "Writing Position 10 at Time 2\n",
-      "Opening subgroup B/2/FOV1\n",
-      "Writing Position 11 at Time 2\n",
-      "Opening subgroup B/2/FOV2\n",
-      "Writing Position 12 at Time 2\n",
-      "Opening subgroup B/3/FOV1\n",
-      "Writing Position 13 at Time 2\n",
-      "Opening subgroup B/3/FOV2\n",
-      "Writing Position 14 at Time 2\n",
-      "Opening subgroup B/4/FOV1\n",
-      "Writing Position 15 at Time 2\n",
-      "Opening subgroup B/4/FOV2\n",
-      "Writing Position 16 at Time 2\n",
-      "Opening subgroup C/1/FOV1\n",
-      "Writing Position 17 at Time 2\n",
-      "Opening subgroup C/1/FOV2\n",
-      "Writing Position 18 at Time 2\n",
-      "Opening subgroup C/2/FOV1\n",
-      "Writing Position 19 at Time 2\n",
-      "Opening subgroup C/2/FOV2\n",
-      "Writing Position 20 at Time 2\n",
-      "Opening subgroup C/3/FOV1\n",
-      "Writing Position 21 at Time 2\n",
-      "Opening subgroup C/3/FOV2\n",
-      "Writing Position 22 at Time 2\n",
-      "Opening subgroup C/4/FOV1\n",
-      "Writing Position 23 at Time 2\n",
-      "Opening subgroup C/4/FOV2\n",
-      "Writing Position 24 at Time 2\n",
-      "Opening subgroup D/1/FOV1\n",
-      "Writing Position 25 at Time 2\n",
-      "Opening subgroup D/1/FOV2\n",
-      "Writing Position 26 at Time 2\n",
-      "Opening subgroup D/2/FOV1\n",
-      "Writing Position 27 at Time 2\n",
-      "Opening subgroup D/2/FOV2\n",
-      "Writing Position 28 at Time 2\n",
-      "Opening subgroup D/3/FOV1\n",
-      "Writing Position 29 at Time 2\n",
-      "Opening subgroup D/3/FOV2\n",
-      "Writing Position 30 at Time 2\n",
-      "Opening subgroup D/4/FOV1\n",
-      "Writing Position 31 at Time 2\n",
-      "Opening subgroup D/4/FOV2\n"
-     ]
-    }
-   ],
-   "source": [
-    "# Must first initialize all of the arrays, since you will be returning to them \n",
-    "# at different times in the data-saving process.\n",
-    "num_pos = 32\n",
-    "\n",
-    "for i in range(num_pos):\n",
-    "    writer.init_array(i, data_shape, chunk_size, chan_names, dtype, clims, overwrite=True)    \n",
-    "# Now begin the rime loop \n",
-    "for t in range(3):\n",
-    "    \n",
-    "    # At each time point, write data for all positions\n",
-    "    for pos in range(num_pos):\n",
-    "        \n",
-    "        print(f'Writing Position {pos} at Time {t}')\n",
-    "        \n",
-    "        #write this position data\n",
-    "        writer.write(data[pos, t], p=pos, t=t, c=slice(0, 2), z=slice(0, 65))"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 19,
-   "id": "8b00e1fb",
-   "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "/\n",
-      " ├── A\n",
-      " │   ├── 1\n",
-      " │   │   ├── FOV1\n",
-      " │   │   │   └── array (3, 2, 11, 256, 256) float32\n",
-      " │   │   └── FOV2\n",
-      " │   │       └── array (3, 2, 11, 256, 256) float32\n",
-      " │   ├── 2\n",
-      " │   │   ├── FOV1\n",
-      " │   │   │   └── array (3, 2, 11, 256, 256) float32\n",
-      " │   │   └── FOV2\n",
-      " │   │       └── array (3, 2, 11, 256, 256) float32\n",
-      " │   ├── 3\n",
-      " │   │   ├── FOV1\n",
-      " │   │   │   └── array (3, 2, 11, 256, 256) float32\n",
-      " │   │   └── FOV2\n",
-      " │   │       └── array (3, 2, 11, 256, 256) float32\n",
-      " │   └── 4\n",
-      " │       ├── FOV1\n",
-      " │       │   └── array (3, 2, 11, 256, 256) float32\n",
-      " │       └── FOV2\n",
-      " │           └── array (3, 2, 11, 256, 256) float32\n",
-      " ├── B\n",
-      " │   ├── 1\n",
-      " │   │   ├── FOV1\n",
-      " │   │   │   └── array (3, 2, 11, 256, 256) float32\n",
-      " │   │   └── FOV2\n",
-      " │   │       └── array (3, 2, 11, 256, 256) float32\n",
-      " │   ├── 2\n",
-      " │   │   ├── FOV1\n",
-      " │   │   │   └── array (3, 2, 11, 256, 256) float32\n",
-      " │   │   └── FOV2\n",
-      " │   │       └── array (3, 2, 11, 256, 256) float32\n",
-      " │   ├── 3\n",
-      " │   │   ├── FOV1\n",
-      " │   │   │   └── array (3, 2, 11, 256, 256) float32\n",
-      " │   │   └── FOV2\n",
-      " │   │       └── array (3, 2, 11, 256, 256) float32\n",
-      " │   └── 4\n",
-      " │       ├── FOV1\n",
-      " │       │   └── array (3, 2, 11, 256, 256) float32\n",
-      " │       └── FOV2\n",
-      " │           └── array (3, 2, 11, 256, 256) float32\n",
-      " ├── C\n",
-      " │   ├── 1\n",
-      " │   │   ├── FOV1\n",
-      " │   │   │   └── array (3, 2, 11, 256, 256) float32\n",
-      " │   │   └── FOV2\n",
-      " │   │       └── array (3, 2, 11, 256, 256) float32\n",
-      " │   ├── 2\n",
-      " │   │   ├── FOV1\n",
-      " │   │   │   └── array (3, 2, 11, 256, 256) float32\n",
-      " │   │   └── FOV2\n",
-      " │   │       └── array (3, 2, 11, 256, 256) float32\n",
-      " │   ├── 3\n",
-      " │   │   ├── FOV1\n",
-      " │   │   │   └── array (3, 2, 11, 256, 256) float32\n",
-      " │   │   └── FOV2\n",
-      " │   │       └── array (3, 2, 11, 256, 256) float32\n",
-      " │   └── 4\n",
-      " │       ├── FOV1\n",
-      " │       │   └── array (3, 2, 11, 256, 256) float32\n",
-      " │       └── FOV2\n",
-      " │           └── array (3, 2, 11, 256, 256) float32\n",
-      " └── D\n",
-      "     ├── 1\n",
-      "     │   ├── FOV1\n",
-      "     │   │   └── array (3, 2, 11, 256, 256) float32\n",
-      "     │   └── FOV2\n",
-      "     │       └── array (3, 2, 11, 256, 256) float32\n",
-      "     ├── 2\n",
-      "     │   ├── FOV1\n",
-      "     │   │   └── array (3, 2, 11, 256, 256) float32\n",
-      "     │   └── FOV2\n",
-      "     │       └── array (3, 2, 11, 256, 256) float32\n",
-      "     ├── 3\n",
-      "     │   ├── FOV1\n",
-      "     │   │   └── array (3, 2, 11, 256, 256) float32\n",
-      "     │   └── FOV2\n",
-      "     │       └── array (3, 2, 11, 256, 256) float32\n",
-      "     └── 4\n",
-      "         ├── FOV1\n",
-      "         │   └── array (3, 2, 11, 256, 256) float32\n",
-      "         └── FOV2\n",
-      "             └── array (3, 2, 11, 256, 256) float32\n"
-     ]
-    }
-   ],
-   "source": [
-    "print(writer.store.tree())"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "de34e590",
-   "metadata": {},
-   "outputs": [],
-   "source": []
-  }
- ],
- "metadata": {
-  "kernelspec": {
-   "display_name": "Python 3",
-   "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.7.10"
-  }
- },
- "nbformat": 4,
- "nbformat_minor": 5
-}
\ No newline at end of file
diff --git a/examples/hcs_ome_zarr_writer.py b/examples/hcs_ome_zarr_writer.py
new file mode 100644
index 00000000..39e6867c
--- /dev/null
+++ b/examples/hcs_ome_zarr_writer.py
@@ -0,0 +1,75 @@
+# %%
+# This script shows how to write a high content screening (HCS) OME-Zarr dataset
+# with a single FOV and a single scaling level per well,
+# and how to add an extra well-position to an existing dataset.
+# It can be run as a plain Python script, or as interactive cells in some IDEs.
+
+import zarr
+import numpy as np
+from iohub.zarrfile import HCSReader
+from iohub.writer import new_zarr, HCSWriter
+
+# %%
+# Create a new zarr store and open it
+
+root = new_zarr("hcs.zarr")
+
+# %%
+# Initialize a writer object
+writer = HCSWriter(root, channel_names=["Retardance", "Orientation"])
+
+writer.append_channel("Phase")
+# %%
+# This function writes some random 5-D data to a given position
+
+
+def write_dummy_5d(writer: HCSWriter, position: zarr.Group):
+    tczyx = np.random.randint(
+        0, np.iinfo(np.uint16).max, size=(5, 3, 3, 32, 32), dtype=np.uint16
+    )
+    for t, czyx in enumerate(tczyx):
+        for c, zyx in enumerate(czyx):
+            writer.write_zstack(zyx, position, t, c)
+
+
+# %%
+# Define some positions (row, column, FOV)
+
+positions = (
+    ("A", "1", "0"),
+    ("A", "12", "0"),
+    ("H", "6", "0"),
+)
+
+# %%
+# Write to the positions
+
+for pos in positions:
+    p = writer.require_position(*pos)
+    write_dummy_5d(writer, p)
+
+writer.close()
+
+# %%
+# The sections below shows how to add a new position
+# to an existing store (that we just created).
+
+#%%
+# Populate a new writer from the reader.
+
+reader = HCSReader("hcs.zarr")
+writer = HCSWriter.from_reader(reader)
+
+#%%
+# Alternatively, use the convenience method
+# if default parameters does not need to be changed
+
+while False:  # FIXME: remove this clause and the following indentation
+    writer = HCSWriter.open("hcs.zarr", mode="r+")
+
+
+# %%
+# Add an FOV in a new well
+
+position = writer.require_position("D", "8", "0")
+write_dummy_5d(writer, position)
diff --git a/examples/ome_zarr_writer.py b/examples/ome_zarr_writer.py
new file mode 100644
index 00000000..adf1b1c4
--- /dev/null
+++ b/examples/ome_zarr_writer.py
@@ -0,0 +1,70 @@
+# %%
+# This script shows how to write a single-FOV, single-scale OME-Zarr dataset,
+# and how to append an extra channel to an existing dataset.
+# It can be run as a plain Python script, or as interactive cells in some IDEs.
+
+import numpy as np
+from iohub.writer import new_zarr, OMEZarrWriter
+
+# %%
+# Create a new zarr store and initialize a writer
+
+root = new_zarr("ome.zarr")
+writer = OMEZarrWriter(root, channel_names=["DAPI", "GFP"])
+
+#%%
+# Alternatively, use the convenience method 
+# if default parameters does not need to be changed
+
+while False: # FIXME: remove this clause and the following indentation
+    writer = OMEZarrWriter.open("ome.zarr", mode="a", channel_names=["DAPI", "GFP"])
+
+# %%
+# Generate some random 5-D data
+
+tczyx = np.random.randint(
+    0, np.iinfo(np.uint16).max, size=(5, 2, 3, 32, 32), dtype=np.uint16
+)
+
+# %%
+# Write Z-stacks with metadata
+
+for t, czyx in enumerate(tczyx):
+    for c, zyx in enumerate(czyx):
+        writer.write_zstack(
+            zyx, writer.root, time_index=t, channel_index=c, auto_meta=True
+        )
+
+# %%
+# Done!
+
+writer.close()
+
+# %%
+# The sections below shows how to add a new channel
+# to an existing store (that we just created).
+# It should run without the results from previous steps in RAM.
+
+import zarr
+import numpy as np
+from iohub.writer import OMEZarrWriter
+from iohub.zarrfile import OMEZarrReader
+
+# %%
+# Populate a new writer from the reader.
+
+reader = OMEZarrReader("ome.zarr")
+writer = OMEZarrWriter.from_reader(reader)
+print(writer.channel_names)
+
+# %%
+# Add a new channel
+
+writer.append_channel("Phase3D")
+
+new_tczyx = np.random.randint(
+    0, np.iinfo(np.uint16).max, size=(5, 1, 3, 32, 32), dtype=np.uint16
+)
+
+for t, new_czyx in enumerate(new_tczyx):
+    writer.write_zstack(new_czyx[0], writer.root, t, 2)
diff --git a/iohub/__init__.py b/iohub/__init__.py
index eb453e27..3a1213c8 100644
--- a/iohub/__init__.py
+++ b/iohub/__init__.py
@@ -1,5 +1,3 @@
 #!/usr/bin/env python
 # -*- coding: utf-8 -*-
 
-from .reader import WaveorderReader
-from .writer import WaveorderWriter
diff --git a/iohub/lf_utils.py b/iohub/lf_utils.py
new file mode 100644
index 00000000..0b0ecf61
--- /dev/null
+++ b/iohub/lf_utils.py
@@ -0,0 +1,58 @@
+# Utility functions for label-free microscopy data
+
+import numpy as np
+from typing import Tuple
+
+from iohub.ngff_meta import WindowDict, ChannelMeta
+
+
+def channel_display_settings(
+    chan_name: str,
+    clim: Tuple[float, float, float, float] = None,
+    first_chan: bool = False,
+):
+    """This will create a dictionary used for OME-zarr metadata.  Allows custom contrast limits and channel
+    names for display. Defaults everything to grayscale.
+
+    Parameters
+    ----------
+    chan_name : str
+        Desired name of the channel for display
+    clim : Tuple[float, float, float, float], optional
+        Contrast limits (start, end, min, max)
+    first_chan : bool, optional
+        Whether or not this is the first channel of the dataset (display will be set to active), by default False
+
+    Returns
+    -------
+    dict
+        Display settings adherent to ome-zarr standards
+    """
+    U16_FMAX = float(np.iinfo(np.uint16).max)
+    channel_settings = {
+        "Retardance": (0.0, 100.0, 0.0, 1000.0),
+        "Orientation": (0.0, np.pi, 0.0, np.pi),
+        "Phase2D": (-0.2, 0.2, -10, 10),
+        "Phase3D": (-0.2, 0.2, -10, 10),
+        "BF": (0.0, 5.0, 0.0, U16_FMAX),
+        "S0": (0.0, 1.0, 0.0, U16_FMAX),
+        "S1": (-0.5, 0.5, -10.0, 10.0),
+        "S2": (-0.5, 0.5, -10.0, 10.0),
+        "S3": (-1.0, 1.0, -10.0, -10.0),
+        "Other": (0, U16_FMAX, 0.0, U16_FMAX),
+    }
+    if not clim:
+        if chan_name in channel_settings.keys():
+            clim = channel_settings[chan_name]
+        else:
+            clim = channel_settings["Other"]
+    window = WindowDict(start=clim[0], end=clim[1], min=clim[2], max=clim[3])
+    return ChannelMeta(
+        active=first_chan,
+        coefficient=1.0,
+        color="FFFFFF",
+        family="linear",
+        inverted=False,
+        label=chan_name,
+        window=window,
+    )
diff --git a/iohub/ngff_meta.py b/iohub/ngff_meta.py
new file mode 100644
index 00000000..7194da53
--- /dev/null
+++ b/iohub/ngff_meta.py
@@ -0,0 +1,462 @@
+"""
+Dataclasses with validation for OME-NGFF metadata.
+Developed against OME-NGFF v0.4 and ome-zarr v0.6
+
+Attributes are 'snake_case' with aliases to match NGFF names in JSON output.
+See https://ngff.openmicroscopy.org/0.4/index.html#naming-style about 'camelCase' inconsistency.
+"""
+
+import re, json
+from pydantic import BaseModel, validator, root_validator, Field
+from pydantic.color import Color, ColorTuple, ColorType
+import pandas as pd
+
+from typing import (
+    Optional,
+    Union,
+    Literal,
+    Any,
+    List,
+    Dict,
+    TypedDict,
+    ClassVar,
+)
+
+
+def unique_validator(
+    data: List[Union[BaseModel, TypedDict]], field: Union[str, List[str]]
+):
+    """Called by validators to ensure the uniqueness of certain fields.
+
+    Parameters
+    ----------
+    data : List[Union[BaseModel, TypedDict]]
+        list of pydantic models or typed dictionaries
+    field : Union[str, List[str]]
+        field(s) of the dataclass that must be unique
+
+    Raises
+    ------
+    ValueError
+        raised if any value is not unique
+    """
+    fields = [field] if isinstance(field, str) else field
+    if not isinstance(data[0], dict):
+        data = [d.dict() for d in data]
+    df = pd.DataFrame(data)
+    for key in fields:
+        if not df[key].is_unique:
+            raise ValueError(f"'{key}' must be unique!")
+
+
+def alpha_numeric_validator(data: str):
+    """Called by validators to ensure that strings are alpha-numeric.
+
+    Parameters
+    ----------
+    data : str
+        string to check
+
+    Raises
+    ------
+    ValueError
+        raised if the string contains characters other than [a-zA-z0-9]
+    """
+    if not (data.isalnum() or data.isnumeric()):
+        raise ValueError(
+            f"The path name must be alphanumerical! Got invalid value: '{data}'."
+        )
+
+
+TO_DICT_SETTINGS = dict(exclude_none=True, by_alias=True)
+
+
+class MetaBase(BaseModel):
+    class Config:
+        allow_population_by_field_name = True
+
+
+class AxisMeta(MetaBase):
+    """https://ngff.openmicroscopy.org/0.4/index.html#axes-md"""
+
+    # MUST
+    name: str
+    # SHOULD
+    type: Optional[Literal["space", "time", "channel"]]
+    # SHOULD
+    unit: Optional[str]
+    # store constants as class variables
+    SPACE_UNITS: ClassVar[set] = {
+        "angstrom",
+        "attometer",
+        "centimeter",
+        "decimeter",
+        "exameter",
+        "femtometer",
+        "foot",
+        "gigameter",
+        "hectometer",
+        "inch",
+        "kilometer",
+        "megameter",
+        "meter",
+        "micrometer",
+        "mile",
+        "millimeter",
+        "nanometer",
+        "parsec",
+        "petameter",
+        "picometer",
+        "terameter",
+        "yard",
+        "yoctometer",
+        "yottameter",
+        "zeptometer",
+        "zettameter",
+    }
+    TIME_UNITS: ClassVar[set] = {
+        "attosecond",
+        "centisecond",
+        "day",
+        "decisecond",
+        "exasecond",
+        "femtosecond",
+        "gigasecond",
+        "hectosecond",
+        "hour",
+        "kilosecond",
+        "megasecond",
+        "microsecond",
+        "millisecond",
+        "minute",
+        "nanosecond",
+        "petasecond",
+        "picosecond",
+        "second",
+        "terasecond",
+        "yoctosecond",
+        "yottasecond",
+        "zeptosecond",
+        "zettasecond",
+    }
+
+    @validator("unit")
+    def valid_unit(cls, v, values: dict):
+        if values.get("type") and v is not None:
+            if values["type"] == "channel":
+                raise ValueError(
+                    f"Channel axis must not have a unit! Got unit: {v}."
+                )
+            if values["type"] == "space" and v not in cls.SPACE_UNITS:
+                raise ValueError(
+                    f"Got invalid space unit: '{v}' not in {cls.SPACE_UNITS}"
+                )
+            if values["type"] == "time" and v not in cls.TIME_UNITS:
+                raise ValueError(
+                    f"Got invalid time unit: '{v}' not in {cls.TIME_UNITS}"
+                )
+        return v
+
+
+class TransformationMeta(MetaBase):
+    """https://ngff.openmicroscopy.org/0.4/index.html#trafo-md"""
+
+    # MUST
+    type: Literal["identity", "translation", "scale"]
+    # MUST? (keyword not found in spec for the fields below)
+    translation: Optional[List[float]] = None
+    scale: Optional[List[float]] = None
+    path: Optional[str] = None
+
+    @root_validator
+    def no_extra_method(cls, values: dict):
+        count = sum([bool(v) for _, v in values.items()])
+        if values["type"] == "identity" and count > 1:
+            raise ValueError(
+                f"Method should not be specified for identity transformation!"
+            )
+        elif count > 2:
+            raise ValueError(
+                "Only one type of transformation method is allowed."
+            )
+        return values
+
+
+class DatasetMeta(MetaBase):
+    """https://ngff.openmicroscopy.org/0.4/index.html#multiscale-md"""
+
+    # MUST
+    path: str
+    # MUST
+    coordinate_transformations: List[TransformationMeta] = Field(
+        alias="coordinateTransformations"
+    )
+
+
+class VersionMeta(MetaBase):
+    """OME-NGFF spec version. Default is the current version (0.4)."""
+
+    # SHOULD
+    version: Optional[Literal["0.1", "0.2", "0.3", "0.4"]]
+
+
+class MultiScaleMeta(VersionMeta):
+    """https://ngff.openmicroscopy.org/0.4/index.html#multiscale-md"""
+
+    # MUST
+    axes: List[AxisMeta]
+    # MUST
+    datasets: List[DatasetMeta]
+    # SHOULD
+    name: Optional[str] = None
+    # MAY
+    coordinate_transformations: Optional[List[TransformationMeta]] = Field(
+        alias="coordinateTransformations"
+    )
+    # SHOULD, describes the downscaling method (e.g. 'gaussian')
+    type: Optional[str] = None
+    # SHOULD, additional information about the downscaling method
+    metadata: Optional[dict] = None
+
+    @validator("axes")
+    def unique_name(cls, v):
+        unique_validator(v, "name")
+        return v
+
+    def get_dataset_paths(self):
+        return [dataset.path for dataset in self.datasets]
+
+
+class WindowDict(TypedDict):
+    """Dictionary of contrast limit settings"""
+
+    start: float
+    end: float
+    min: float
+    max: float
+
+
+class ChannelMeta(MetaBase):
+    """Channel display settings without clear documentation from the NGFF spec.
+    https://docs.openmicroscopy.org/omero/5.6.1/developers/Web/WebGateway.html#imgdata"""
+
+    active: bool = False
+    coefficient: float = 1.0
+    color: ColorType = Color("FFFFFF").as_hex()
+    family: str = "linear"
+    inverted: bool = False
+    label: str = None
+    window: WindowDict = None
+
+    class Config:
+        json_encoders = {ColorType: lambda c: Color(c).as_hex()}
+
+
+class RDefsMeta(MetaBase):
+    """Rendering settings without clear documentation from the NGFF spec.
+    https://docs.openmicroscopy.org/omero/5.6.1/developers/Web/WebGateway.html#imgdata"""
+
+    default_t: int = Field(alias="defaultT")
+    default_z: int = Field(alias="defaultZ")
+    model: str = "color"
+    projection: Optional[str] = "normal"
+
+
+class OMEROMeta(VersionMeta):
+    """https://ngff.openmicroscopy.org/0.4/index.html#omero-md"""
+
+    id: int
+    name: Optional[str]
+    channels: Optional[List[ChannelMeta]]
+    rdefs: Optional[RDefsMeta]
+
+
+class ImagesMeta(MetaBase):
+    """Metadata needed for 'Images' (or positions/FOVs) in an OME-NGFF dataset.
+    https://ngff.openmicroscopy.org/0.4/index.html#image-layout"""
+
+    multiscales: List[MultiScaleMeta]
+    omero: OMEROMeta
+
+
+class LabelsMeta(MetaBase):
+    """https://ngff.openmicroscopy.org/0.4/index.html#labels-md"""
+
+    # SHOULD? (keyword not found in spec)
+    labels: str
+    # unlisted groups MAY be labels
+
+
+class LabelColorMeta(MetaBase):
+    """https://ngff.openmicroscopy.org/0.4/index.html#label-md"""
+
+    # MUST
+    label_value: int = Field(alias="label-value")
+    # MAY
+    rgba: Optional[ColorTuple] = None
+
+    class Config:
+        json_encoders = {
+            ColorTuple: lambda c: Color(c).as_rgb_tuple(alpha=True)
+        }
+
+
+class ImageLabelMeta(VersionMeta):
+    """https://ngff.openmicroscopy.org/0.4/index.html#label-md"""
+
+    # SHOULD
+    colors: List[LabelColorMeta]
+    # MAY
+    properties: List[Dict[str, Any]]
+    # MAY
+    source: Dict[str, Any]
+
+    @validator("colors", "properties")
+    def unique_label_value(cls, v):
+        # MUST
+        unique_validator(v, "label_value")
+        return v
+
+
+class AcquisitionMeta(MetaBase):
+    """https://ngff.openmicroscopy.org/0.4/index.html#plate-md"""
+
+    # MUST
+    id: int
+    # SHOULD
+    name: Optional[str] = None
+    # SHOULD
+    maximum_field_count: Optional[int] = Field(alias="maximumfieldcount")
+    # MAY
+    description: Optional[str] = None
+    # MAY
+    start_time: Optional[int] = Field(alias="starttime")
+    # MAY
+    end_time: Optional[int] = Field(alias="endtime")
+
+    @validator("id", "maximum_field_count", "start_time", "end_time")
+    def geq_zero(cls, v):
+        # MUST
+        if v < 0:
+            raise ValueError(
+                "The integer value must be equal or greater to zero!"
+            )
+        return v
+
+    @validator("end_time")
+    def end_after_start(cls, v: int, values: dict):
+        # CUSTOM
+        st = values.get("start_time")
+        if st:
+            if st > v:
+                raise ValueError(
+                    f"The start timestamp {st} should not be larger than the end timestamp {v}"
+                )
+        return v
+
+
+class PlateAxisMeta(MetaBase):
+    """OME-NGFF metadata for a row or a column on a multi-well plate.
+    https://ngff.openmicroscopy.org/0.4/index.html#plate-md"""
+
+    # MUST
+    name: str
+
+    @validator("name")
+    def alpha_numeric(cls, v: str):
+        # MUST
+        alpha_numeric_validator(v)
+        return v
+
+
+class WellIndexMeta(MetaBase):
+    """OME-NGFF metadata for a well on a multi-well plate.
+    https://ngff.openmicroscopy.org/0.4/index.html#plate-md"""
+
+    path: str
+    row_index: int = Field(alias="rowIndex")
+    column_index: int = Field(alias="columnIndex")
+
+    @validator("path")
+    def row_slash_column(cls, v: str):
+        # MUST
+        # regex: one line that is exactly two words separated by one forward slash
+        if len(re.findall(r"^\w+\/\w+$", v)) != 1:
+            raise ValueError(
+                f"The well path '{v}' is not in the form of 'row/column'!"
+            )
+        return v
+
+    @validator("row_index", "column_index")
+    def geq_zero(cls, v: int):
+        # MUST
+        if v < 0:
+            raise ValueError("Well position indices must not be negative!")
+        return v
+
+
+class PlateMeta(VersionMeta):
+    """OME-NGFF high-content screening plate metadata.
+    https://ngff.openmicroscopy.org/0.4/index.html#plate-md"""
+
+    # SHOULD
+    name: Optional[str]
+    # MAY
+    acquisitions: Optional[List[AcquisitionMeta]]
+    # MUST
+    rows: List[PlateAxisMeta]
+    # MUST
+    columns: List[PlateAxisMeta]
+    # MUST
+    wells: List[WellIndexMeta]
+    # SHOULD
+    field_count: Optional[int]
+
+    @validator("acquisitions")
+    def unique_id(cls, v):
+        # MUST
+        unique_validator(v, "id")
+        return v
+
+    @validator("rows", "columns")
+    def unique_name(cls, v):
+        # MUST
+        unique_validator(v, "name")
+        return v
+
+    @validator("wells")
+    def unique_well(cls, v):
+        # CUSTOM
+        unique_validator(v, "path")
+        return v
+
+    @validator("field_count")
+    def positive(cls, v):
+        # MUST
+        if v <= 0:
+            raise ValueError("Field count must be a positive integer!")
+        return v
+
+
+class ImageMeta(MetaBase):
+    """Image metadata field under an HCS well group.
+    https://ngff.openmicroscopy.org/0.4/index.html#well-md"""
+
+    # MUST if `PlateMeta.acquisitions` contains multiple acquisitions
+    acquisition: Optional[int]
+    # MUST
+    path: str
+
+    @validator("path")
+    def alpha_numeric(cls, v):
+        # MUST
+        alpha_numeric_validator(v)
+        return v
+
+
+class WellGroupMeta(VersionMeta):
+    """OME-NGFF high-content screening well group metadata.
+    https://ngff.openmicroscopy.org/0.4/index.html#well-md"""
+
+    # MUST
+    images: List[ImageMeta]
diff --git a/iohub/writer.py b/iohub/writer.py
index a76f8dd6..271200d5 100644
--- a/iohub/writer.py
+++ b/iohub/writer.py
@@ -1,198 +1,709 @@
-import zarr
-import os
-from iohub.writer_structures import HCSZarr, DefaultZarr
+# TODO: remove this in the future (PEP deferred for 3.11, now 3.12?)
+from __future__ import annotations
 
+import os, logging
+from numcodecs import Blosc
+import zarr
+from ome_zarr.format import format_from_version
 
-class WaveorderWriter:
-    """
-    given stokes or physical data, construct a standard hierarchy in zarr for output
-        should conform to the ome-zarr standard as much as possible
+from iohub.zarrfile import OMEZarrReader, HCSReader, _DEFAULT_AXES
+from iohub.ngff_meta import *
+from iohub.lf_utils import channel_display_settings
 
-    """
+from typing import TYPE_CHECKING, Union, Tuple, List, Dict, Literal
+from numpy.typing import NDArray, DTypeLike
 
-    __builder = None
-    __save_dir = None
-    __root_store_path = None
-    __builder_name = None
-    __current_zarr_group = None
-    store = None
+if TYPE_CHECKING:
+    from _typeshed import StrOrBytesPath
 
-    current_group_name = None
-    current_position = None
 
-    def __init__(
-        self,
-        save_dir: str = None,
-        hcs: bool = False,
-        hcs_meta: dict = None,
-        verbose: bool = False,
-    ):
+def new_zarr(
+    store_path: StrOrBytesPath, mode: Literal["r", "r+", "a", "w", "w-"] = "a"
+):
+    """Open the root group of a new OME-NGFF-compatible Zarr store if not present.
 
-        self.verbose = verbose
-        self.use_hcs = hcs
-        self.hcs_meta = hcs_meta
+    Parameters
+    ----------
+    store_path : StrOrBytesPath
+        Path of the new store
+    mode : Literal["r", "r+", "a", "w", "w-"], optional
+        File mode (passed to `zarr.open()`) used to open the root group,
+        by default "a" (read/write, create if not present)
 
-        if os.path.exists(save_dir) and save_dir.endswith(".zarr"):
-            print(f"Opening existing store at {save_dir}")
-            self._open_zarr_root(save_dir)
-        else:
-            self._check_is_dir(save_dir)
-
-        # initialize the subwriter based upon HCS or Default
-        if self.use_hcs:
-            if not self.hcs_meta:
-                raise ValueError(
-                    "No HCS Metadata provided. If HCS format is to be used you must specify the HCS Metadata"
-                )
-
-            self.sub_writer = HCSZarr(
-                self.store, self.__root_store_path, self.hcs_meta
-            )
-        else:
-            self.sub_writer = DefaultZarr(self.store, self.__root_store_path)
+    Returns
+    -------
+    Group
+        Root Zarr group of a directory store with `/` as the dimension separator
+    """
+    if os.path.exists(store_path):
+        raise FileExistsError(f"{store_path} already exists.")
+    store = zarr.DirectoryStore(str(store_path), dimension_separator="/")
+    return zarr.open(store, mode=mode)
+
+
+class OMEZarrWriter:
+    """Generic OME-Zarr writer instance for an existing Zarr store.
+
+    Parameters
+    ----------
+    root : zarr.Group
+        The root group of the Zarr store, dimension separator should be '/'
+    channel_names: List[str]
+        Names of all the channels present in data ordered according to channel indices
+    version : Literal["0.1", "0.4"], optional
+        OME-NGFF version, default to "0.4" if not provided
+    arr_name : str, optional
+        Base name of the arrays, by default '0'
+    axes : List[AxisMeta], optional
+        OME axes metadata, by default:
+        ```
+        [AxisMeta(name='T', type='time', unit='second'),
+        AxisMeta(name='C', type='channel', unit=None),
+        AxisMeta(name='Z', type='space', unit='micrometer'),
+        AxisMeta(name='Y', type='space', unit='micrometer'),
+        AxisMeta(name='X', type='space', unit='micrometer')]
+        ````
+    """
 
-        if self.verbose:
-            self.sub_writer.set_verbosity(self.verbose)
+    _READER_TYPE = OMEZarrReader
 
-    def _check_is_dir(self, path):
-        """
-        directory verification
-        assigns self.__save_dir
+    @classmethod
+    def from_reader(cls, reader: OMEZarrReader):
+        reader.store.close()
+        root = zarr.open(
+            zarr.DirectoryStore(
+                str(reader.store.path), dimension_separator="/"
+            ),
+            mode="a",
+        )
+        writer = cls(
+            root,
+            reader.channel_names,
+            version=reader.version,
+            arr_name=reader.array_keys[0],
+        )
+        writer.images_meta = ImagesMeta(
+            multiscales=root.attrs["multiscales"], omero=root.attrs["omero"]
+        )
+        writer.axes = reader.axes
+        return writer
+
+    @classmethod
+    def open(
+        cls,
+        store_path: StrOrBytesPath,
+        mode: Literal["r+", "a", "w-"] = "a",
+        channel_names: List[str] = None,
+        version: Literal["0.1", "0.4"] = "0.4",
+    ):
+        """Convenience method to open Zarr stores.
+        Uses default parameters to initiate readers/writers. Initiate manually to alter them.
 
         Parameters
         ----------
-        path (str):
+        store_path : StrOrBytesPath
+            Path to the Zarr store to open
+        mode : Literal["r+", "a", "w-"], optional
+            Persistence mode: 'r+' means read/write (must exist); 'a' means read/write (create if doesn't exist); 'w-' means create (fail if exists),
+            by default 'a'
+        channel_names : List[str], optional
+            Channel names used to create a new data store, ignored for existing stores,
+            by default None
+        version : Literal["0.1", "0.4"], optional
+            OME-NGFF version, by default "0.4"
 
         Returns
         -------
-
+        OMEZarrWriter
+            writer instance
         """
-        if os.path.isdir(path) and os.path.exists(path):
-            self.__save_dir = path
-        else:
-            print(
-                f"No existing directory found. Creating new directory at {path}"
+        if mode == "w-" and os.path.exists(store_path):
+            raise FileExistsError(
+                f"Persistence mode 'w-' does not allow overwriting data at {store_path}."
             )
-            os.mkdir(path)
-            self.__save_dir = path
+        try:
+            reader = cls._READER_TYPE(store_path, version=version)
+            logging.info(f"Found existing OME-NGFF dataset at {store_path}")
+            return cls.from_reader(reader)
+        except:
+            not_found_msg = f"OME-NGFF dataset not found at {store_path}"
+            if mode == "r+":
+                raise FileNotFoundError(not_found_msg)
+            else:
+                logging.info(not_found_msg)
+                if not channel_names:
+                    raise ValueError(
+                        "Cannot initiate writer without channel names."
+                    )
+                try:
+                    root = new_zarr(store_path)
+                except FileExistsError:
+                    raise ValueError(
+                        f"Existing data at {store_path} is not a compatible store."
+                    )
+                logging.info(f"Creating new data store at {store_path}")
+                return cls(root, channel_names, version=version)
 
-    def _open_zarr_root(self, path):
-
-        # TODO: Use case where user opens an already HCS-store?
-        """
-        Change current zarr to an existing store
-        if zarr doesn't exist, raise error
+    def __init__(
+        self,
+        root: zarr.Group,
+        channel_names: List[str],
+        version: Literal["0.1", "0.4"] = None,
+        arr_name: str = "0",
+        axes: List[AxisMeta] = None,
+    ):
+        self.root = root
+        self._channel_names = channel_names
+        self.version = version
+        self.fmt = format_from_version(str(version)) if version else "0.4"
+        self.arr_name = arr_name
+        self.axes = axes if axes else _DEFAULT_AXES
+        self._overwrite = False
+
+    @property
+    def channel_names(self):
+        return self._channel_names
+
+    @property
+    def _storage_options(self):
+        return {
+            "compressor": Blosc(
+                cname="zstd", clevel=1, shuffle=Blosc.BITSHUFFLE
+            ),
+            "overwrite": self._overwrite,
+        }
+
+    def _rel_keys(self, rel_path: str) -> List[str]:
+        return [name for name in self.root[rel_path].group_keys()]
+
+    def require_array(
+        self,
+        group: zarr.Group,
+        name: str,
+        zyx_shape: Tuple[int],
+        dtype: DTypeLike,
+        chunks: Tuple[int] = None,
+    ):
+        """Create a new array filled with zeros if it does not exist.
 
         Parameters
         ----------
-        path:       (str) path to store. Must end in .zarr
+        group : zarr.Group
+            Parent zarr group
+        name : str
+            name key of the array
+        zyx_shape : Tuple[int]
+            Shape of the z-stack (Z, Y, X) in the array
+        dtype : DTypeLike
+            Data type
+        chunks : Tuple[int], optional
+            Chunk size for the new array if not present, by default a z-stack (1, 1, Z, Y, X)
 
         Returns
         -------
+        Array
+            Zarr array. Zero-filled with shape (1, 1, Z, Y, X) if created.
 
+        Raises
+        ------
+        ValueError
+        """
+        if len(zyx_shape) != 3:
+            raise ValueError("Shape {shape} does not have 3 dimensions.")
+        value = group.get(name)
+        if value and not isinstance(value, zarr.Array):
+            raise FileExistsError(
+                f"Name '{name}' maps to a non-array object of type {type(value)}."
+            )
+        elif isinstance(value, zarr.Array):
+            if value.shape[-3:] == zyx_shape:
+                return value
+            else:
+                raise FileExistsError(
+                    f"An array exists with incompatible shape {value.shape}."
+                )
+        elif value is None:
+            if not chunks:
+                chunks = (1, 1) + zyx_shape
+            return group.zeros(
+                name,
+                shape=(1, 1, *zyx_shape),
+                chunks=chunks,
+                dtype=dtype,
+                **self._storage_options,
+            )
+
+    def write_zstack(
+        self,
+        data: NDArray,
+        group: zarr.Group,
+        time_index: int,
+        channel_index: int,
+        transform: List[TransformationMeta] = None,
+        name: str = None,
+        auto_meta=True,
+        additional_meta: dict = None,
+    ):
+        """Write a z-stack with OME-NGFF metadata.
+
+        Parameters
+        ----------
+        data : NDArray
+            Image data with shape (Z, Y, X)
+        group : zarr.Group
+            Parent Zarr group
+        time_index : int
+            Time index
+        channel_index : int
+            Channel index
+        transform: List[TransformationMeta], optional
+            Coordinate transformations (see iohub.ngff_meta.TransformationMeta) for the array, by default identity
+        name : str, optional
+            Name key of the 5D array, by default None
+        auto_meta : bool, optional
+            Whether to track and store metadata automatically, by default True
+        additional_meta : dict, optional
+            Additional metadata, by default None
         """
+        if len(data.shape) < 3:
+            raise ValueError("Data has less than 3 dimensions.")
+        if time_index < 0 or channel_index < 0:
+            raise ValueError("Time and channel indices must be non-negative.")
+        if not name:
+            name = self.arr_name
+        zyx_shape = data.shape[-3:]
+        # get 5D array
+        array = self.require_array(
+            group, name, zyx_shape=zyx_shape, dtype=data.dtype
+        )
+        if time_index >= array.shape[0] or channel_index >= array.shape[1]:
+            array.resize(
+                max(time_index + 1, array.shape[0]),
+                max(channel_index + 1, array.shape[1]),
+                *zyx_shape,
+            )
+        # write data
+        array[time_index, channel_index] = data
+        if auto_meta:
+            self._dump_zstack_meta(name, transform, additional_meta)
 
-        if os.path.exists(path):
-            self.store = zarr.open(path)
-            self.__root_store_path = path
+    def _dump_zstack_meta(
+        self,
+        name: str,
+        transform: List[TransformationMeta],
+        additional_meta: dict,
+    ):
+        dataset_meta = self._dataset_meta(name, transform=transform)
+        if "multiscales" not in self.root.attrs:
+            multiscales = [
+                MultiScaleMeta(
+                    version=self.version,
+                    axes=self.axes,
+                    datasets=[dataset_meta],
+                    metadata=additional_meta,
+                )
+            ]
+            omero = self._omero_meta(0, self.root)
+            self.images_meta = ImagesMeta(multiscales=multiscales, omero=omero)
         else:
-            raise FileNotFoundError(
-                f"No store found at {path}, check spelling or create new store with create_zarr"
+            if (
+                dataset_meta.path
+                not in self.images_meta.multiscales[0].get_dataset_paths()
+            ):
+                self.images_meta.multiscales[0].datasets.append(dataset_meta)
+        self.root.attrs.put(self.images_meta.dict(**TO_DICT_SETTINGS))
+
+    def _dataset_meta(
+        self,
+        name: str,
+        transform: List[TransformationMeta] = None,
+    ):
+        if not transform:
+            transform = [TransformationMeta(type="identity")]
+        dataset_meta = DatasetMeta(
+            path=name, coordinate_transformations=transform
+        )
+        return dataset_meta
+
+    def _omero_meta(
+        self,
+        id: int,
+        group: zarr.Group,
+        clims: List[Tuple[float, float, float, float]] = None,
+    ):
+        if not clims:
+            clims = [None] * 4
+        channels = []
+        for i, (channel_name, clim) in enumerate(
+            zip(self.channel_names, clims)
+        ):
+            if i == 0:
+                first_chan = True
+            channels.append(
+                channel_display_settings(
+                    channel_name, clim=clim, first_chan=first_chan
+                )
             )
+        omero_meta = OMEROMeta(
+            version=self.version,
+            id=id,
+            name=group.name,
+            channels=channels,
+            rdefs=RDefsMeta(default_t=0, default_z=0),
+        )
+        return omero_meta
 
-    def create_zarr_root(self, name):
-        """
-        Method for creating the root zarr store.
-        If the store already exists, it will raise an error.
-        Name corresponds to the root directory name (highest level) zarr store.
+    def append_channel(self, chan_name: str):
+        """Append a channel to the end of the channel list.
 
         Parameters
         ----------
-        name:       (string) Name of the zarr store.
-
+        chan_name : str
+            Name of the new channel
         """
+        if chan_name in self._channel_names:
+            raise ValueError(
+                f"Channel name {chan_name} already exists in the dataset."
+            )
+        self._channel_names.append(chan_name)
+        if "omero" in self.root.attrs:
+            self.images_meta.omero.channels.append(
+                channel_display_settings(chan_name)
+            )
+            self.root.attrs["omero"] = self.images_meta.omero.dict(
+                **TO_DICT_SETTINGS
+            )
 
-        if not name.endswith(".zarr"):
-            name = name + ".zarr"
+    def close(self):
+        self.root.store.close()
+
+
+class HCSWriter(OMEZarrWriter):
+    """High-content screening OME-Zarr writer instance for an existing Zarr store.
+
+    Parameters
+    ----------
+    root : zarr.Group
+        The root group of the Zarr store, dimension separator should be '/'
+    channel_names: List[str]
+        Names of all the channels present in data ordered according to channel indices
+    plate_name: str, optional
+        Name of the plate, by default None
+    version : Literal["0.1", "0.4"], optional
+        OME-NGFF version, by default 0.4
+    arr_name : str, optional
+        Base name of the arrays, by default '0'
+    axes : List[AxisMeta], optional
+        OME axes metadata, by default:
+        `
+        [{'name': 'T', 'type': 'time', 'unit': 'second'},
+        {'name': 'C', 'type': 'channel'},
+        {'name': 'Z', 'type': 'space', 'unit': 'micrometer'},
+        {'name': 'Y', 'type': 'space', 'unit': 'micrometer'},
+        {'name': 'X', 'type': 'space', 'unit': 'micrometer'}]
+        `
+    """
 
-        zarr_path = os.path.join(self.__save_dir, name)
-        if os.path.exists(zarr_path):
-            raise FileExistsError("A zarr store with this name already exists")
+    _READER_TYPE = HCSReader
 
-        print(f"Creating new zarr store at {zarr_path}")
-        self.store = zarr.open(zarr_path)
-        self.__root_store_path = zarr_path
-        self.sub_writer.set_store(self.store)
-        self.sub_writer.set_root(self.__root_store_path)
-        self.sub_writer.init_hierarchy()
+    @classmethod
+    def from_reader(
+        cls,
+        reader: HCSReader,
+        detect_arr_name: bool = True,
+        detect_layout: bool = True,
+    ):
+        # TODO: update reader API to use `reader.store` instead of `reader.root.store`
+        reader.root.store.close()
+        root = zarr.open(
+            zarr.DirectoryStore(
+                str(reader.root.store.path), dimension_separator="/"
+            ),
+            mode="a",
+        )
+        writer = cls(
+            root=root,
+            channel_names=reader.channel_names,
+            plate_name=root.attrs.get("name"),
+            version=reader.plate_meta.version,
+            axes=reader.axes,
+            acquisitions=reader.plate_meta.acquisitions,
+        )
+        if detect_arr_name:
+            writer.arr_name = reader.arr_name
+        if detect_layout:
+            writer.rows = reader.rows_meta
+            writer.columns = reader.columns_meta
+            writer.wells = reader.wells_meta
+            writer.positions = reader.positions_meta
+            writer.plate_meta = reader.plate_meta
+        return writer
 
-    def init_array(
+    def __init__(
         self,
-        position,
-        data_shape,
-        chunk_size,
-        chan_names,
-        dtype="float32",
-        clims=None,
-        position_name=None,
-        overwrite=False,
+        root: zarr.Group,
+        channel_names: List[str],
+        plate_name: str = None,
+        version: Literal["0.1", "0.4"] = "0.4",
+        arr_name: str = "0",
+        axes: List[AxisMeta] = None,
+        acquisitions: List[AcquisitionMeta] = None,
     ):
-        """
+        super().__init__(root, channel_names, version, arr_name, axes)
+        self.plate_name = plate_name
+        self.plate_meta: PlateMeta = None
+        self.acquisitions = (
+            [AcquisitionMeta(id=0)] if not acquisitions else acquisitions
+        )
+        self.rows: Dict[str, Union[PlateAxisMeta, int, float]] = {}
+        self.columns: Dict[str, Union[PlateAxisMeta, int, float]] = {}
+        self.wells: Dict[str, Union[WellIndexMeta, int, float]] = {}
+        self.positions: Dict[str, Dict[str, Any]] = {}
 
-        Creates a subgroup structure based on position index.  Then initializes the zarr array under the
-        current position subgroup.  Array level is called 'array' in the hierarchy.
+    @property
+    def row_names(self) -> List[str]:
+        """Row names in the plate (sub-groups under root)"""
+        return self._rel_keys("")
+
+    def get_cols_in_row(self, row_name: str) -> List[str]:
+        """Get column names in a row (wells).
 
         Parameters
         ----------
-        position:           (int) Position index upon which to initialize array
-        data_shape:         (tuple)  Desired Shape of your data (T, C, Z, Y, X).  Must match data
-        chunk_size:         (tuple) Desired Chunk Size (T, C, Z, Y, X).  Chunking each image would be (1, 1, 1, Y, X)
-        dtype:              (str) Data Type, i.e. 'uint16'
-        chan_names:         (list) List of strings corresponding to your channel names.  Used for OME-zarr metadata
-        clims:              (list) list of tuples corresponding to contrast limtis for channel.  OME-Zarr metadata
-        overwrite:          (bool) Whether or not to overwrite the existing data that may be present.
+        row_name : str
+            Name path of the parent row, e.g. `"A"`, `"H"`
 
         Returns
         -------
+        List[str]
+            list of the names of the columns
+        """
+        return self._rel_keys(row_name)
 
+    def get_pos_in_well(self, well_name: str) -> List[str]:
+        """Get column names in a row (wells).
+
+        Parameters
+        ----------
+        well_name : str
+            Name path of the well relative to the root, e.g. `"A/1"`, `"H/12"`
+
+        Returns
+        -------
+        List[str]
+            list of the names of the positions
         """
+        return self._rel_keys(well_name)
 
-        pos_name = position_name if position_name else f"Pos_{position:03d}"
+    @property
+    def plate_layout(self) -> Dict[str, List[str]]:
+        """Names of rows and columns of non-empty wells.
 
-        # Make sure data matches OME zarr structure
-        if len(data_shape) != 5:
-            raise ValueError("Data shape must be (T, C, Z, Y, X)")
+        Returns
+        -------
+        Dict[str, List[str]]
+            {"A": ["1", ...],...}
+        """
+        return {
+            row_name: self.get_cols_in_row() for row_name in self.row_names
+        }
 
-        self.sub_writer.create_position(position, pos_name)
-        self.sub_writer.init_array(
-            data_shape, chunk_size, dtype, chan_names, clims, overwrite
-        )
+    def require_row(
+        self, name: str, index: int = None, overwrite: bool = False
+    ):
+        """Creates a row in the hierarchy (first level below zarr root) if it does not exist.
 
-    def write(self, data, p, t=None, c=None, z=None):
+        Parameters
+        ----------
+        name : str
+            Name of the row
+        index : int, optional
+            Unique index of the new row, by default incremented by 1
+        overwrite : bool, optional
+            Delete all existing data in the row group, by default false
+
+        Returns
+        -------
+        Group
+            Zarr group for the required row
         """
-        Wrapper that calls the builder's write function.
-        Will write to existing array of zeros and place
-        data over the specified indicies.
+        if name in self.rows:
+            if index:
+                existing_id = self.rows[name]["id"]
+                if existing_id != index:
+                    raise ValueError(
+                        f"Requested index {index} conflicts with existing row {name} index {existing_id}."
+                    )
+        else:
+            if not index:
+                index = len(self.rows)
+            self.rows[name] = {
+                "id": index,
+                "meta": PlateAxisMeta(name=name),
+            }
+        return self.root.require_group(name, overwrite=overwrite)
+
+    def require_well(
+        self,
+        row_name: str,
+        col_name: str,
+        row_index: int = None,
+        col_index: int = None,
+        overwrite: bool = False,
+    ):
+        """Creates a well ('row_name/col_name') in the hierarchy if it does not exist.
+        Will also create the parent row if it does not exist
 
         Parameters
         ----------
-        data:   (nd-array), data to be saved. Must be the shape that matches indices (T, C, Z, Y, X)
-        p:      (int), Position index in which to write the data into
-        t:      (slice or int), index or index range of the time dimension
-        c:      (slice or int), index or index range of the channel dimension
-        z:      (slice or int), index or index range of the Z-slice dimension
+        row_name : str
+            Name of the parent row
+        col_name : str
+            Name of the column
+        row_index: int, optional
+            Unique index of the new row, by default incremented by 1
+        col_index: int, optional
+            Unique index of the new column, by default incremented by 1
+        overwrite : bool, optional
+            Delete all existing data in the row group, by default false
+
         Returns
         -------
+        Group
+            Zarr group for the required well
         """
-        self.sub_writer.open_position(p)
+        well_name = os.path.join(row_name, col_name)
+        row = self.require_row(row_name, index=row_index)
+        if col_name in self.columns and col_index:
+            _id = self.columns[col_name]["id"]
+            if _id != col_index:
+                raise ValueError(
+                    f"Requested index {col_index} conflicts with existing column {col_name} index {_id}."
+                )
+        else:
+            if not col_index:
+                col_index = len(self.columns)
+            self.columns[col_name] = {
+                "id": col_index,
+                "meta": PlateAxisMeta(name=col_name),
+            }
+        well = row.require_group(col_name, overwrite=overwrite)
+        if well.name not in self.wells:
+            self.wells[well.name] = {
+                "meta": WellIndexMeta(
+                    path=well_name,
+                    rowIndex=self.rows[row_name]["id"],
+                    columnIndex=self.columns[col_name]["id"],
+                ),
+                "positions": [],
+                "image_meta_list": [],
+            }
+        return well
+
+    def require_position(
+        self, row: str, column: str, fov: str, acq_id: int = 0, **kwargs
+    ):
+        """Create a row, a column, and a FOV/position if they do not exist.
 
-        if t is None:
-            t = slice(0, data.shape[0])
+        Parameters
+        ----------
+        row : str
+            Name key of the row, e.g. 'A'
+        column : str
+            Name key of the column, e.g. '12'
+        fov : str
+            Name key of the FOV/position, e.g. '0'
+        acq_id : int, optional
+            Acquisition ID, by default 0
+        **kwargs :
+            Keyword arguments for `require_well()` and `require_position`
+        """
+        well = self.require_well(row, column, **kwargs)
+        position = well.require_group(fov)
+        self.positions[position.name] = {"id": len(self.positions)}
+        if position.name not in self.wells[well.name]["positions"]:
+            image_meta = ImageMeta(acquisition=acq_id, path=position.basename)
+            self.wells[well.name]["image_meta_list"].append(image_meta)
+            self.wells[well.name]["positions"].append(position.name)
+        return position
+
+    def write_zstack(
+        self,
+        data: NDArray,
+        group: zarr.Group,
+        time_index: int,
+        channel_index: int,
+        transform: List[TransformationMeta] = None,
+        name: str = None,
+        auto_meta=True,
+        additional_meta: dict = None,
+    ):
+        if not name:
+            name = self.arr_name
+        super().write_zstack(
+            data,
+            group,
+            time_index,
+            channel_index,
+            name=name,
+            auto_meta=False,
+        )
+        if auto_meta:
+            self._dump_zstack_meta(group, name, transform, additional_meta)
+        well = self.root.get(os.path.dirname(group.name))
+        self._dump_well_meta(well)
+        self._dump_plate_meta()
 
-        if c is None:
-            c = slice(0, data.shape[1])
+    def _dump_zstack_meta(
+        self,
+        position: zarr.Group,
+        name: str,
+        transform: List[TransformationMeta],
+        additional_meta: dict,
+    ):
+        dataset_meta = self._dataset_meta(name, transform=transform)
+        if "attrs" not in self.positions[position.name]:
+            multiscales = [
+                MultiScaleMeta(
+                    version=self.version,
+                    axes=self.axes,
+                    datasets=[dataset_meta],
+                    metadata=additional_meta,
+                )
+            ]
+            id = self.positions[position.name]["id"]
+            omero = self._omero_meta(id, position)
+            images_meta = ImagesMeta(multiscales=multiscales, omero=omero)
+            self.positions[position.name]["attrs"] = images_meta
+        else:
+            if (
+                dataset_meta.path
+                not in self.positions[position.name]["attrs"]
+                .multiscales[0]
+                .get_dataset_paths()
+            ):
+                self.positions[position.name]["attrs"].multiscales[
+                    0
+                ].datasets.append(dataset_meta)
+        position.attrs.put(
+            self.positions[position.name]["attrs"].dict(**TO_DICT_SETTINGS)
+        )
 
-        if z is None:
-            z = slice(0, data.shape[2])
+    def _dump_plate_meta(self):
+        self.plate_meta = PlateMeta(
+            version=self.version,
+            name=self.plate_name,
+            acquisitions=self.acquisitions,
+            rows=[row["meta"] for _, row in self.rows.items()],
+            columns=[col["meta"] for _, col in self.columns.items()],
+            wells=[well["meta"] for _, well in self.wells.items()],
+            field_count=len(self.positions),
+        )
+        self.root.attrs["plate"] = self.plate_meta.dict(**TO_DICT_SETTINGS)
 
-        self.sub_writer.write(data, t, c, z)
+    def _dump_well_meta(self, well: zarr.Group):
+        well_group_meta = WellGroupMeta(
+            version=self.version,
+            images=self.wells[well.name]["image_meta_list"],
+        )
+        well.attrs["well"] = well_group_meta.dict(**TO_DICT_SETTINGS)
diff --git a/iohub/writer_structures.py b/iohub/writer_structures.py
deleted file mode 100644
index f692799e..00000000
--- a/iohub/writer_structures.py
+++ /dev/null
@@ -1,740 +0,0 @@
-import os
-from numcodecs import Blosc
-import numpy as np
-
-
-class WriterBase:
-    """
-    ABC for all writer types
-    """
-
-    def __init__(self, store, root_path):
-
-        # init common attributes
-        self.store = store
-        self.root_path = root_path
-        self.current_pos_group = None
-        self.current_position = None
-        self.current_well_group = None
-        self.verbose = False
-        self.dtype = None
-
-        # set hardcoded compressor
-        self.__compressor = Blosc(
-            cname="zstd", clevel=1, shuffle=Blosc.BITSHUFFLE
-        )
-
-        # maps to keep track of hierarchies
-        self.rows = dict()
-        self.columns = dict()
-        self.positions = dict()
-
-    # Silence print statements
-    def set_verbosity(self, verbose: bool):
-        self.verbose = verbose
-
-    # Initialize zero array
-    def init_array(
-        self, data_shape, chunk_size, dtype, chan_names, clims, overwrite=False
-    ):
-        """
-
-        Initializes the zarr array under the current position subgroup.
-        array level is called 'arr_0' in the hierarchy.  Sets omero/multiscales metadata based upon
-        chan_names and clims
-
-        Parameters
-        ----------
-        data_shape:         (tuple)  Desired Shape of your data (T, C, Z, Y, X).  Must match data
-        chunk_size:         (tuple) Desired Chunk Size (T, C, Z, Y, X).  Chunking each image would be (1, 1, 1, Y, X)
-        dtype:              (str or np.dtype) Data Type, i.e. 'uint16' or np.uint16
-        chan_names:         (list) List of strings corresponding to your channel names.  Used for OME-zarr metadata
-        clims:              (list) list of tuples corresponding to contrast limtis for channel.  OME-Zarr metadata
-                                    tuple can be of (start, end, min, max) or (start, end)
-        overwrite:          (bool) Whether or not to overwrite the existing data that may be present.
-
-        Returns
-        -------
-
-        """
-
-        self.dtype = np.dtype(dtype)
-
-        self.set_channel_attributes(chan_names, clims)
-        self.current_pos_group.zeros(
-            "arr_0",
-            shape=data_shape,
-            chunks=chunk_size,
-            dtype=dtype,
-            compressor=self.__compressor,
-            overwrite=overwrite,
-        )
-
-    def write(self, data, t, c, z):
-        """
-        Write data to specified index of initialized zarr array
-
-        :param data: (nd-array), data to be saved. Must be the shape that matches indices (T, C, Z, Y, X)
-        :param t: (list), index or index slice of the time dimension
-        :param c: (list), index or index slice of the channel dimension
-        :param z: (list), index or index slice of the z dimension
-
-        """
-
-        shape = np.shape(data)
-
-        if self.current_pos_group.__len__() == 0:
-            raise ValueError("Array not initialized")
-
-        if not isinstance(t, int) and not isinstance(t, slice):
-            raise TypeError("t specification must be either int or slice")
-
-        if not isinstance(c, int) and not isinstance(c, slice):
-            raise TypeError("c specification must be either int or slice")
-
-        if not isinstance(z, int) and not isinstance(z, slice):
-            raise TypeError("z specification must be either int or slice")
-
-        if isinstance(t, int) and isinstance(c, int) and isinstance(z, int):
-
-            if len(shape) > 2:
-                raise ValueError("Index dimensions exceed data dimensions")
-            else:
-                self.current_pos_group["arr_0"][t, c, z] = data
-
-        else:
-            self.current_pos_group["arr_0"][t, c, z] = data
-
-    def create_channel_dict(self, chan_name, clim=None, first_chan=False):
-        """
-        This will create a dictionary used for OME-zarr metadata.  Allows custom contrast limits and channel
-        names for display.  Defaults everything to grayscale.
-
-        Parameters
-        ----------
-        chan_name:          (str) Desired name of the channel for display
-        clim:               (tuple) contrast limits (start, end, min, max)
-        first_chan:         (bool) whether or not this is the first channel of the dataset (display will be set to active)
-
-        Returns
-        -------
-        dict_:              (dict) dictionary adherent to ome-zarr standards
-
-        """
-
-        if chan_name == "Retardance":
-            min = clim[2] if clim else 0.0
-            max = clim[3] if clim else 1000.0
-            start = clim[0] if clim else 0.0
-            end = clim[1] if clim else 100.0
-        elif chan_name == "Orientation":
-            min = clim[2] if clim else 0.0
-            max = clim[3] if clim else np.pi
-            start = clim[0] if clim else 0.0
-            end = clim[1] if clim else np.pi
-
-        elif chan_name == "Phase3D":
-            min = clim[2] if clim else -10.0
-            max = clim[3] if clim else 10.0
-            start = clim[0] if clim else -0.2
-            end = clim[1] if clim else 0.2
-
-        elif chan_name == "BF":
-            min = clim[2] if clim else 0.0
-            max = clim[3] if clim else 65535.0
-            start = clim[0] if clim else 0.0
-            end = clim[1] if clim else 5.0
-
-        elif chan_name == "S0":
-            min = clim[2] if clim else 0.0
-            max = clim[3] if clim else 65535.0
-            start = clim[0] if clim else 0.0
-            end = clim[1] if clim else 1.0
-
-        elif chan_name == "S1":
-            min = clim[2] if clim else 10.0
-            max = clim[3] if clim else -10.0
-            start = clim[0] if clim else -0.5
-            end = clim[1] if clim else 0.5
-
-        elif chan_name == "S2":
-            min = clim[2] if clim else -10.0
-            max = clim[3] if clim else 10.0
-            start = clim[0] if clim else -0.5
-            end = clim[1] if clim else 0.5
-
-        elif chan_name == "S3":
-            min = clim[2] if clim else -10
-            max = clim[3] if clim else 10
-            start = clim[0] if clim else -1.0
-            end = clim[1] if clim else 1.0
-
-        else:
-            min = clim[2] if clim else 0.0
-            max = clim[3] if clim else 65535.0
-            start = clim[0] if clim else 0.0
-            end = clim[1] if clim else 65535.0
-
-        dict_ = {
-            "active": first_chan,
-            "coefficient": 1.0,
-            "color": "FFFFFF",
-            "family": "linear",
-            "inverted": False,
-            "label": chan_name,
-            "window": {"end": end, "max": max, "min": min, "start": start},
-        }
-
-        return dict_
-
-    def create_row(self, idx, name=None):
-        """
-        Creates a row in the hierarchy (first level below zarr store). Option to name
-        this row.  Default is Row_{idx}.  Keeps track of the row name + row index for later
-        metadata creation
-
-        Parameters
-        ----------
-        idx:            (int) Index of the row (order in which it is placed)
-        name:           (str) Optional name to replace default row name
-
-        Returns
-        -------
-
-        """
-
-        row_name = f"Row_{idx}" if not name else name
-        row_path = os.path.join(self.root_path, row_name)
-
-        # check if the user is trying to create a row that already exsits
-        if os.path.exists(row_path):
-            raise FileExistsError(
-                f"A row subgroup with the name {row_name} already exists"
-            )
-        else:
-            self.store.create_group(row_name)
-            self.rows[idx] = row_name
-
-    def create_column(self, row_idx, idx, name=None):
-        """
-        Creates a column in the hierarchy (second level below zarr store, one below row). Option to name
-        this column.  Default is Col_{idx}.  Keeps track of the column name + column index for later
-        metadata creation
-
-        Parameters
-        ----------
-        row_idx:        (int) Index of the row to place the column underneath
-        idx:            (int) Index of the column (order in which it is placed)
-        name:           (str) Optional name to replace default column name
-
-        Returns
-        -------
-
-        """
-
-        col_name = f"Col_{idx}" if not name else name
-        row_name = self.rows[row_idx]
-        col_path = os.path.join(
-            os.path.join(self.root_path, row_name), col_name
-        )
-
-        # check to see if the user is trying to create a row that already exists
-        if os.path.exists(col_path):
-            raise FileExistsError(
-                f"A column subgroup with the name {col_name} already exists"
-            )
-        else:
-            self.store[self.rows[row_idx]].create_group(col_name)
-            self.columns[idx] = col_name
-
-    def open_position(self, position: int):
-        """
-        Opens a position based upon the position index.  It will navigate the rows/column to
-        find where this position is based off of the generation position map which keeps track
-        of this information.  It will set current_pos_group to this position for writing the data
-
-        Parameters
-        ----------
-        position:           (int) Index of the position you wish to open
-
-        Returns
-        -------
-
-        """
-
-        # get row, column, and path to the well
-        row_name = self.positions[position]["row"]
-        col_name = self.positions[position]["col"]
-        well_path = os.path.join(
-            os.path.join(self.root_path, row_name), col_name
-        )
-
-        # check to see if this well exists (row/column)
-        if os.path.exists(well_path):
-            pos_name = self.positions[position]["name"]
-            pos_path = os.path.join(well_path, pos_name)
-
-            # check to see if the position exists
-            if os.path.exists(pos_path):
-
-                if self.verbose:
-                    print(f"Opening subgroup {row_name}/{col_name}/{pos_name}")
-
-                # update trackers to note the current status of the writer
-                self.current_pos_group = self.store[row_name][col_name][
-                    pos_name
-                ]
-                self.current_well_group = self.store[row_name][col_name]
-                self.current_position = position
-
-            else:
-                raise FileNotFoundError(
-                    f"Could not find zarr position subgroup at {row_name}/{col_name}/{pos_name}\
-                                                    Check spelling or create position subgroup with create_position"
-                )
-        else:
-            raise FileNotFoundError(
-                f"Could not find zarr position subgroup at {row_name}/{col_name}/\
-                                                Check spelling or create column/position subgroup with create_position"
-            )
-
-    def set_root(self, root):
-        """
-        set the root path of the zarr store.  Used in the main writer class.
-
-        Parameters
-        ----------
-        root:               (str) path to the zarr store (folder ending in .zarr)
-
-        Returns
-        -------
-
-        """
-        self.root_path = root
-
-    def set_store(self, store):
-        """
-        Sets the zarr store.  Used in the main writer class
-
-        Parameters
-        ----------
-        store:              (Zarr StoreObject) Opened zarr store at the highest level
-
-        Returns
-        -------
-
-        """
-        self.store = store
-
-    def get_zarr(self):
-        return self.current_pos_group
-
-    def set_channel_attributes(self, chan_names, clims=None):
-        """
-        A method for creating ome-zarr metadata dictionary.
-        Channel names are defined by the user, everything else
-        is pre-defined.
-
-        Parameters
-        ----------
-        chan_names:     (list) List of channel names in the order of the channel dimensions
-                                i.e. if 3D Phase is C = 0, list '3DPhase' first.
-
-        clims:          (list of tuples) contrast limits to display for every channel
-
-        """
-
-        rdefs = {
-            "defaultT": 0,
-            "model": "color",
-            "projection": "normal",
-            "defaultZ": 0,
-        }
-
-        multiscale_dict = [{"datasets": [{"path": "arr_0"}], "version": "0.1"}]
-        dict_list = []
-
-        if clims and len(chan_names) < len(clims):
-            raise ValueError(
-                "Contrast Limits specified exceed the number of channels given"
-            )
-
-        for i in range(len(chan_names)):
-            if clims:
-                if len(clims[i]) == 2:
-                    if "float" in self.dtype.name:
-                        clim = (
-                            float(clims[i][0]),
-                            float(clims[i][1]),
-                            -1000,
-                            1000,
-                        )
-                    else:
-                        info = np.iinfo(self.dtype)
-                        clim = (
-                            float(clims[i][0]),
-                            float(clims[i][1]),
-                            info.min,
-                            info.max,
-                        )
-                elif len(clims[i]) == 4:
-                    clim = (
-                        float(clims[i][0]),
-                        float(clims[i][1]),
-                        float(clims[i][2]),
-                        float(clims[i][3]),
-                    )
-                else:
-                    raise ValueError(
-                        "clim specification must a tuple of length 2 or 4"
-                    )
-
-            first_chan = True if i == 0 else False
-            if not clims or i >= len(clims):
-                dict_list.append(
-                    self.create_channel_dict(
-                        chan_names[i], first_chan=first_chan
-                    )
-                )
-            else:
-                dict_list.append(
-                    self.create_channel_dict(
-                        chan_names[i], clim, first_chan=first_chan
-                    )
-                )
-
-        full_dict = {
-            "multiscales": multiscale_dict,
-            "omero": {"channels": dict_list, "rdefs": rdefs, "version": 0.1},
-        }
-
-        self.current_pos_group.attrs.put(full_dict)
-
-    def init_hierarchy(self):
-        pass
-
-    def create_position(self, position: int, name: str):
-        pass
-
-
-class DefaultZarr(WriterBase):
-    """
-    This writer is based off creating a default HCS hierarchy for non-hcs datasets.
-    Currently, we decide that all positions will live under individual columns under
-    a single row.  i.e. this produces the following structure:
-
-    Dataset.zarr
-        ____> Row_0
-            ---> Col_0
-                ---> Pos_000
-            ...
-
-            --> Col_N
-                ---> Pos_N
-
-    We assume this structure in the metadata updating/position creation
-    """
-
-    def __init__(self, store, root_path):
-
-        super().__init__(store, root_path)
-
-        self.dataset_name = None
-        self.plate_meta = dict()
-        self.well_meta = dict()
-
-    def init_hierarchy(self):
-        """
-        method to init the default hierarchy.
-        Will create the first row and initialize metadata fields
-
-        Returns
-        -------
-
-        """
-        self.create_row(0)
-        self.dataset_name = os.path.basename(self.root_path).strip(".zarr")
-
-        self.plate_meta["plate"] = {
-            "acquisitions": [
-                {
-                    "id": 1,
-                    "maximumfieldcount": 1,
-                    "name": "Dataset",
-                    "starttime": 0,
-                }
-            ],
-            "columns": [],
-            "field_count": 1,
-            "name": self.dataset_name,
-            "rows": [],
-            "version": "0.1",
-            "wells": [],
-        }
-
-        self.plate_meta["plate"]["rows"].append({"name": self.rows[0]})
-
-        self.well_meta["well"] = {"images": [], "version": "0.1"}
-        self.well_meta = dict(self.well_meta)
-
-    def create_position(self, position, name):
-        """
-        Creates a column and position subgroup given the index and name.  Name is
-        provided by the main writer class
-
-        Parameters
-        ----------
-        position:           (int) Index of the position to create
-        name:               (str) name of the position subgroup
-
-        Returns
-        -------
-
-        """
-
-        # get row name and create a column
-        row_name = self.rows[0]
-        self.create_column(0, position)
-        col_name = self.columns[position]
-
-        if self.verbose:
-            print(
-                f"Creating and opening subgroup {row_name}/{col_name}/{name}"
-            )
-
-        # create position subgroup
-        self.store[row_name][col_name].create_group(name)
-
-        # update trackers
-        self.current_pos_group = self.store[row_name][col_name][name]
-        self.current_well_group = self.store[row_name][col_name]
-        self.current_position = position
-
-        # update ome-metadata
-        self.positions[position] = {
-            "name": name,
-            "row": row_name,
-            "col": col_name,
-        }
-        self._update_plate_meta(position)
-        self._update_well_meta(position)
-
-    def _update_plate_meta(self, pos):
-        """
-        Updates the plate metadata which lives at the highest level (top store).
-        This metadata carries information on the rows/columns and their paths.
-
-        Parameters
-        ----------
-        pos:            (int) Position index to update the metadata
-
-        Returns
-        -------
-
-        """
-
-        self.plate_meta["plate"]["columns"].append({"name": self.columns[pos]})
-        self.plate_meta["plate"]["wells"].append(
-            {"path": f"{self.rows[0]}/{self.columns[pos]}"}
-        )
-        self.store.attrs.put(self.plate_meta)
-
-    def _update_well_meta(self, pos):
-        """
-        Updates the well metadata which lives at the column level.
-        This metadata carries information about the positions underneath.
-        Assumes only one position will ever be underneath this level.
-
-        Parameters
-        ----------
-        pos:            (int) Index of the position to update
-
-        Returns
-        -------
-
-        """
-
-        self.well_meta["well"]["images"] = [
-            {"path": self.positions[pos]["name"]}
-        ]
-        self.store[self.rows[0]][self.columns[pos]].attrs.put(self.well_meta)
-
-
-class HCSZarr(WriterBase):
-    """
-    This writer version will write new data based upon provided HCS metadata.
-    Useful for when we are reconstructing data that already has a specific HCS heirarchy.
-    allows for the process data structure to match the raw data structure for consistent
-    user browsing.  Also works for copying the structure created by the DefaultZarr writer.
-
-    Note, this doesn't make any new structural changes but strictly copies the one provided
-    """
-
-    def __init__(self, store, root_path, hcs_meta):
-
-        super().__init__(store, root_path)
-        self.hcs_meta = hcs_meta
-
-    def init_hierarchy(self):
-        """
-        Creates the entire Row/Col/Pos hierarchy structure based upon the accepted
-        HCS metadata. Creates everything row by row
-
-        Returns
-        -------
-
-        """
-
-        # check to make sure HCS metadata can be parsed / isn't missing critical info
-        self._check_HCS_meta()
-
-        plate_meta = self.hcs_meta["plate"]
-        row_count = 0
-        col_count = 0
-        well_count = 0
-        pos_count = 0
-
-        # go through rows
-        for row in plate_meta["rows"]:
-
-            self.create_row(row_count, row["name"])
-
-            # go through columns under this row
-            for col in plate_meta["columns"]:
-
-                self.create_column(row_count, col_count, col["name"])
-
-                well_meta = self.hcs_meta["well"][well_count]
-
-                # go through positions under this column
-                for image in well_meta["images"]:
-                    self.positions[pos_count] = {
-                        "name": image["path"],
-                        "row": row["name"],
-                        "col": col["name"],
-                    }
-                    self.store[row["name"]][col["name"]].create_group(
-                        image["path"]
-                    )
-                    pos_count += 1
-
-                # place well_metadata provided from user
-                self.store[row["name"]][col["name"]].attrs.put(
-                    {"well": well_meta}
-                )
-
-                col_count += 1
-                well_count += 1
-
-            col_count = 0
-
-        # place plate metadata as provided by the user
-        self.store.attrs.put({"plate": plate_meta})
-
-    def _check_HCS_meta(self):
-        """
-        checks to make sure the HCS metadata is formatted properly. HCS metadata should be passed in the following
-        structure:
-
-        {'plate':       {<plate_metadata>},
-        'well':         [{well_meta}] # length must be equal to number of wells
-        }
-
-        Only the 'plate' metadata is required here.  All others are optional and will be used if present.
-        Specifications on the specific metadata structures can be found here https://ngff.openmicroscopy.org/0.1/
-
-        Returns
-        -------
-
-        """
-
-        HCS_Defaults = {
-            "acquisitions": [
-                {
-                    "id": 1,
-                    "maxmimumfieldcount": 1,
-                    "name": "Dataset",
-                    "starttime": 0,
-                }
-            ],
-            "field_count": 1,
-            "name": os.path.basename(self.root_path).strip(".zarr"),
-            "version": "0.1",
-        }
-
-        # Check to see if all of the required keys are present
-        if "plate" not in self.hcs_meta.keys():
-            raise KeyError(f"HCS metadata missing plate metadata")
-
-        if "well" not in self.hcs_meta.keys():
-            raise KeyError(f"HCS metadata missing well metadata")
-
-        plate_keys = self.hcs_meta["plate"].keys()
-        for key in plate_keys:
-            if key in HCS_Defaults and key not in plate_keys:
-                self.hcs_meta["plate"][key] = HCS_Defaults[key]
-
-        if "rows" not in plate_keys:
-            raise KeyError("rows key is missing from plate metadata")
-
-        if "columns" not in plate_keys:
-            raise KeyError("columns key is missing from plate metadata")
-
-        # create wells data based on rows, columns if not present
-        if "wells" not in self.hcs_meta["plate"].keys():
-            self.hcs_meta["plate"]["wells"] = [
-                {"path": f"{row}/{col}"}
-                for row in self.hcs_meta["plate"]["rows"]
-                for col in self.hcs_meta["plate"]["columns"]
-            ]
-
-        # Check Plate Meta
-        n_col = len(self.hcs_meta["plate"]["columns"])
-        n_row = len(self.hcs_meta["plate"]["rows"])
-        n_well = len(self.hcs_meta["plate"]["wells"])
-        if n_col * n_row != n_well:
-            raise ValueError(
-                "Plate metadata error: Numbers of rows/columns does not match number of wells"
-            )
-
-        cnt = 0
-        for row in self.hcs_meta["plate"]["rows"]:
-            for col in self.hcs_meta["plate"]["columns"]:
-                row_name = row["name"]
-                col_name = col["name"]
-                well_path = self.hcs_meta["plate"]["wells"][cnt]["path"]
-                if well_path != f"{row_name}/{col_name}":
-                    raise ValueError(
-                        f"Plate metadata error: well path {well_path} \
-                                        does not match row {row_name}, col {col_name}"
-                    )
-                cnt += 1
-
-        # check to make sure well metadata is correct
-        if "well" in self.hcs_meta.keys():
-            for well in self.hcs_meta["well"]:
-                if len(well["images"]) > self.hcs_meta["plate"]["field_count"]:
-                    raise ValueError(
-                        "Well metadata error: number of FOV exceeds maximum field count in plate metadata"
-                    )
-
-    def create_position(self, position, name=None):
-        """
-        Assumes all of the positions were created upon _init_hierarchy called before this function.
-        Will raise an error if the user tries to create a position with this writer version
-
-        Parameters
-        ----------
-        position:           (int) position index to create
-        name:               (str) or None.  Not used.
-
-        Returns
-        -------
-
-        """
-        try:
-            self.open_position(position)
-        except:
-            raise ValueError(
-                "HCS Writer already initialized positions. Cannot create a position that already exists"
-            )
diff --git a/iohub/zarrfile.py b/iohub/zarrfile.py
index 5884e8a8..ea7d16bf 100644
--- a/iohub/zarrfile.py
+++ b/iohub/zarrfile.py
@@ -1,39 +1,147 @@
-import os
+# TODO: remove this in the future (PEP deferred for 3.11, now 3.12?)
+from __future__ import annotations
+
+import os, logging
 import zarr
+import numpy as np
 from copy import copy
+from ome_zarr.io import parse_url
+from ome_zarr.format import format_from_version
+
 from iohub.reader_base import ReaderBase
+from iohub.ngff_meta import *
+
+from typing import TYPE_CHECKING, Literal
+
+if TYPE_CHECKING:
+    from _typeshed import StrOrBytesPath
+
+
+_DEFAULT_AXES = [
+    AxisMeta(name="T", type="time", unit="second"),
+    AxisMeta(name="C", type="channel"),
+    *[
+        AxisMeta(name=i, type="space", unit="micrometer")
+        for i in ("Z", "Y", "X")
+    ],
+]
+
+
+class OMEZarrReader(ReaderBase):
+    """Reader for non-HCS OME-Zarr stores.
+    Zarr array containing a single position/FOV under root is anticipated.
+
+    Parameters
+    ----------
+    store_path : StrOrBytesPath
+        Path to the data store
+    version : Literal["0.1", "0.4"], optional
+        OME-NGFF specification version, by default "0.4"
+
+    Attributes
+    ----------
+    version : Literal["0.1", "0.4"]
+        OME-NGFF specification version
+    store : FSStore
+        Zarr file system store
+    root : Group
+        Root Zarr group holding arrays
+    array_keys : List[str]
+        Name keys of arrays under root (not recursive)
+    channel_names : List[str]
+        Name of the channels
+    axes : List[AxisMeta]
+        Axes metadata
+    """
+
+    def __init__(
+        self,
+        store_path: StrOrBytesPath,
+        version: Literal["0.1", "0.4"] = "0.4",
+    ):
+        super().__init__()
+        # check zarr store
+        if not version == "0.4":
+            logging.warn(
+                "\n".join(
+                    "The OMEZarrReader is only tested against OME-NGFF v0.4.",
+                    f"Requested version {version} may not work properly.",
+                )
+            )
+        location = parse_url(
+            store_path, mode="r", fmt=format_from_version(version)
+        )
+        if not location:
+            raise FileNotFoundError(
+                f"OME-Zarr store not found at {store_path}."
+            )
+        if not location.exists():
+            raise FileNotFoundError(
+                "Array and group metadata not found. Is it an empty store?"
+            )
+        self.version = version
+        self.store = location.store
+        self.root = zarr.open(self.store, mode="r")
+        array_keys = list(self.root.array_keys())
+        if array_keys:
+            self.array_keys = array_keys
+        else:
+            raise FileNotFoundError(
+                "Array not found at top level. Is this an HCS store?"
+            )
+        try:
+            channels: list = self.root.attrs.get("omero").get("channels")
+            self.channel_names = [c["label"] for c in channels]
+        except KeyError:
+            logging.warn(
+                "OMERO channel metadata not found. Channel names cannot be determined."
+            )
+        try:
+            self.axes = [
+                AxisMeta(**ax)
+                for ax in self.root.attrs["multiscales"][0]["axes"]
+            ]
+        except KeyError:
+            logging.warn("Axes meta data not found.")
 
 
 class ZarrReader(ReaderBase):
 
     """
+    .. deprecated:: 0.0.1
+          `ZarrReader` will be removed in future iohub releases, it is replaced by
+          `HCSReader` to enforce upgrade to version 0.4 of the OME-Zarr specification.
+
     Reader for HCS ome-zarr arrays.  OME-zarr structure can be found here: https://ngff.openmicroscopy.org/0.1/
     Also collects the HCS metadata so it can be later copied.
     """
 
-    def __init__(self, zarrfile: str):
+    def __init__(
+        self, store_path: str, version: Literal["0.1", "0.4"] = "0.4"
+    ):
         super().__init__()
 
         # zarr files (.zarr) are directories
-        if not os.path.isdir(zarrfile):
+        if not os.path.isdir(store_path):
             raise ValueError("file does not exist")
-
-        self.zf = zarrfile
-
         try:
-            self.store = zarr.open(self.zf, "r")
+            if version == "0.4":
+                dimension_separator = "/"
+            self.store = zarr.DirectoryStore(
+                store_path, dimension_separator=dimension_separator
+            )
+            self.root = zarr.open(self.store, "r")
         except:
-            raise FileNotFoundError("Supplies path is not a valid zarr store")
-
+            raise FileNotFoundError("Supplies path is not a valid zarr root")
         try:
-            row = self.store[list(self.store.group_keys())[0]]
+            row = self.root[list(self.root.group_keys())[0]]
             col = row[list(row.group_keys())[0]]
             pos = col[list(col.group_keys())[0]]
             self.arr_name = list(pos.array_keys())[0]
         except IndexError:
             raise IndexError("Incompatible zarr format")
 
-        self.plate_meta = self.store.attrs.get("plate")
+        self.plate_meta = self.root.attrs.get("plate")
         self._get_rows()
         self._get_columns()
         self._get_wells()
@@ -41,17 +149,21 @@ def __init__(self, zarrfile: str):
         self._get_positions()
 
         # structure of zarr array
+        first_arr_shape = self.root[self.position_map[0]["well"]][
+            self.position_map[0]["name"]
+        ][self.arr_name].shape
         (
             self.frames,
             self.channels,
             self.slices,
             self.height,
             self.width,
-        ) = self.store[self.position_map[0]["well"]][
-            self.position_map[0]["name"]
-        ][
-            self.arr_name
-        ].shape
+        ) = np.pad(
+            first_arr_shape,
+            (5 - len(first_arr_shape), 0),
+            "constant",
+            constant_values=(1),
+        )
         self.positions = len(self.position_map)
         self.channel_names = []
         self.stage_positions = 0
@@ -122,7 +234,7 @@ def _get_positions(self):
         idx = 0
         # Assumes that the positions are indexed in the order of Row-->Well-->FOV
         for well in self.wells:
-            for pos in self.store[well].attrs.get("well").get("images"):
+            for pos in self.root[well].attrs.get("well").get("images"):
                 name = pos["path"]
                 self.position_map[idx] = {"name": name, "well": well}
                 idx += 1
@@ -141,7 +253,7 @@ def _generate_hcs_meta(self):
 
         well_metas = []
         for well in self.wells:
-            meta = self.store[well].attrs.get("well")
+            meta = self.root[well].attrs.get("well")
             well_metas.append(meta)
 
         self.hcs_meta["well"] = well_metas
@@ -154,7 +266,7 @@ def _set_mm_meta(self):
         -------
 
         """
-        self.mm_meta = self.store.attrs.get("Summary")
+        self.mm_meta = self.root.attrs.get("Summary")
         mm_version = self.mm_meta["MicroManagerVersion"]
 
         if mm_version != "pycromanager":
@@ -191,7 +303,7 @@ def _get_channel_names(self):
         well = self.hcs_meta["plate"]["wells"][0]["path"]
         pos = self.hcs_meta["well"][0]["images"][0]["path"]
 
-        omero_meta = self.store[well][pos].attrs.asdict()["omero"]
+        omero_meta = self.root[well][pos].attrs.asdict()["omero"]
 
         for chan in omero_meta["channels"]:
             self.channel_names.append(chan["label"])
@@ -265,7 +377,7 @@ def get_image_plane_metadata(self, p, c, z):
 
         """
         coord_str = f"({p}, 0, {c}, {z})"
-        return self.store.attrs.get("ImagePlaneMetadata").get(coord_str)
+        return self.root.attrs.get("ImagePlaneMetadata").get(coord_str)
 
     def get_zarr(self, position):
         """
@@ -283,7 +395,7 @@ def get_zarr(self, position):
         pos_info = self.position_map[position]
         well = pos_info["well"]
         pos = pos_info["name"]
-        return self.store[well][pos][self.arr_name]
+        return self.root[well][pos][self.arr_name]
 
     def get_array(self, position):
         """
@@ -323,3 +435,73 @@ def get_image(self, p, t, c, z):
 
     def get_num_positions(self) -> int:
         return self.positions
+
+
+class HCSReader(ZarrReader):
+    def __init__(
+        self,
+        store_path: StrOrBytesPath,
+        version: Literal["0.1", "0.4"] = "0.4",
+    ):
+        super().__init__(store_path, version)
+        self._get_axes_meta()
+        self.plate_meta = PlateMeta(**self.root.attrs["plate"])
+
+    def _get_rows(self):
+        self.rows = []
+        self.rows_meta = {}
+        for i, row in enumerate(self.plate_meta["rows"]):
+            self.rows.append(row["name"])
+            self.rows_meta[row["name"]] = {
+                "id": i,
+                "meta": PlateAxisMeta(**row),
+            }
+
+    def _get_columns(self):
+        self.columns = []
+        self.columns_meta = {}
+        for i, column in enumerate(self.plate_meta["columns"]):
+            self.columns.append(column["name"])
+            self.columns_meta[column["name"]] = {
+                "id": i,
+                "meta": PlateAxisMeta(**column),
+            }
+
+    def _get_wells(self):
+        self.wells = []
+        self.wells_meta = {}
+        self.positions_meta = {}
+        for well in self.plate_meta["wells"]:
+            well_name = well["path"]
+            self.wells.append(well_name)
+            self.wells_meta[well_name] = {
+                "meta": WellIndexMeta(**well),
+                "positions": [],
+                "image_meta_list": [
+                    ImageMeta(**image_meta)
+                    for image_meta in self.root[well_name].attrs["well"][
+                        "images"
+                    ]
+                ],
+            }
+            for _, position in self.root[well_name].groups():
+                self.wells_meta[well_name]["positions"].append(position.name)
+                pos_attrs = ImagesMeta(**position.attrs)
+                self.positions_meta[position.name] = {
+                    "attrs": pos_attrs,
+                    "id": pos_attrs.omero.id,
+                }
+
+    def _get_axes_meta(self):
+        first_position = self.root[next(iter(self.positions_meta))]
+        ms = first_position.attrs.get("multiscales")[0]
+        warning = "Axes metadata not found. Using default."
+        if ms:
+            try:
+                self.axes = MultiScaleMeta(**ms).axes
+                return
+            except KeyError:
+                logging.warn(warning)
+        else:
+            logging.warn(warning)
+        self.axes = _DEFAULT_AXES
diff --git a/setup.cfg b/setup.cfg
index c83d80b5..c8aa350a 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -12,7 +12,6 @@ classifiers =
     Intended Audience :: Science/Research
     License :: OSI Approved :: BSD License
     Programming Language :: Python :: 3 :: Only
-    Programming Language :: Python :: 3.7
     Topic :: Scientific/Engineering
     Topic :: Scientific/Engineering :: Image Processing
     Topic :: Scientific/Engineering :: Visualization
@@ -31,20 +30,17 @@ project_urls =
 [options]
 package = find:
 include_package_data = True
-python_requires = >3.8
+python_requires = >3.9
 setup_requires = setuptools_scm
 # add your package requirements here
 install_requires =
-    numpy>=1.17.4
-    matplotlib>=3.1.1
-    scipy>=1.3.0
-    pywavelets>=1.1.1
-    ipywidgets>=7.5.1
-    opencv-python>=3.4.1
+    pandas>=1.5.2
+    pydantic>=1.10.2
     tifffile>=2021.11.2
-    zarr>=2.6.1
+    ome-zarr>=0.6.1
     natsort>=7.1.1
-    pycromanager>=0.13.2
+    pycromanager>=0.24.0
+    hypothesis>=6.61.0
 
 [options.extras_require]
 dev =


Test patch:
diff --git a/.github/workflows/pytests.yml b/.github/workflows/pytests.yml
index 5713b0d9..8f24f327 100644
--- a/.github/workflows/pytests.yml
+++ b/.github/workflows/pytests.yml
@@ -11,7 +11,7 @@ jobs:
     runs-on: ubuntu-latest
     strategy:
       matrix:
-        python-version: ['3.8', '3.9', '3.10']
+        python-version: ['3.9', '3.10', '3.11']
 
     steps:
     - uses: actions/checkout@v3
diff --git a/tests/conftest.py b/tests/conftest.py
index bf33378e..d4d950d0 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -1,6 +1,7 @@
 import pytest
 import shutil
 import os
+from os.path import join as pjoin
 import random
 from wget import download
 
@@ -24,14 +25,12 @@
     "mm2.0-20201209_4p_20t_5z_3c_512k_1_multipositions",
 }
 
-join = os.path.join
-
 
 @pytest.fixture(scope="session")
 def setup_test_data():
 
-    temp_folder = join(os.getcwd(), "pytest_temp")
-    test_data = os.path.join(temp_folder, "test_data")
+    temp_folder = pjoin(os.getcwd(), "pytest_temp")
+    test_data = pjoin(temp_folder, "test_data")
     if not os.path.isdir(temp_folder):
         os.mkdir(temp_folder)
         print("\nsetting up temp folder")
@@ -39,49 +38,39 @@ def setup_test_data():
         os.mkdir(test_data)
 
     # Zenodo URL
-    url = "https://zenodo.org/record/6983916/files/waveOrder_test_data.zip?download=1"
+    custom_url = (
+        "https://zenodo.org/record/6983916/files/waveOrder_test_data.zip"
+    )
+    # Reference v0.4 HCS dataset from OME
+    # See the last line of https://github.com/ome/ngff/issues/140#issuecomment-1309972511
+    ome_hcs_url = "https://zenodo.org/record/7274533/files/20200812-CardiomyocyteDifferentiation14-Cycle1.zarr.zip"
 
     # download files to temp folder
-    output = join(test_data, "waveOrder_test_data.zip")
     if not os.listdir(test_data):
         print("Downloading test files...")
-        download(url, out=output)
-        shutil.unpack_archive(output, extract_dir=test_data)
+        for url in (custom_url, ome_hcs_url):
+            output = pjoin(test_data, os.path.basename(url))
+            download(url, out=output)
+            shutil.unpack_archive(output, extract_dir=test_data)
 
     yield test_data
 
 
-@pytest.fixture(scope="function")
-def setup_writer_folder():
-    temp_folder = join(os.getcwd(), "pytest_temp", "writer_temp")
-    if not os.path.isdir(temp_folder):
-        os.mkdir(temp_folder)
-        print("\nsetting up temp folder")
-
-    yield temp_folder
-
-    try:
-        # remove temp folder
-        shutil.rmtree(temp_folder)
-    except OSError as e:
-        print(f"Error while deleting temp folder: {e.strerror}")
-
-
 @pytest.fixture(scope="function")
 def setup_mm2gamma_ome_tiffs():
 
-    test_data = os.path.join(
+    test_data = pjoin(
         os.getcwd(), "pytest_temp", "test_data", "MM20_ome-tiffs"
     )
 
     subfolders = [
-        f for f in os.listdir(test_data) if os.path.isdir(join(test_data, f))
+        f for f in os.listdir(test_data) if os.path.isdir(pjoin(test_data, f))
     ]
 
     # specific folder
-    one_folder = join(test_data, subfolders[0])
+    one_folder = pjoin(test_data, subfolders[0])
     # random folder
-    rand_folder = join(test_data, random.choice(subfolders))
+    rand_folder = pjoin(test_data, random.choice(subfolders))
     # return path to unzipped folder containing test images as well as specific folder paths
     yield test_data, one_folder, rand_folder
 
@@ -92,11 +81,11 @@ def setup_mm2gamma_ome_tiffs_incomplete():
     This fixture returns a dataset with 11 timepoints
     """
 
-    test_data = os.path.join(
+    test_data = pjoin(
         os.getcwd(), "pytest_temp", "test_data", "MM20_ometiff_incomplete"
     )
 
-    src = os.path.join(test_data, "mm2.0-20201209_20t_5z_3c_512k_incomplete_1")
+    src = pjoin(test_data, "mm2.0-20201209_20t_5z_3c_512k_incomplete_1")
 
     yield src
 
@@ -104,18 +93,18 @@ def setup_mm2gamma_ome_tiffs_incomplete():
 @pytest.fixture(scope="function")
 def setup_mm2gamma_singlepage_tiffs():
 
-    test_data = os.path.join(
+    test_data = pjoin(
         os.getcwd(), "pytest_temp", "test_data", "MM20_singlepage-tiffs"
     )
 
     subfolders = [
-        f for f in os.listdir(test_data) if os.path.isdir(join(test_data, f))
+        f for f in os.listdir(test_data) if os.path.isdir(pjoin(test_data, f))
     ]
 
     # specific folder
-    one_folder = join(test_data, subfolders[0])
+    one_folder = pjoin(test_data, subfolders[0])
     # random folder
-    rand_folder = join(test_data, random.choice(subfolders))
+    rand_folder = pjoin(test_data, random.choice(subfolders))
     # return path to unzipped folder containing test images as well as specific folder paths
     yield test_data, one_folder, rand_folder
 
@@ -127,13 +116,11 @@ def setup_mm2gamma_singlepage_tiffs_incomplete():
     The MDA definition at start of the experiment specifies 20 timepoints
     """
 
-    test_data = os.path.join(
+    test_data = pjoin(
         os.getcwd(), "pytest_temp", "test_data", "MM20_singlepage_incomplete"
     )
 
-    src = os.path.join(
-        test_data, "mm2.0-20201209_20t_5z_3c_512k_incomplete_1 2"
-    )
+    src = pjoin(test_data, "mm2.0-20201209_20t_5z_3c_512k_incomplete_1 2")
 
     yield src
 
@@ -141,18 +128,18 @@ def setup_mm2gamma_singlepage_tiffs_incomplete():
 @pytest.fixture(scope="function")
 def setup_mm1422_ome_tiffs():
 
-    test_data = os.path.join(
+    test_data = pjoin(
         os.getcwd(), "pytest_temp", "test_data", "MM1422_ome-tiffs"
     )
 
     subfolders = [
-        f for f in os.listdir(test_data) if os.path.isdir(join(test_data, f))
+        f for f in os.listdir(test_data) if os.path.isdir(pjoin(test_data, f))
     ]
 
     # specific folder
-    one_folder = join(test_data, subfolders[0])
+    one_folder = pjoin(test_data, subfolders[0])
     # random folder
-    rand_folder = join(test_data, random.choice(subfolders))
+    rand_folder = pjoin(test_data, random.choice(subfolders))
     # return path to unzipped folder containing test images as well as specific folder paths
     yield test_data, one_folder, rand_folder
 
@@ -160,18 +147,18 @@ def setup_mm1422_ome_tiffs():
 @pytest.fixture(scope="function")
 def setup_mm1422_singlepage_tiffs():
 
-    test_data = os.path.join(
+    test_data = pjoin(
         os.getcwd(), "pytest_temp", "test_data", "MM1422_singlepage-tiffs"
     )
 
     subfolders = [
-        f for f in os.listdir(test_data) if os.path.isdir(join(test_data, f))
+        f for f in os.listdir(test_data) if os.path.isdir(pjoin(test_data, f))
     ]
 
     # specific folder
-    one_folder = join(test_data, subfolders[0])
+    one_folder = pjoin(test_data, subfolders[0])
     # random folder
-    rand_folder = join(test_data, random.choice(subfolders))
+    rand_folder = pjoin(test_data, random.choice(subfolders))
     # return path to unzipped folder containing test images as well as specific folder paths
     yield test_data, one_folder, rand_folder
 
@@ -179,19 +166,27 @@ def setup_mm1422_singlepage_tiffs():
 @pytest.fixture(scope="function")
 def setup_mm2gamma_zarr():
 
-    test_data = os.path.join(
-        os.getcwd(), "pytest_temp", "test_data", "MM20_zarr"
-    )
+    test_data = pjoin(os.getcwd(), "pytest_temp", "test_data", "MM20_zarr")
 
-    zp = os.path.join(test_data, "mm2.0-20201209_4p_2t_5z_1c_512k_1.zarr")
+    zp = pjoin(test_data, "mm2.0-20201209_4p_2t_5z_1c_512k_1.zarr")
 
     # return path to unzipped folder containing test images as well as specific folder paths
     yield zp
 
 
+@pytest.fixture(scope="session")
+def setup_hcs_ref():
+    yield pjoin(
+        os.getcwd(),
+        "pytest_temp",
+        "test_data",
+        "20200812-CardiomyocyteDifferentiation14-Cycle1.zarr",
+    )
+
+
 @pytest.fixture(scope="function")
 def setup_pycromanager_test_data():
-    test_data = os.path.join(
+    test_data = pjoin(
         os.getcwd(), "pytest_temp", "test_data", "MM20_pycromanager"
     )
     datasets = [
@@ -206,7 +201,7 @@ def setup_pycromanager_test_data():
         "mm2.0-20210713_pm0.13.2_7z_1",
     ]
 
-    dataset_dirs = [os.path.join(test_data, ds) for ds in datasets]
+    dataset_dirs = [pjoin(test_data, ds) for ds in datasets]
     first_dir, rand_dir, ptcz_dir = (
         dataset_dirs[0],
         random.choice(dataset_dirs),
diff --git a/tests/writer/test_default_writer.py b/tests/writer/test_default_writer.py
deleted file mode 100644
index 67feefad..00000000
--- a/tests/writer/test_default_writer.py
+++ /dev/null
@@ -1,255 +0,0 @@
-import pytest
-import os
-import zarr
-import numpy as np
-from iohub.writer import WaveorderWriter
-from iohub.writer_structures import DefaultZarr, HCSZarr, WriterBase
-
-
-def test_constructor(setup_writer_folder):
-    """
-    Test that constructor finds correct save directory
-
-    Returns
-    -------
-
-    """
-    folder = setup_writer_folder
-    writer_def = WaveorderWriter(
-        os.path.join(folder, "Test"), hcs=False, hcs_meta=None, verbose=False
-    )
-
-    assert isinstance(writer_def.sub_writer, DefaultZarr)
-    assert isinstance(writer_def.sub_writer, WriterBase)
-
-
-def test_constructor_existing(setup_writer_folder):
-    """
-    Test isntantiating the writer into an existing zarr directory
-
-    Parameters
-    ----------
-    setup_folder
-
-    Returns
-    -------
-
-    """
-
-    folder = setup_writer_folder
-
-    writer = WaveorderWriter(os.path.join(folder, "Test"))
-    writer.create_zarr_root("existing.zarr")
-
-    writer_existing = WaveorderWriter(
-        os.path.join(folder, "Test", "existing.zarr")
-    )
-
-    assert writer_existing.sub_writer.root_path == os.path.join(
-        folder, "Test", "existing.zarr"
-    )
-    assert writer_existing.sub_writer.store is not None
-
-
-def test_create_functions(setup_writer_folder):
-    """
-    Test create root zarr, create position subfolders, and switching between
-    position substores
-
-    Parameters
-    ----------
-    setup_folder
-
-    Returns
-    -------
-
-    """
-
-    folder = setup_writer_folder
-    writer = WaveorderWriter(
-        os.path.join(folder, "Test"), hcs=False, hcs_meta=None, verbose=False
-    )
-
-    writer.create_zarr_root("test_zarr_root")
-
-    assert writer.sub_writer.root_path == os.path.join(
-        folder, "Test", "test_zarr_root.zarr"
-    )
-    assert writer.sub_writer.store is not None
-    assert isinstance(writer.sub_writer.store["Row_0"], zarr.Group)
-
-    # Check Plate Metadata
-    assert "plate" in writer.sub_writer.plate_meta
-    assert "rows" in writer.sub_writer.plate_meta.get("plate").keys()
-    assert "columns" in writer.sub_writer.plate_meta.get("plate").keys()
-    assert "wells" in writer.sub_writer.plate_meta.get("plate").keys()
-    assert len(writer.sub_writer.plate_meta.get("plate").get("wells")) == 0
-    assert len(writer.sub_writer.plate_meta.get("plate").get("columns")) == 0
-    assert len(writer.sub_writer.plate_meta.get("plate").get("rows")) == 1
-
-    # Check Well metadata
-    assert "well" in writer.sub_writer.well_meta
-    assert len(writer.sub_writer.well_meta.get("well").get("images")) == 0
-
-
-def test_init_array(setup_writer_folder):
-    """
-    Test the correct initialization of desired array and the associated
-    metadata
-
-    Parameters
-    ----------
-    setup_folder
-
-    Returns
-    -------
-
-    """
-
-    folder = setup_writer_folder
-    writer = WaveorderWriter(
-        os.path.join(folder, "Test"), hcs=False, hcs_meta=None, verbose=False
-    )
-    writer.create_zarr_root("test_zarr_root")
-
-    data_shape = (3, 3, 21, 128, 128)  # T, C, Z, Y, X
-    chunk_size = (1, 1, 1, 128, 128)
-    chan_names = ["State0", "State1", "State3"]
-    clims = [(-0.5, 0.5), (0, 25), (0, 10000)]
-    dtype = "uint16"
-
-    writer.init_array(
-        0,
-        data_shape,
-        chunk_size,
-        chan_names,
-        dtype,
-        clims,
-        position_name=None,
-        overwrite=False,
-    )
-    writer.init_array(
-        1,
-        data_shape,
-        chunk_size,
-        chan_names,
-        dtype,
-        clims,
-        position_name="Test",
-        overwrite=False,
-    )
-
-    assert isinstance(
-        writer.sub_writer.store["Row_0"]["Col_0"]["Pos_000"], zarr.Group
-    )
-    meta_folder = writer.store["Row_0"]["Col_0"]["Pos_000"]
-    meta = meta_folder.attrs.asdict()
-    array = meta_folder["arr_0"]
-
-    assert meta_folder is not None
-    assert array is not None
-    assert array.shape == data_shape
-    assert array.chunks == chunk_size
-    assert array.dtype == dtype
-
-    assert meta is not None
-    assert "multiscales" in meta
-    assert "omero" in meta
-    assert "rdefs" in meta["omero"]
-
-    print(meta["omero"]["channels"])
-    # Test Chan Names and clims
-    for i in range(len(meta["omero"]["channels"])):
-        assert meta["omero"]["channels"][i]["label"] == chan_names[i]
-        assert meta["omero"]["channels"][i]["window"]["start"] == clims[i][0]
-        assert meta["omero"]["channels"][i]["window"]["end"] == clims[i][1]
-
-    assert isinstance(
-        writer.sub_writer.store["Row_0"]["Col_1"]["Test"], zarr.Group
-    )
-    meta_folder = writer.store["Row_0"]["Col_1"]["Test"]
-    meta = meta_folder.attrs.asdict()
-    array = meta_folder["arr_0"]
-
-    assert meta_folder is not None
-    assert array is not None
-    assert array.shape == data_shape
-    assert array.chunks == chunk_size
-    assert array.dtype == dtype
-
-    assert meta is not None
-    assert "multiscales" in meta
-    assert "omero" in meta
-    assert "rdefs" in meta["omero"]
-
-    # Test Chan Names and clims
-    for i in range(len(meta["omero"]["channels"])):
-        assert meta["omero"]["channels"][i]["label"] == chan_names[i]
-        assert meta["omero"]["channels"][i]["window"]["start"] == clims[i][0]
-        assert meta["omero"]["channels"][i]["window"]["end"] == clims[i][1]
-
-
-def test_write(setup_writer_folder):
-    """
-    Test the write function of the writer
-
-    Parameters
-    ----------
-    setup_folder
-
-    Returns
-    -------
-
-    """
-
-    folder = setup_writer_folder
-    writer = WaveorderWriter(
-        os.path.join(folder, "Test"), hcs=False, hcs_meta=None, verbose=False
-    )
-    writer.create_zarr_root("test_zarr_root")
-
-    data = np.random.randint(
-        1, 60000, size=(3, 3, 11, 128, 128), dtype="uint16"
-    )
-
-    data_shape = data.shape
-    chunk_size = (1, 1, 1, 128, 128)
-    chan_names = ["State0", "State1", "State3"]
-    clims = [(-0.5, 0.5), (0, 25), (0, 10000)]
-    dtype = "uint16"
-
-    writer.init_array(
-        0,
-        data_shape,
-        chunk_size,
-        chan_names,
-        dtype,
-        clims,
-        position_name=None,
-        overwrite=True,
-    )
-
-    # Write single index for each channel
-    writer.write(data[0, 0, 0], p=0, t=0, c=0, z=0)
-    assert np.array_equal(
-        writer.sub_writer.store["Row_0"]["Col_0"]["Pos_000"]["arr_0"][0, 0, 0],
-        data[0, 0, 0],
-    )
-
-    # Write full data
-    writer.write(data, p=0, t=slice(0, 3), c=slice(0, 3), z=slice(0, 11))
-    assert np.array_equal(
-        writer.sub_writer.store["Row_0"]["Col_0"]["Pos_000"]["arr_0"][
-            :, :, :, :, :
-        ],
-        data,
-    )
-
-    # Write full data with alt method
-    writer.write(data, p=0)
-    assert np.array_equal(
-        writer.sub_writer.store["Row_0"]["Col_0"]["Pos_000"]["arr_0"][
-            :, :, :, :, :
-        ],
-        data,
-    )
diff --git a/tests/writer/test_hcs_writer.py b/tests/writer/test_hcs_writer.py
deleted file mode 100644
index 1b31a83c..00000000
--- a/tests/writer/test_hcs_writer.py
+++ /dev/null
@@ -1,287 +0,0 @@
-import pytest
-import os
-import zarr
-import numpy as np
-from iohub.writer import WaveorderWriter
-from iohub.writer_structures import DefaultZarr, HCSZarr, WriterBase
-
-hcs_meta = {
-    "plate": {
-        "acquisitions": [
-            {
-                "id": 1,
-                "maximumfieldcount": 2,
-                "name": "Dataset",
-                "starttime": 0,
-            }
-        ],
-        "columns": [
-            {"name": "1"},
-            {"name": "2"},
-            {"name": "3"},
-            {"name": "4"},
-        ],
-        "field_count": 2,
-        "name": "MultiWell_Plate_Example",
-        "rows": [{"name": "A"}, {"name": "B"}, {"name": "C"}, {"name": "D"}],
-        "version": "0.1",
-        "wells": [
-            {"path": "A/1"},
-            {"path": "A/2"},
-            {"path": "A/3"},
-            {"path": "A/4"},
-            {"path": "B/1"},
-            {"path": "B/2"},
-            {"path": "B/3"},
-            {"path": "B/4"},
-            {"path": "C/1"},
-            {"path": "C/2"},
-            {"path": "C/3"},
-            {"path": "C/4"},
-            {"path": "D/1"},
-            {"path": "D/2"},
-            {"path": "D/3"},
-            {"path": "D/4"},
-        ],
-    },
-    "well": [
-        {"images": [{"path": "FOV1"}, {"path": "FOV2"}]},
-        {"images": [{"path": "FOV1"}, {"path": "FOV2"}]},
-        {"images": [{"path": "FOV1"}, {"path": "FOV2"}]},
-        {"images": [{"path": "FOV1"}, {"path": "FOV2"}]},
-        {"images": [{"path": "FOV1"}, {"path": "FOV2"}]},
-        {"images": [{"path": "FOV1"}, {"path": "FOV2"}]},
-        {"images": [{"path": "FOV1"}, {"path": "FOV2"}]},
-        {"images": [{"path": "FOV1"}, {"path": "FOV2"}]},
-        {"images": [{"path": "FOV1"}, {"path": "FOV2"}]},
-        {"images": [{"path": "FOV1"}, {"path": "FOV2"}]},
-        {"images": [{"path": "FOV1"}, {"path": "FOV2"}]},
-        {"images": [{"path": "FOV1"}, {"path": "FOV2"}]},
-        {"images": [{"path": "FOV1"}, {"path": "FOV2"}]},
-        {"images": [{"path": "FOV1"}, {"path": "FOV2"}]},
-        {"images": [{"path": "FOV1"}, {"path": "FOV2"}]},
-        {"images": [{"path": "FOV1"}, {"path": "FOV2"}]},
-    ],
-}
-
-
-def test_constructor(setup_writer_folder):
-    """
-    Test that constructor finds correct save directory
-
-    Returns
-    -------
-
-    """
-    folder = setup_writer_folder
-    writer_def = WaveorderWriter(
-        os.path.join(folder, "Test"),
-        hcs=True,
-        hcs_meta=hcs_meta,
-        verbose=False,
-    )
-
-    assert isinstance(writer_def.sub_writer, HCSZarr)
-    assert isinstance(writer_def.sub_writer, WriterBase)
-
-
-def test_create_functions(setup_writer_folder):
-    """
-    Test create root zarr, create position subfolders, and switching between
-    position substores
-
-    Parameters
-    ----------
-    setup_folder
-
-    Returns
-    -------
-
-    """
-
-    folder = setup_writer_folder
-    writer = WaveorderWriter(
-        os.path.join(folder, "Test"),
-        hcs=True,
-        hcs_meta=hcs_meta,
-        verbose=False,
-    )
-
-    writer.create_zarr_root("test_zarr_root")
-
-    assert writer.sub_writer.root_path == os.path.join(
-        folder, "Test", "test_zarr_root.zarr"
-    )
-    assert writer.sub_writer.store is not None
-
-    # Check that the correct hierarchy was initialized
-    cnt = 0
-    for row in hcs_meta["plate"]["rows"]:
-        for col in hcs_meta["plate"]["columns"]:
-            for fov in range(1, 3):
-                path = f'{row["name"]}/{col["name"]}/FOV{fov}'
-                assert isinstance(writer.sub_writer.store[path], zarr.Group)
-                assert writer.sub_writer.positions[cnt] == {
-                    "name": f"FOV{fov}",
-                    "row": row["name"],
-                    "col": col["name"],
-                }
-                cnt += 1
-
-
-def test_init_array(setup_writer_folder):
-    """
-    Test the correct initialization of desired array and the associated
-    metadata
-
-    Parameters
-    ----------
-    setup_folder
-
-    Returns
-    -------
-
-    """
-
-    folder = setup_writer_folder
-    writer = WaveorderWriter(
-        os.path.join(folder, "Test"),
-        hcs=True,
-        hcs_meta=hcs_meta,
-        verbose=False,
-    )
-    writer.create_zarr_root("test_zarr_root")
-
-    data_shape = (3, 3, 21, 128, 128)  # T, C, Z, Y, X
-    chunk_size = (1, 1, 1, 128, 128)
-    chan_names = ["State0", "State1", "State3"]
-    clims = [(-0.5, 0.5), (0, 25), (0, 10000)]
-    dtype = "uint16"
-
-    writer.init_array(
-        0,
-        data_shape,
-        chunk_size,
-        chan_names,
-        dtype,
-        clims,
-        position_name=None,
-        overwrite=False,
-    )
-    writer.init_array(
-        11,
-        data_shape,
-        chunk_size,
-        chan_names,
-        dtype,
-        clims,
-        position_name="Test",
-        overwrite=False,
-    )
-
-    assert isinstance(writer.sub_writer.store["A"]["1"]["FOV1"], zarr.Group)
-    meta_folder = writer.store["A"]["1"]["FOV1"]
-    meta = meta_folder.attrs.asdict()
-    array = meta_folder["arr_0"]
-
-    assert meta_folder is not None
-    assert array is not None
-    assert array.shape == data_shape
-    assert array.chunks == chunk_size
-    assert array.dtype == dtype
-
-    assert meta is not None
-    assert "multiscales" in meta
-    assert "omero" in meta
-    assert "rdefs" in meta["omero"]
-
-    # Test Chan Names and clims
-    for i in range(len(meta["omero"]["channels"])):
-        assert meta["omero"]["channels"][i]["label"] == chan_names[i]
-        assert meta["omero"]["channels"][i]["window"]["start"] == clims[i][0]
-        assert meta["omero"]["channels"][i]["window"]["end"] == clims[i][1]
-
-    assert isinstance(writer.sub_writer.store["B"]["2"]["FOV2"], zarr.Group)
-    meta_folder = writer.store["B"]["2"]["FOV2"]
-    meta = meta_folder.attrs.asdict()
-    array = meta_folder["arr_0"]
-
-    assert meta_folder is not None
-    assert array is not None
-    assert array.shape == data_shape
-    assert array.chunks == chunk_size
-    assert array.dtype == dtype
-
-    assert meta is not None
-    assert "multiscales" in meta
-    assert "omero" in meta
-    assert "rdefs" in meta["omero"]
-
-    # Test Chan Names and clims
-    for i in range(len(meta["omero"]["channels"])):
-        assert meta["omero"]["channels"][i]["label"] == chan_names[i]
-        assert meta["omero"]["channels"][i]["window"]["start"] == clims[i][0]
-        assert meta["omero"]["channels"][i]["window"]["end"] == clims[i][1]
-
-
-def test_write(setup_writer_folder):
-    """
-    Test the write function of the writer
-
-    Parameters
-    ----------
-    setup_folder
-
-    Returns
-    -------
-
-    """
-
-    folder = setup_writer_folder
-    writer = WaveorderWriter(
-        os.path.join(folder, "Test"),
-        hcs=True,
-        hcs_meta=hcs_meta,
-        verbose=False,
-    )
-    writer.create_zarr_root("test_zarr_root")
-
-    data = np.random.randint(
-        1, 60000, size=(3, 3, 11, 128, 128), dtype="uint16"
-    )
-
-    data_shape = data.shape
-    chunk_size = (1, 1, 1, 128, 128)
-    chan_names = ["State0", "State1", "State3"]
-    clims = [(-0.5, 0.5), (0, 25), (0, 10000)]
-    dtype = "uint16"
-
-    writer.init_array(
-        0,
-        data_shape,
-        chunk_size,
-        chan_names,
-        dtype,
-        clims,
-        position_name=None,
-        overwrite=True,
-    )
-
-    # Write single index for each channel
-    writer.write(data[0, 0, 0], p=0, t=0, c=0, z=0)
-    assert np.array_equal(
-        writer.sub_writer.store["A"]["1"]["FOV1"]["arr_0"][0, 0, 0],
-        data[0, 0, 0],
-    )
-
-    # Write full data
-    writer.write(data, p=0)
-    assert np.array_equal(
-        writer.sub_writer.store["A"]["1"]["FOV1"]["arr_0"][:, :, :, :, :], data
-    )
-
-    # Write full data with alt method
-    writer.write(data, p=0, t=slice(0, 3), c=slice(0, 3), z=slice(0, 11))
-    assert np.array_equal(
-        writer.sub_writer.store["A"]["1"]["FOV1"]["arr_0"][:, :, :, :, :], data
-    )
diff --git a/tests/writer/test_writer.py b/tests/writer/test_writer.py
new file mode 100644
index 00000000..78feb301
--- /dev/null
+++ b/tests/writer/test_writer.py
@@ -0,0 +1,201 @@
+from __future__ import annotations
+
+import pytest
+import os
+import string
+import shutil
+from contextlib import contextmanager
+from tempfile import TemporaryDirectory
+import zarr
+import numpy as np
+from ome_zarr.io import parse_url
+from ome_zarr.reader import Reader
+from numpy.testing import assert_array_almost_equal
+from hypothesis import given, settings, strategies as st
+import hypothesis.extra.numpy as npst
+from typing import TYPE_CHECKING
+from numpy.typing import NDArray
+
+if TYPE_CHECKING:
+    from _typeshed import StrPath
+
+from iohub.writer import new_zarr, OMEZarrWriter, HCSWriter
+
+
+short_text_st = st.text(min_size=1, max_size=16)
+t_dim_st = st.integers(1, 4)
+c_dim_st = st.integers(1, 4)
+z_dim_st = st.integers(1, 4)
+y_dim_st = st.integers(1, 32)
+x_dim_st = st.integers(1, 32)
+channel_names_st = c_dim_st.flatmap(
+    (
+        lambda c_dim: st.lists(
+            short_text_st, min_size=c_dim, max_size=c_dim, unique=True
+        )
+    )
+)
+short_alpha_numeric = st.text(
+    alphabet=list(
+        string.ascii_lowercase + string.ascii_uppercase + string.digits
+    ),
+    min_size=1,
+    max_size=16,
+)
+row_names_st = st.lists(
+    short_alpha_numeric, min_size=1, max_size=8, unique=True
+)
+col_names_st = st.lists(
+    short_alpha_numeric, min_size=1, max_size=8, unique=True
+)
+
+
+@st.composite
+def _random_5d_with_channels(draw, c_dim: int):
+    arr_shape = (
+        draw(t_dim_st),
+        c_dim,
+        draw(z_dim_st),
+        draw(y_dim_st),
+        draw(x_dim_st),
+    )
+    dtype = draw(
+        st.one_of(
+            npst.integer_dtypes(),
+            npst.unsigned_integer_dtypes(),
+            npst.floating_dtypes(),
+            npst.boolean_dtypes(),
+        )
+    )
+    return draw(npst.arrays(dtype, shape=arr_shape))
+
+
+@st.composite
+def _channels_and_random_5d(draw):
+    channel_names = draw(channel_names_st)
+    random_5d = draw(_random_5d_with_channels(c_dim=len(channel_names)))
+    return channel_names, random_5d
+
+
+def test_new_zarr():
+    """Test `iohub.writer.new_zarr()"""
+    with TemporaryDirectory() as temp_dir:
+        store_path = os.path.join(temp_dir, "new.zarr")
+        root = new_zarr(store_path)
+        assert isinstance(root, zarr.Group)
+        assert isinstance(root.store, zarr.DirectoryStore)
+        assert root.store._dimension_separator == "/"
+        assert root.store.path == store_path
+
+
+def test_new_zarr_same_path():
+    """Test `iohub.writer.new_zarr()"""
+    with TemporaryDirectory() as temp_dir:
+        store_path = os.path.join(temp_dir, "new.zarr")
+        os.mkdir(store_path)
+        with pytest.raises(FileExistsError):
+            _ = new_zarr(store_path)
+
+
+@given(channel_names=channel_names_st, arr_name=short_text_st)
+@settings(max_examples=16)
+def test_init_ome_zarr(channel_names, arr_name):
+    """Test `iohub.writer.OMEZarrWriter.__call__()`"""
+    with TemporaryDirectory() as temp_dir:
+        root = new_zarr(os.path.join(temp_dir, "ome.zarr"))
+        writer = OMEZarrWriter(root, channel_names, arr_name=arr_name)
+        assert writer.channel_names == channel_names
+        assert writer.arr_name == arr_name
+        assert writer._rel_keys("/") == []
+
+
+@contextmanager
+def _temp_ome_zarr_writer(image_5d: NDArray, channel_names: list[str]):
+    try:
+        temp_dir = TemporaryDirectory()
+        writer = OMEZarrWriter.open(
+            os.path.join(temp_dir.name, "ome.zarr"),
+            mode="a",
+            channel_names=channel_names,
+        )
+        for t, time_point in enumerate(image_5d):
+            for c, zstack in enumerate(time_point):
+                writer.write_zstack(zstack, writer.root, t, c)
+        yield writer
+    finally:
+        writer.close()
+        temp_dir.cleanup()
+
+
+@given(channels_and_random_5d=_channels_and_random_5d())
+@settings(max_examples=16, deadline=2000)
+def test_write_ome_zarr(channels_and_random_5d):
+    """Test `iohub.writer.OMEZarrWriter.write_zstack()`"""
+    channel_names, random_5d = channels_and_random_5d
+    with _temp_ome_zarr_writer(random_5d, channel_names) as writer:
+        assert_array_almost_equal(writer.root["0"][:], random_5d)
+        # round-trip test with the offical reader implementation
+        ext_reader = Reader(parse_url(writer.root.store.path))
+        node = list(ext_reader())[0]
+        assert node.metadata["name"] == channel_names
+        assert node.specs[0].datasets == [writer.arr_name]
+        assert node.data[0].shape == random_5d.shape
+        assert node.data[0].dtype == random_5d.dtype
+
+
+@given(channel_names=channel_names_st, arr_name=short_text_st)
+@settings(max_examples=16)
+def test_init_hcs(channel_names, arr_name):
+    with TemporaryDirectory() as temp_dir:
+        root = new_zarr(os.path.join(temp_dir, "hcs.zarr"))
+        writer = HCSWriter(root, channel_names, arr_name=arr_name)
+        assert writer.channel_names == channel_names
+        assert writer.arr_name == arr_name
+        assert writer._rel_keys("/") == []
+
+
+def test_open_hcs_create_empty():
+    """Test `iohub.writer.HCSWriter.open()`"""
+    with TemporaryDirectory() as temp_dir:
+        store_path = os.path.join(temp_dir, "hcs.zarr")
+        writer = HCSWriter.open(store_path, mode="a", channel_names=["GFP"])
+        assert writer.root.store.path == store_path
+        writer.close()
+        with pytest.raises(FileExistsError):
+            _ = HCSWriter.open(store_path, mode="w-", channel_names=["GFP"])
+        with pytest.raises(ValueError):
+            _ = HCSWriter.open(store_path, mode="a", channel_names=["GFP"])
+        with pytest.raises(FileNotFoundError):
+            _ = HCSWriter.open(store_path, mode="r+")
+        with pytest.raises(FileNotFoundError):
+            _ = HCSWriter.open("do-not-exist.zarr", mode="r+")
+
+
+@contextmanager
+def _temp_copy(src: StrPath):
+    """Create a temporary copy of data on disk."""
+    try:
+        temp_dir = TemporaryDirectory()
+        yield shutil.copytree(src, temp_dir.name, dirs_exist_ok=True)
+    finally:
+        temp_dir.cleanup()
+
+
+def test_modify_hcs_ref(setup_hcs_ref):
+    """Test `iohub.writer.HCSWriter.open()`"""
+    with _temp_copy(setup_hcs_ref) as store_path:
+        writer = HCSWriter.open(store_path, mode="r+")
+        assert writer.axes[0].name == "c"
+        assert writer.channel_names == ["DAPI"]
+        writer.append_channel("GFP")
+        assert writer.channel_names == ["DAPI", "GFP"]
+
+
+@given(row_names=row_names_st)
+def test_require_row(row_names: list[str]):
+    with TemporaryDirectory() as temp_dir:
+        store_path = os.path.join(temp_dir, "hcs.zarr")
+        writer = HCSWriter.open(store_path, mode="a", channel_names=["GFP"])
+        for row_name in row_names:
+            writer.require_row(row_name)
+        assert list(writer.rows.keys()) == row_names