diff --git a/tests/conftest.py b/tests/conftest.py
index 987e906..2a5770a 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -27,7 +27,7 @@
 import rasterio as rio
 from click.testing import CliRunner
 from rasterio.transform import from_bounds, from_origin
-from rasterio.warp import transform_bounds
+from rasterio.warp import transform_bounds, array_bounds
 
 from orthority.camera import (
     BrownCamera,
@@ -42,7 +42,10 @@
 
 _dem_resolution = (30.0, 30.0)
 """Default DEM resolution (m)."""
-
+_dem_offset = 825.0
+"""Default DEM average height (m)."""
+_dem_gain = 25.0
+"""Default DEM absolute height variation (m)."""
 
 if '__file__' in globals():
     root_path = Path(__file__).absolute().parents[1]
@@ -50,8 +53,8 @@
     root_path = Path(os.getcwd())
 
 
-def checkerboard(shape, square: int = 25, vals: np.ndarray = None):
-    """Return a checkerboard image given an image shape."""
+def checkerboard(shape: tuple[int, int], square: int = 25, vals: np.ndarray = None) -> np.ndarray:
+    """Return a checkerboard image given an image ``shape``."""
     # adapted from https://stackoverflow.com/questions/2169478/how-to-make-a-checkerboard-in-numpy
     vals = np.array([127, 255], dtype=np.uint8) if vals is None else vals
     coords = np.ogrid[0 : shape[0], 0 : shape[1]]
@@ -59,103 +62,72 @@ def checkerboard(shape, square: int = 25, vals: np.ndarray = None):
     return vals[idx]
 
 
-def sinusoidal(shape: tuple):
-    """Return a sinusoidal surface with z vals 0..1, given an array shape."""
+def sinusoid(shape: tuple[int, int]) -> np.ndarray:
+    """Return a sinusoidal surface with z vals -1..1, given an array ``shape``."""
     # adapted from https://docs.enthought.com/mayavi/mayavi/auto/mlab_helper_functions.html#surf
     x = np.linspace(-4 * np.pi, 4 * np.pi, shape[1])
     y = np.linspace(-4 * np.pi, 4 * np.pi, shape[0]) * shape[0] / shape[1]
     x, y = np.meshgrid(x, y)
 
     array = np.sin(x + y) + np.sin(2 * x - y) + np.cos(3 * x + 4 * y)
-    array -= array.min()
-    array /= array.max()
+    array -= array.mean()
+    array /= np.max(np.abs((array.min(), array.max())))
     return array
 
 
-def ortho_bounds(
-    camera: Camera, dem_min: float = Ortho._egm_minmax[0], include_camera: bool = False
-) -> tuple:
-    """Return ortho bounds for the given ``camera`` at z=``dem_min``."""
-    size = camera._im_size
-    ji = np.array([[0, 0], [0, size[1]], [*size], [size[0], 0]]).T
-    xyz = camera.pixel_to_world_z(ji, dem_min)
-    if include_camera and camera.pos:
+def ortho_bounds(camera: Camera, z: float = Ortho._egm_minmax[0]) -> tuple:
+    """Return (left, bottom, right, top) ortho bounds for the given ``camera`` at z=``z``."""
+    w, h = np.array(camera.im_size) - 1
+    ji = np.array(
+        [[0, 0], [w / 2, 0], [w, 0], [w, h / 2], [w, h], [w / 2, h], [0, h], [0, h / 2]]
+    ).T
+    xyz = camera.pixel_to_world_z(ji, z)
+    if isinstance(camera, FrameCamera) and camera.pos:
         xyz = np.column_stack((xyz, camera.pos))
     return *xyz[:2].min(axis=1), *xyz[:2].max(axis=1)
 
 
-def create_dem(
-    camera: FrameCamera,
-    camera_crs: str,
-    dem_crs: str = None,
-    resolution: tuple = _dem_resolution,
-    dtype: str = 'float32',
-    include_camera=False,
-) -> tuple[np.ndarray, dict]:
+def create_zsurf(
+    bounds: tuple[float],
+    z_off: float = _dem_offset,
+    z_gain: float = _dem_gain,
+    resolution: tuple[float, float] = _dem_resolution,
+) -> tuple[np.ndarray, rio.Affine]:
     """
-    Create a 2 band DEM file that covers the ortho bounds of the given ``camera``.
+    Return a z surface (DEM array) and corresponding transform, such that the surface covers the
+    given ``bounds``, at the given ``resolution``.
 
-    Band 1 is a sinusoidal surface, and band 2, a planar surface.
+    The returned array has 2 bands.  Band 1 is a sinusoidal surface with average height ``z_off``
+    and ~amplitude ``z_gain`` (approx.).  Band 2 is a horizontal plane at height ``z_off``.
     """
-    # TODO: remove dependency on camera.pos?  instead specify mean height and sin amplitude
-    bounds = np.array(ortho_bounds(camera, include_camera=include_camera))
-    size = 1 + np.ceil((bounds[2:] - bounds[:2]) / resolution).astype('int')
-    array = np.stack(
-        (
-            sinusoidal(size[::-1]) * 50 + camera.pos[2] - 200,
-            np.ones(size[::-1]) * (50 / 2) + camera.pos[2] - 200,
-        ),
-        axis=0,
-    ).astype(dtype)
-
-    if dem_crs is not None:
-        camera_crs = rio.CRS.from_string(camera_crs)
-        dem_crs = rio.CRS.from_string(dem_crs)
-        bounds = transform_bounds(camera_crs, dem_crs, *bounds)
-        transform = from_bounds(*bounds, *size)
-    else:
-        dem_crs = camera_crs
-        transform = from_origin(bounds[0], bounds[3], *resolution)
-
-    profile = dict(
-        crs=dem_crs,
-        transform=transform,
-        dtype=dtype,
-        width=size[0],
-        height=size[1],
-        count=2,
-        nodata=float('nan'),
-    )
-    return array, profile
-
-
-def create_src(
-    filename: Path,
-    size: tuple,
-    dtype: str = 'uint8',
-    count: int = 3,
-    camera: FrameCamera = None,
-    crs: str = None,
-):
-    """Create a source checkerboard file with optional CRS where ``camera`` & ``crs`` are
-    specified.
-    """
-    profile = dict(
-        crs=None, transform=None, dtype=dtype, width=size[0], height=size[1], count=count
-    )
+    bounds = np.array(bounds)
+    shape = np.ceil((bounds[2:] - bounds[:2]) / resolution).astype('int')[::-1]
 
-    if camera is not None:
-        # find bounds & transform
-        size = camera._im_size
-        bounds = ortho_bounds(camera, dem_min=camera.pos[2] - 100)
-        transform = from_bounds(*bounds, *size)
-        profile.update(crs=crs, transform=transform)
+    array = np.stack((sinusoid(shape) * z_gain + z_off, np.ones(shape) * z_off), axis=0)
+    array = array.astype('float32')
+    transform = from_origin(bounds[0], bounds[3], *resolution)
+    return array, transform
 
-    array = checkerboard(size[::-1])
-    array = np.stack((array,) * count, axis=0).astype(dtype)
 
-    with rio.open(filename, 'w', **profile) as src_im:
-        src_im.write(array)
+def create_profile(
+    array: np.ndarray,
+    transform: rio.Affine = None,
+    crs: str | rio.CRS = None,
+    nodata: int | float = None,
+) -> dict:
+    """Return a Rasterio profile for the given parameters."""
+    if array.ndim != 2 and array.ndim != 3:
+        raise ValueError("'array' should be 2D or 3D.")
+    shape = (1, *array.shape) if array.ndim == 2 else array.shape
+    return dict(
+        crs=crs,
+        transform=transform,
+        dtype=array.dtype,
+        width=shape[2],
+        height=shape[1],
+        count=shape[0],
+        nodata=nodata,
+    )
 
 
 def oty_to_osfm_int_param(int_param_dict: dict) -> dict:
@@ -281,6 +253,8 @@ def pinhole_camera(frame_args: dict) -> FrameCamera:
 @pytest.fixture(scope='session')
 def pinhole_camera_und(frame_args: dict) -> FrameCamera:
     """Example ``PinholeCamera`` object with near-nadir orientation and ``distort=False``."""
+    # TODO: use parameterised fixtures for these and similar?
+    #  https://docs.pytest.org/en/latest/example/parametrize.html#indirect-parametrization
     return PinholeCamera(**frame_args, distort=False)
 
 
@@ -320,6 +294,20 @@ def fisheye_camera_und(frame_args: dict, fisheye_dist_param: dict) -> FrameCamer
     return FisheyeCamera(**frame_args, **fisheye_dist_param, distort=False)
 
 
+@pytest.fixture(scope='session')
+def xyz_grids(
+    pinhole_camera: Camera,
+) -> tuple[tuple[np.ndarray, np.ndarray, np.ndarray], rio.Affine]:
+    """(x, y, z) coordinate grids and transform at 5m resolution.  Z grid has 2 bands: band 1 is
+    a sinusoidal surface, and band 2, a horizontal plane.
+    """
+    bounds = ortho_bounds(pinhole_camera)
+    z, transform = create_zsurf(bounds, resolution=(5, 5))
+    j, i = np.meshgrid(range(0, z.shape[2]), range(0, z.shape[1]), indexing='xy')
+    x, y = (transform * rio.Affine.translation(0.5, 0.5)) * [j, i]
+    return (x, y, z), transform
+
+
 @pytest.fixture(scope='session')
 def utm34n_crs() -> str:
     """CRS string for UTM zone 34N with no vertical CRS."""
@@ -401,51 +389,58 @@ def wgs84_egm2008_crs() -> str:
 @pytest.fixture(scope='session')
 def rgb_byte_src_file(tmp_path_factory: pytest.TempPathFactory, im_size: tuple) -> Path:
     """An RGB byte checkerboard image with no CRS."""
+    array = checkerboard(im_size[::-1]).astype('uint8')
+    array = np.stack((array,) * 3, axis=0)
+    profile = create_profile(array)
     src_filename = tmp_path_factory.mktemp('data').joinpath('rgb_byte_src.tif')
-    create_src(src_filename, im_size, dtype='uint8', count=3)
+    with rio.open(src_filename, 'w', **profile) as im:
+        im.write(array)
     return src_filename
 
 
 @pytest.fixture(scope='session')
 def float_src_file(tmp_path_factory: pytest.TempPathFactory, im_size: tuple) -> Path:
     """A single band float64 checkerboard image with no CRS."""
+    array = np.expand_dims(checkerboard(im_size[::-1]).astype('float64'), axis=0)
+    profile = create_profile(array)
     src_filename = tmp_path_factory.mktemp('data').joinpath('float_src.tif')
-    create_src(src_filename, im_size, dtype='float32', count=1)
+    with rio.open(src_filename, 'w', **profile) as im:
+        im.write(array)
     return src_filename
 
 
 @pytest.fixture(scope='session')
 def rgb_byte_utm34n_src_file(
-    tmp_path_factory: pytest.TempPathFactory, pinhole_camera, utm34n_crs
+    tmp_path_factory: pytest.TempPathFactory, pinhole_camera: Camera, utm34n_crs: str
 ) -> Path:
     """An RGB byte checkerboard image with UTM zone 34N CRS and bounds 100m below
     ``pinhole_camera``.
     """
-    src_filename = tmp_path_factory.mktemp('data').joinpath('rgb_byte_src.tif')
-    create_src(
-        src_filename,
-        pinhole_camera._im_size,
-        dtype='uint8',
-        count=3,
-        camera=pinhole_camera,
-        crs=utm34n_crs,
-    )
+    array = checkerboard(pinhole_camera.im_size[::-1]).astype('uint8')
+    array = np.stack((array,) * 3, axis=0)
+    bounds = ortho_bounds(pinhole_camera)
+    transform = from_bounds(*bounds, *pinhole_camera.im_size)
+    profile = create_profile(array, crs=utm34n_crs, transform=transform, nodata=0)
+
+    src_filename = tmp_path_factory.mktemp('data').joinpath('rgb_byte_utm34n_src.tif')
+    with rio.open(src_filename, 'w', **profile) as im:
+        im.write(array)
     return src_filename
 
 
 @pytest.fixture(scope='session')
 def float_utm34n_dem_file(
-    tmp_path_factory: pytest.TempPathFactory, pinhole_camera, utm34n_crs
+    tmp_path_factory: pytest.TempPathFactory, pinhole_camera: Camera, utm34n_crs: str
 ) -> Path:
     """
     A 2 band float DEM file in UTM zone 34N with no vertical CRS.
 
-    Band 1 is a sinusoidal surface, and band 2, a planar surface.
+    Band 1 is a sinusoidal surface, and band 2, a horizontal plane.
     """
+    bounds = ortho_bounds(pinhole_camera)
+    array, transform = create_zsurf(bounds)
+    profile = create_profile(array, transform=transform, crs=utm34n_crs, nodata=float('nan'))
     filename = tmp_path_factory.mktemp('data').joinpath('float_utm34n_dem.tif')
-    array, profile = create_dem(
-        pinhole_camera, utm34n_crs, resolution=_dem_resolution, dtype='float32'
-    )
     with rio.open(filename, 'w', **profile) as im:
         im.write(array)
     return filename
@@ -453,17 +448,18 @@ def float_utm34n_dem_file(
 
 @pytest.fixture(scope='session')
 def float_utm34n_wgs84_dem_file(
-    tmp_path_factory: pytest.TempPathFactory, pinhole_camera, utm34n_wgs84_crs
+    tmp_path_factory: pytest.TempPathFactory, pinhole_camera: Camera, utm34n_wgs84_crs: str
 ) -> Path:
     """
     A 2 band float DEM file in UTM zone 34N with WGS84 ellipsoid vertical CRS.
 
-    Band 1 is a sinusoidal surface, and band 2, a planar surface.
+    Band 1 is a sinusoidal surface, and band 2, a horizontal plane.
     """
+    bounds = ortho_bounds(pinhole_camera)
+    array, transform = create_zsurf(bounds)
+    profile = create_profile(array, transform=transform, crs=utm34n_wgs84_crs, nodata=float('nan'))
+
     filename = tmp_path_factory.mktemp('data').joinpath('float_utm34n_wgs84_dem.tif')
-    array, profile = create_dem(
-        pinhole_camera, utm34n_wgs84_crs, resolution=_dem_resolution, dtype='float32'
-    )
     with rio.open(filename, 'w', **profile) as im:
         im.write(array)
     return filename
@@ -471,17 +467,18 @@ def float_utm34n_wgs84_dem_file(
 
 @pytest.fixture(scope='session')
 def float_utm34n_egm96_dem_file(
-    tmp_path_factory: pytest.TempPathFactory, pinhole_camera, utm34n_egm96_crs
+    tmp_path_factory: pytest.TempPathFactory, pinhole_camera: Camera, utm34n_egm96_crs: str
 ) -> Path:
     """
     A 2 band float DEM file in UTM zone 34N with EGM96 geoid vertical CRS.
 
-    Band 1 is a sinusoidal surface, and band 2, a planar surface.
+    Band 1 is a sinusoidal surface, and band 2, a horizontal plane.
     """
+    bounds = ortho_bounds(pinhole_camera)
+    array, transform = create_zsurf(bounds)
+    profile = create_profile(array, transform=transform, crs=utm34n_egm96_crs, nodata=float('nan'))
+
     filename = tmp_path_factory.mktemp('data').joinpath('float_utm34n_egm96_dem.tif')
-    array, profile = create_dem(
-        pinhole_camera, utm34n_egm96_crs, resolution=_dem_resolution, dtype='float32'
-    )
     with rio.open(filename, 'w', **profile) as im:
         im.write(array)
     return filename
@@ -489,17 +486,20 @@ def float_utm34n_egm96_dem_file(
 
 @pytest.fixture(scope='session')
 def float_utm34n_egm2008_dem_file(
-    tmp_path_factory: pytest.TempPathFactory, pinhole_camera, utm34n_egm2008_crs
+    tmp_path_factory: pytest.TempPathFactory, pinhole_camera: Camera, utm34n_egm2008_crs: str
 ) -> Path:
     """
     A 2 band float DEM file in UTM zone 34N with EGM2008 geoid vertical CRS.
 
-    Band 1 is a sinusoidal surface, and band 2, a planar surface.
+    Band 1 is a sinusoidal surface, and band 2, a horizontal plane.
     """
-    filename = tmp_path_factory.mktemp('data').joinpath('float_utm34n_egm96_dem.tif')
-    array, profile = create_dem(
-        pinhole_camera, utm34n_egm2008_crs, resolution=_dem_resolution, dtype='float32'
+    bounds = ortho_bounds(pinhole_camera)
+    array, transform = create_zsurf(bounds)
+    profile = create_profile(
+        array, transform=transform, crs=utm34n_egm2008_crs, nodata=float('nan')
     )
+
+    filename = tmp_path_factory.mktemp('data').joinpath('float_utm34n_egm96_dem.tif')
     with rio.open(filename, 'w', **profile) as im:
         im.write(array)
     return filename
@@ -507,21 +507,24 @@ def float_utm34n_egm2008_dem_file(
 
 @pytest.fixture(scope='session')
 def float_wgs84_wgs84_dem_file(
-    tmp_path_factory: pytest.TempPathFactory, pinhole_camera, utm34n_crs
+    tmp_path_factory: pytest.TempPathFactory,
+    pinhole_camera: Camera,
+    utm34n_crs: str,
+    wgs84_wgs84_crs: str,
 ) -> Path:
     """
     A 2 band float DEM file in WGS84 with WGS84 ellipsoid vertical CRS.
 
-    Band 1 is a sinusoidal surface, and band 2, a planar surface.
+    Band 1 is a sinusoidal surface, and band 2, a horizontal plane.
     """
+    bounds = ortho_bounds(pinhole_camera)
+    array, transform = create_zsurf(bounds)
+    bounds = array_bounds(*array.shape[-2:], transform)
+    bounds = transform_bounds(utm34n_crs, wgs84_wgs84_crs, *bounds)
+    transform = from_bounds(*bounds, *array.shape[-2::-1])
+    profile = create_profile(array, transform=transform, crs=wgs84_wgs84_crs, nodata=float('nan'))
+
     filename = tmp_path_factory.mktemp('data').joinpath('float_wgs84_wgs84_dem.tif')
-    array, profile = create_dem(
-        pinhole_camera,
-        utm34n_crs,
-        resolution=_dem_resolution,
-        dtype='float32',
-        dem_crs='EPSG:4326+4326',
-    )
     with rio.open(filename, 'w', **profile) as im:
         im.write(array)
     return filename
@@ -529,20 +532,22 @@ def float_wgs84_wgs84_dem_file(
 
 @pytest.fixture(scope='session')
 def float_utm34n_partial_dem_file(
-    tmp_path_factory: pytest.TempPathFactory, pinhole_camera, utm34n_crs
+    tmp_path_factory: pytest.TempPathFactory, pinhole_camera: Camera, utm34n_crs: str
 ) -> Path:
     """
     A 2 band float DEM file in UTM zone 34N with no vertical CRS.
 
-    Pixels above the diagonal are nodata. Band 1 is a sinusoidal surface, and band 2, a planar
-    surface.
+    Pixels above the diagonal are nodata. Band 1 is a sinusoidal surface, and band 2,
+    a horizontal plane.
     """
-    filename = tmp_path_factory.mktemp('data').joinpath('float_utm34n_dem.tif')
-    array, profile = create_dem(
-        pinhole_camera, utm34n_crs, resolution=_dem_resolution, dtype='float32'
-    )
+    nodata = float('nan')
+    bounds = ortho_bounds(pinhole_camera)
+    array, transform = create_zsurf(bounds)
     mask = np.fliplr(np.tril(np.ones(array.shape, dtype='bool'), k=1))
-    array[mask] = profile['nodata']
+    array[mask] = nodata
+    profile = create_profile(array, transform=transform, crs=utm34n_crs, nodata=nodata)
+
+    filename = tmp_path_factory.mktemp('data').joinpath('float_utm34n_dem.tif')
     with rio.open(filename, 'w', **profile) as im:
         im.write(array)
     return filename
@@ -552,7 +557,7 @@ def float_utm34n_partial_dem_file(
 def rgb_pinhole_utm34n_ortho(
     rgb_byte_src_file: Path,
     float_utm34n_dem_file: Path,
-    pinhole_camera: FrameCamera,
+    pinhole_camera: Camera,
     utm34n_crs: str,
 ) -> Ortho:
     """An Ortho object initialised with RGB byte source image, float DEM in UTM zone 34N (no