diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml
index fc7ab0b..05f9d2f 100644
--- a/.github/workflows/build.yml
+++ b/.github/workflows/build.yml
@@ -23,7 +23,7 @@ jobs:
       - name: Set up Java
         uses: actions/setup-java@v4
         with:
-          java-version: '8'
+          java-version: '11'
           distribution: 'zulu'
           cache: 'maven'
       - name: Set up CI environment
diff --git a/README.md b/README.md
index 29cafee..8d1bf39 100644
--- a/README.md
+++ b/README.md
@@ -133,7 +133,7 @@ GeffNode node0 = new GeffNode.Builder()
     .segmentId(0)
     .color(new double[]{1.0, 0.0, 0.0, 1.0}) // Red color
     .radius(2.5)
-    .covariance2d(new double[]{1.0, 0.2, 0.2, 1.5}) // 2x2 covariance matrix flattened
+    .covariance2d(new double[]{1.0, 0.2, 0.2, 1.5}) // 2x2 covariance matrix, row-major
     .polygonX(new double[]{1.0, 2.0, 3.0, 4.0}) // Polygon X coordinates
     .polygonY(new double[]{5.0, 6.0, 7.0, 8.0}) // Polygon Y coordinates
     .build();
@@ -283,8 +283,8 @@ dataset.zarr/
     │   │   ├── color/values        # RGBA colors [N, 4] (optional)
     │   │   ├── radius/values       # Node radii [N] (optional)
     │   │   ├── <tracklet>/values   # Track IDs [N] (name from track_node_props, optional)
-    │   │   ├── covariance2d/values # Flattened 2D covariance [N, 4] (optional)
-    │   │   ├── covariance3d/values # Flattened 3D covariance [N, 6] (optional)
+    │   │   ├── covariance2d/values # 2D covariance matrices [N, 2, 2] (optional)
+    │   │   ├── covariance3d/values # 3D covariance matrices [N, 3, 3] (optional)
     │   │   └── <varlength_prop>/   # Variable-length property, e.g. polygon (optional)
     │   │       ├── data            # Flattened values [V]
     │   │       ├── values          # Offsets and shapes [N, ndim+1]
diff --git a/src/main/java/org/mastodon/geff/Geff.java b/src/main/java/org/mastodon/geff/Geff.java
index 4ed4a76..e45ece5 100644
--- a/src/main/java/org/mastodon/geff/Geff.java
+++ b/src/main/java/org/mastodon/geff/Geff.java
@@ -82,7 +82,7 @@ public static void main( String[] args )
 
 			// Demonstrate reading edges
 			System.out.println( "\n=== Reading Edges ===" );
-			List< GeffEdge > edges = GeffEdge.readFromZarr( zarrPath, metadata.getGeffVersion() );
+			List< GeffEdge > edges = GeffEdge.readFromZarr( zarrPath, metadata );
 			System.out.println( "Read " + edges.size() + " edges:" );
 			for ( int i = 0; i < Math.min( 5, edges.size() ); i++ )
 			{
diff --git a/src/main/java/org/mastodon/geff/GeffEdge.java b/src/main/java/org/mastodon/geff/GeffEdge.java
index b6cfc2b..7d0f832 100644
--- a/src/main/java/org/mastodon/geff/GeffEdge.java
+++ b/src/main/java/org/mastodon/geff/GeffEdge.java
@@ -313,15 +313,26 @@ public static List< GeffEdge > readFromN5( final N5Reader reader, final String g
 		}
 		final int numEdges = edgeIds.size()[ 1 ];
 
+		// Determine the group names of the standard props from
+		// edge_props_metadata, falling back to the standard names for datasets
+		// that do not declare them.
+		final String distancePropName = metadata != null ? metadata.getEdgePropIdentifier( "distance" ) : "distance";
+		final String scorePropName = metadata != null ? metadata.getEdgePropIdentifier( "score" ) : "score";
+
 		// Read distances from chunks
-		final double[] distances = GeffUtils.readAsDoubleArray( reader, path + "/edges/props/distance/values", "distances" );
-		verifyLength( distances, numEdges, "/edges/props/distance/values" );
+		final String distancePropPath = path + "/edges/props/" + distancePropName + "/values";
+		final double[] distances = GeffUtils.readAsDoubleArray( reader, distancePropPath, "distances" );
+		verifyLength( distances, numEdges, distancePropPath );
 
 		// Read scores from chunks
-		final double[] scores = GeffUtils.readAsDoubleArray( reader, path + "/edges/props/score/values", "scores" );
-		verifyLength( scores, numEdges, "/edges/props/score/values" );
+		final String scorePropPath = path + "/edges/props/" + scorePropName + "/values";
+		final double[] scores = GeffUtils.readAsDoubleArray( reader, scorePropPath, "scores" );
+		verifyLength( scores, numEdges, scorePropPath );
 
 		// Read custom non-standard, non-varlength edge props from metadata
+		final Set< String > standardNames = new HashSet<>( STANDARD_EDGE_PROP_NAMES );
+		standardNames.add( distancePropName );
+		standardNames.add( scorePropName );
 		final Map< String, Object[] > customPropData = new HashMap<>();
 		final Map< String, VarlengthProperty[] > varlengthPropsMap = new HashMap<>();
 		if ( metadata != null && metadata.getEdgePropsMetadata() != null )
@@ -330,7 +341,7 @@ public static List< GeffEdge > readFromN5( final N5Reader reader, final String g
 			{
 				final String propName = entry.getKey();
 				final PropMetadata propMeta = entry.getValue();
-				if ( STANDARD_EDGE_PROP_NAMES.contains( propName ) )
+				if ( standardNames.contains( propName ) )
 					continue;
 				if ( GeffUtils.shouldSkipProperty( propName, propMeta ) )
 					continue;
@@ -375,7 +386,7 @@ public static List< GeffEdge > readFromN5( final N5Reader reader, final String g
 							}
 							customPropData.put( propName, edgeVals );
 						}
-						else if ( ndim == 2 )
+						else if ( ndim >= 2 )
 						{
 							if ( isFloat )
 							{
@@ -515,13 +526,19 @@ public static void writeToN5(
 		// there are no edge properties (mirrors what Python geff writes).
 		writer.createGroup( path + "/edges/props" );
 
+		// Determine the group names of the standard props from
+		// edge_props_metadata, falling back to the standard names for datasets
+		// that do not declare them.
+		final String distancePropName = metadata != null ? metadata.getEdgePropIdentifier( "distance" ) : "distance";
+		final String scorePropName = metadata != null ? metadata.getEdgePropIdentifier( "score" ) : "score";
+
 		// Write distances
-		if ( writeAllProps || edgePropsMetadata.containsKey( "distance" ) )
-			GeffUtils.writeDoubleArray( edges, GeffEdge::getDistance, writer, path + "/edges/props/distance/values", chunkSize );
+		if ( writeAllProps || metadata.getEdgePropMetadata( "distance" ) != null )
+			GeffUtils.writeDoubleArray( edges, GeffEdge::getDistance, writer, path + "/edges/props/" + distancePropName + "/values", chunkSize );
 
 		// Write scores
-		if ( writeAllProps || edgePropsMetadata.containsKey( "score" ) )
-			GeffUtils.writeDoubleArray( edges, GeffEdge::getScore, writer, path + "/edges/props/score/values", chunkSize );
+		if ( writeAllProps || metadata.getEdgePropMetadata( "score" ) != null )
+			GeffUtils.writeDoubleArray( edges, GeffEdge::getScore, writer, path + "/edges/props/" + scorePropName + "/values", chunkSize );
 
 		// When writeAllProps=true (no edgePropsMetadata provided), populate metadata
 		// with the standard props so the output zarr passes Python structural
@@ -529,8 +546,8 @@ public static void writeToN5(
 		if ( writeAllProps && metadata != null )
 		{
 			final Map< String, PropMetadata > edgePropsMap = new HashMap<>();
-			edgePropsMap.put( "distance", new PropMetadata( "distance", "float64", false, null, null, null ) );
-			edgePropsMap.put( "score", new PropMetadata( "score", "float64", false, null, null, null ) );
+			edgePropsMap.put( distancePropName, new PropMetadata( distancePropName, "float64", false, null, null, null ) );
+			edgePropsMap.put( scorePropName, new PropMetadata( scorePropName, "float64", false, null, null, null ) );
 			metadata.setEdgePropsMetadata( edgePropsMap );
 		}
 
diff --git a/src/main/java/org/mastodon/geff/GeffMetadata.java b/src/main/java/org/mastodon/geff/GeffMetadata.java
index 8af14b2..0f953c1 100644
--- a/src/main/java/org/mastodon/geff/GeffMetadata.java
+++ b/src/main/java/org/mastodon/geff/GeffMetadata.java
@@ -353,6 +353,165 @@ private String resolveSpaceAxisName( final String hint, final String standardNam
 		return ( spaceAxes.length == 0 || indexFromLast < 2 ) ? standardName : null;
 	}
 
+	/**
+	 * Get the metadata of a standard node property, as declared in
+	 * {@code node_props_metadata}.
+	 *
+	 * @param standardName
+	 *            the conventional name of the property, e.g. {@code "radius"}
+	 *            or {@code "covariance3d"}
+	 * @return the property metadata, or null if the property is not declared
+	 * @see #findPropMetadata(Map, String)
+	 */
+	public PropMetadata getNodePropMetadata( final String standardName )
+	{
+		return findPropMetadata( nodePropsMetadata, standardName );
+	}
+
+	/**
+	 * Get the metadata of a standard edge property, as declared in
+	 * {@code edge_props_metadata}.
+	 *
+	 * @param standardName
+	 *            the conventional name of the property, e.g. {@code "distance"}
+	 *            or {@code "score"}
+	 * @return the property metadata, or null if the property is not declared
+	 * @see #findPropMetadata(Map, String)
+	 */
+	public PropMetadata getEdgePropMetadata( final String standardName )
+	{
+		return findPropMetadata( edgePropsMetadata, standardName );
+	}
+
+	/**
+	 * Get the identifier of a standard node property, i.e. the name of the
+	 * group holding it in {@code nodes/props}.
+	 *
+	 * @param standardName
+	 *            the conventional name of the property, e.g. {@code "radius"}
+	 *            or {@code "covariance3d"}
+	 * @return the property identifier, never null
+	 * @see #resolvePropIdentifier(Map, String)
+	 */
+	public String getNodePropIdentifier( final String standardName )
+	{
+		return resolvePropIdentifier( nodePropsMetadata, standardName );
+	}
+
+	/**
+	 * Get the identifier of a standard edge property, i.e. the name of the
+	 * group holding it in {@code edges/props}.
+	 *
+	 * @param standardName
+	 *            the conventional name of the property, e.g. {@code "distance"}
+	 *            or {@code "score"}
+	 * @return the property identifier, never null
+	 * @see #resolvePropIdentifier(Map, String)
+	 */
+	public String getEdgePropIdentifier( final String standardName )
+	{
+		return resolvePropIdentifier( edgePropsMetadata, standardName );
+	}
+
+	/**
+	 * Look up a property in a props metadata map. The map is keyed by property
+	 * identifier, so a property is matched by the {@code identifier} of its
+	 * metadata first, and by the map key second, which covers maps built
+	 * without setting the identifier.
+	 *
+	 * @param propsMetadata
+	 *            {@code node_props_metadata} or {@code edge_props_metadata},
+	 *            may be null
+	 * @param standardName
+	 *            the conventional name of the property
+	 * @return the property metadata, or null if the property is not declared
+	 */
+	private static PropMetadata findPropMetadata( final Map< String, PropMetadata > propsMetadata, final String standardName )
+	{
+		if ( propsMetadata == null )
+			return null;
+
+		for ( final PropMetadata propMetadata : propsMetadata.values() )
+			if ( propMetadata != null && standardName.equals( propMetadata.getIdentifier() ) )
+				return propMetadata;
+
+		return propsMetadata.get( standardName );
+	}
+
+	/**
+	 * Resolve the identifier naming the group that holds a property, by looking
+	 * the property up in a props metadata map.
+	 *
+	 * @param propsMetadata
+	 *            {@code node_props_metadata} or {@code edge_props_metadata},
+	 *            may be null
+	 * @param standardName
+	 *            the conventional name of the property
+	 * @return the declared identifier, falling back to {@code standardName} for
+	 *         properties that are not declared, e.g. in datasets written before
+	 *         props metadata existed
+	 */
+	private static String resolvePropIdentifier( final Map< String, PropMetadata > propsMetadata, final String standardName )
+	{
+		final PropMetadata propMetadata = findPropMetadata( propsMetadata, standardName );
+		final String identifier = propMetadata != null ? propMetadata.getIdentifier() : null;
+		if ( identifier != null && !identifier.trim().isEmpty() )
+			return identifier;
+
+		LOG.debug( "no props metadata for '{}', assuming it is stored under that name", standardName );
+		return standardName;
+	}
+
+	/**
+	 * Get the positions that the x, y (and z) axes occupy in a property whose
+	 * dimensions follow the declaration order of the spatial axes, such as a
+	 * covariance matrix: the spec stores those "in the same coordinate system
+	 * as the {@code space} type properties", and axes are declared slowest to
+	 * fastest, i.e. typically {@code z, y, x}.
+	 * <p>
+	 * For axes declared as {@code z, y, x} this returns {@code { 2, 1, 0 }},
+	 * i.e. the x axis is the last dimension of the stored matrix. For axes
+	 * declared as {@code x, y, z} it returns the identity {@code { 0, 1, 2 }}.
+	 *
+	 * @param numDimensions
+	 *            2 for a 2D property (x, y), 3 for a 3D one (x, y, z)
+	 * @return the positions of the x, y (and z) axes in the stored order, or
+	 *         null if the dataset does not declare the spatial axes needed to
+	 *         tell, in which case the stored order is assumed to be x, y, z
+	 */
+	public int[] getSpaceAxisPositions( final int numDimensions )
+	{
+		final String[] spaceAxes = getAxisNamesByType( GeffAxis.TYPE_SPACE );
+		if ( spaceAxes.length < numDimensions )
+			return null;
+
+		final String[] axisNames = numDimensions == 3
+				? new String[] { getHorizontalAxisName(), getVerticalAxisName(), getDepthAxisName() }
+				: new String[] { getHorizontalAxisName(), getVerticalAxisName() };
+
+		final int[] declarationIndices = new int[ numDimensions ];
+		for ( int j = 0; j < numDimensions; ++j )
+		{
+			declarationIndices[ j ] = Arrays.asList( spaceAxes ).indexOf( axisNames[ j ] );
+			if ( declarationIndices[ j ] < 0 )
+			{
+				LOG.debug( "axis '{}' is not among the spatial axes {}, assuming the x, y, z order",
+						axisNames[ j ], Arrays.toString( spaceAxes ) );
+				return null;
+			}
+		}
+
+		// The position of an axis is its rank among the axes taken into
+		// account, so that a covariance2d of a z, y, x dataset is read as
+		// y, x rather than as the declaration indices 1, 2.
+		final int[] positions = new int[ numDimensions ];
+		for ( int j = 0; j < numDimensions; ++j )
+			for ( int k = 0; k < numDimensions; ++k )
+				if ( declarationIndices[ k ] < declarationIndices[ j ] )
+					++positions[ j ];
+		return positions;
+	}
+
 	public RelatedObjects getRelatedObjects()
 	{
 		return relatedObjects;
diff --git a/src/main/java/org/mastodon/geff/GeffNode.java b/src/main/java/org/mastodon/geff/GeffNode.java
index 33d2d41..b14f32a 100644
--- a/src/main/java/org/mastodon/geff/GeffNode.java
+++ b/src/main/java/org/mastodon/geff/GeffNode.java
@@ -93,9 +93,11 @@ public class GeffNode
 
 	public static final double DEFAULT_RADIUS = 1.0;
 
+	/** Identity, as the 4 elements of a row-major 2x2 matrix. */
 	public static final double[] DEFAULT_COVARIANCE_2D = { 1.0, 0.0, 0.0, 1.0 };
 
-	public static final double[] DEFAULT_COVARIANCE_3D = { 1.0, 0.0, 0.0, 1.0, 0.0, 1.0 };
+	/** Identity, as the 9 elements of a row-major 3x3 matrix. */
+	public static final double[] DEFAULT_COVARIANCE_3D = { 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0 };
 
 	/**
 	 * Default constructor
@@ -126,9 +128,11 @@ public GeffNode()
 	 * @param radius
 	 *            The radius of the node.
 	 * @param covariance2d
-	 *            The 2D covariance matrix of the node.
+	 *            The 2D covariance matrix of the node, as the 4 elements of a
+	 *            row-major 2x2 matrix.
 	 * @param covariance3d
-	 *            The 3D covariance matrix of the node.
+	 *            The 3D covariance matrix of the node, as the 9 elements of a
+	 *            row-major 3x3 matrix.
 	 * @param polygonX
 	 *            The x-coordinates of the polygon vertices.
 	 * @param polygonY
@@ -146,7 +150,7 @@ public GeffNode( int id, int timepoint, double x, double y, double z, double[] c
 		this.segmentId = segmentId;
 		this.radius = radius;
 		this.covariance2d = covariance2d != null ? covariance2d : DEFAULT_COVARIANCE_2D;
-		this.covariance3d = covariance3d != null ? covariance3d : DEFAULT_COVARIANCE_3D;
+		this.covariance3d = covariance3d != null ? asFullCovariance3d( covariance3d ) : DEFAULT_COVARIANCE_3D;
 		this.polygonX = polygonX != null ? polygonX : new double[ 0 ];
 		this.polygonY = polygonY != null ? polygonY : new double[ 0 ];
 		this.varlengthProps = new HashMap<>();
@@ -331,7 +335,10 @@ public void setRadius( double radius )
 	/**
 	 * Get the 2D covariance matrix of the node.
 	 *
-	 * @return The 2D covariance matrix as a 4-element array.
+	 * @return The 2D covariance matrix as the 4 elements of a row-major 2x2
+	 *         matrix in x, y order, i.e. {@code { xx, xy, yx, yy }}. On disk
+	 *         the matrix follows the declaration order of the spatial axes
+	 *         instead.
 	 */
 	public double[] getCovariance2d()
 	{
@@ -342,7 +349,10 @@ public double[] getCovariance2d()
 	 * Set the 2D covariance matrix of the node.
 	 *
 	 * @param covariance2d
-	 *            The 2D covariance matrix to set as a 4-element array.
+	 *            The 2D covariance matrix to set, as the 4 elements of a
+	 *            row-major 2x2 matrix, i.e. {@code { xx, xy, yx, yy }}. This
+	 *            matches the {@code [N, 2, 2]} shape the GEFF spec defines for
+	 *            the {@code covariance2d} property.
 	 *
 	 * @throws IllegalArgumentException
 	 *             if the covariance2d array is not of length 4.
@@ -362,7 +372,11 @@ public void setCovariance2d( double[] covariance2d )
 	/**
 	 * Get the 3D covariance matrix of the node.
 	 *
-	 * @return The 3D covariance matrix as a 6-element array.
+	 * @return The 3D covariance matrix as the 9 elements of a row-major 3x3
+	 *         matrix in x, y, z order, i.e.
+	 *         {@code { xx, xy, xz, yx, yy, yz, zx, zy, zz }}. On disk the
+	 *         matrix follows the declaration order of the spatial axes
+	 *         instead.
 	 */
 	public double[] getCovariance3d()
 	{
@@ -373,21 +387,97 @@ public double[] getCovariance3d()
 	 * Set the 3D covariance matrix of the node.
 	 *
 	 * @param covariance3d
-	 *            The 3D covariance matrix to set as a 6-element array.
+	 *            The 3D covariance matrix to set, as the 9 elements of a
+	 *            row-major 3x3 matrix, i.e.
+	 *            {@code { xx, xy, xz, yx, yy, yz, zx, zy, zz }}. This matches
+	 *            the {@code [N, 3, 3]} shape the GEFF spec defines for the
+	 *            {@code covariance3d} property. A 6-element array is also
+	 *            accepted and read as the upper triangle
+	 *            {@code { xx, xy, xz, yy, yz, zz }}, the layout geff-java
+	 *            wrote before the spec shape was implemented.
 	 *
 	 * @throws IllegalArgumentException
-	 *             if the covariance3d array is not of length 6.
+	 *             if the covariance3d array is neither of length 9 nor 6.
 	 */
 	public void setCovariance3d( double[] covariance3d )
 	{
+		this.covariance3d = asFullCovariance3d( covariance3d );
+	}
+
+	/**
+	 * Convert a 3D covariance matrix to the 9 elements of a row-major 3x3
+	 * matrix.
+	 *
+	 * @param covariance3d
+	 *            either the 9 elements of a row-major 3x3 matrix, or the 6
+	 *            upper-triangular elements
+	 *            {@code { xx, xy, xz, yy, yz, zz }} of a symmetric matrix
+	 * @return the 9 elements of the row-major 3x3 matrix
+	 * @throws IllegalArgumentException
+	 *             if the array is neither of length 9 nor 6
+	 */
+	static double[] asFullCovariance3d( final double[] covariance3d )
+	{
+		if ( covariance3d != null && covariance3d.length == 9 )
+			return covariance3d;
+
 		if ( covariance3d != null && covariance3d.length == 6 )
 		{
-			this.covariance3d = covariance3d;
+			final double xx = covariance3d[ 0 ];
+			final double xy = covariance3d[ 1 ];
+			final double xz = covariance3d[ 2 ];
+			final double yy = covariance3d[ 3 ];
+			final double yz = covariance3d[ 4 ];
+			final double zz = covariance3d[ 5 ];
+			return new double[] { xx, xy, xz, xy, yy, yz, xz, yz, zz };
 		}
-		else
+
+		throw new IllegalArgumentException( "Covariance3D must be a 9-element array (a row-major 3x3 matrix) "
+				+ "or a 6-element array (its upper triangle)" );
+	}
+
+	/**
+	 * Permute the dimensions of a covariance matrix between the order the
+	 * spatial axes are declared in, which the GEFF spec uses for the stored
+	 * matrix, and the x, y, z order this class exposes.
+	 *
+	 * @param values
+	 *            the elements of a row-major matrix, i.e. 4 values for a 2x2
+	 *            or 9 values for a 3x3 matrix
+	 * @param positions
+	 *            the positions of the x, y (and z) axes in the stored matrix,
+	 *            or null to leave the matrix untouched
+	 * @param storedToXyz
+	 *            true to convert a stored matrix into x, y, z order, false for
+	 *            the other direction
+	 * @return the permuted matrix, or {@code values} itself if nothing has to
+	 *         be permuted
+	 * @see GeffMetadata#getSpaceAxisPositions(int)
+	 */
+	static double[] permuteCovariance( final double[] values, final int[] positions, final boolean storedToXyz )
+	{
+		if ( values == null || positions == null )
+			return values;
+
+		final int n = positions.length;
+		if ( values.length != n * n )
 		{
-			throw new IllegalArgumentException( "Covariance3D must be a 6-element array" );
+			LOG.debug( "not permuting a {}-element covariance matrix of {} dimensions", values.length, n );
+			return values;
 		}
+
+		final double[] permuted = new double[ values.length ];
+		for ( int a = 0; a < n; ++a )
+			for ( int b = 0; b < n; ++b )
+			{
+				final int xyzIndex = a * n + b;
+				final int storedIndex = positions[ a ] * n + positions[ b ];
+				if ( storedToXyz )
+					permuted[ xyzIndex ] = values[ storedIndex ];
+				else
+					permuted[ storedIndex ] = values[ xyzIndex ];
+			}
+		return permuted;
 	}
 
 	/**
@@ -648,16 +738,15 @@ public Builder covariance2d( double[] covariance2d )
 			return this;
 		}
 
+		/**
+		 * @param covariance3d
+		 *            the 9 elements of a row-major 3x3 matrix, or the 6
+		 *            upper-triangular elements of a symmetric matrix
+		 * @see GeffNode#setCovariance3d(double[])
+		 */
 		public Builder covariance3d( double[] covariance3d )
 		{
-			if ( covariance3d != null && covariance3d.length == 6 )
-			{
-				this.covariance3d = covariance3d;
-			}
-			else
-			{
-				throw new IllegalArgumentException( "Covariance3D must be a 6-element array" );
-			}
+			this.covariance3d = asFullCovariance3d( covariance3d );
 			return this;
 		}
 
@@ -770,9 +859,18 @@ public static List< GeffNode > readFromN5( final N5Reader reader, final String g
 			zCoords = null;
 		}
 
+		// Determine the group names of the non-axis standard props from
+		// node_props_metadata, falling back to the standard names for datasets
+		// that do not declare them.
+		final String colorPropName = metadata.getNodePropIdentifier( "color" );
+		final String radiusPropName = metadata.getNodePropIdentifier( "radius" );
+		final String covariance2dPropName = metadata.getNodePropIdentifier( "covariance2d" );
+		final String covariance3dPropName = metadata.getNodePropIdentifier( "covariance3d" );
+
 		// Read color from chunks
-		final FlattenedDoubles colors = GeffUtils.readAsDoubleMatrix( reader, path + "/nodes/props/color/values", "color" );
-		verifyLength( colors, numNodes, "/nodes/props/color/values" );
+		final String colorPropPath = path + "/nodes/props/" + colorPropName + "/values";
+		final FlattenedDoubles colors = GeffUtils.readAsDoubleMatrix( reader, colorPropPath, "color" );
+		verifyLength( colors, numNodes, colorPropPath );
 
 		// Read track IDs from chunks
 		final String trackletProp = metadata.getTrackNodeProps() != null && metadata.getTrackNodeProps().containsKey( "tracklet" )
@@ -794,16 +892,35 @@ public static List< GeffNode > readFromN5( final N5Reader reader, final String g
 		}
 
 		// Read radius from chunks
-		double[] radius = GeffUtils.readAsDoubleArray( reader, path + "/nodes/props/radius/values", "radius" );
-		verifyLength( radius, numNodes, "/nodes/props/radius/values" );
-
-		// Read covariance2d from chunks
-		final FlattenedDoubles covariance2ds = GeffUtils.readAsDoubleMatrix( reader, path + "/nodes/props/covariance2d/values", "covariance2d" );
-		verifyLength( covariance2ds, numNodes, "/nodes/props/covariance2d/values" );
-
-		// Read covariance3d from chunks
-		final FlattenedDoubles covariance3ds = GeffUtils.readAsDoubleMatrix( reader, path + "/nodes/props/covariance3d/values", "covariance3d" );
-		verifyLength( covariance3ds, numNodes, "/nodes/props/covariance3d/values" );
+		final String radiusPropPath = path + "/nodes/props/" + radiusPropName + "/values";
+		final double[] radius = GeffUtils.readAsDoubleArray( reader, radiusPropPath, "radius" );
+		verifyLength( radius, numNodes, radiusPropPath );
+
+		// Read covariance2d from chunks. The spec shape [N, 2, 2] flattens to
+		// the same 4 values per node as the [N, 4] arrays geff-java wrote
+		// before.
+		final String covariance2dPropPath = path + "/nodes/props/" + covariance2dPropName + "/values";
+		final FlattenedDoubles covariance2ds = GeffUtils.readAsDoubleMatrix( reader, covariance2dPropPath, "covariance2d" );
+		verifyLength( covariance2ds, numNodes, covariance2dPropPath );
+		GeffUtils.verifyNumValues( covariance2ds, covariance2dPropPath, 4 );
+
+		// Read covariance3d from chunks. The spec shape [N, 3, 3] gives 9
+		// values per node; the [N, 6] arrays geff-java wrote before give the 6
+		// upper-triangular values, which are expanded below.
+		final String covariance3dPropPath = path + "/nodes/props/" + covariance3dPropName + "/values";
+		final FlattenedDoubles covariance3ds = GeffUtils.readAsDoubleMatrix( reader, covariance3dPropPath, "covariance3d" );
+		verifyLength( covariance3ds, numNodes, covariance3dPropPath );
+		GeffUtils.verifyNumValues( covariance3ds, covariance3dPropPath, 9, 6 );
+
+		// The stored matrices follow the declaration order of the spatial axes
+		// (typically z, y, x), while this class exposes them in x, y, z order.
+		// Only the [N, 2, 2] and [N, 3, 3] shapes of the spec are permuted: the
+		// flat [N, 4] and [N, 6] arrays written by earlier versions of
+		// geff-java hold whatever order the caller passed in, i.e. x, y, z.
+		final int[] covariance2dPositions = covariance2ds != null && covariance2ds.size().length == 3
+				? metadata.getSpaceAxisPositions( 2 ) : null;
+		final int[] covariance3dPositions = covariance3ds != null && covariance3ds.size().length == 3
+				? metadata.getSpaceAxisPositions( 3 ) : null;
 
 		// Read polygon from chunks
 		double[][] polygonsX = null;
@@ -832,6 +949,10 @@ public static List< GeffNode > readFromN5( final N5Reader reader, final String g
 		// Read custom non-standard, non-varlength node props from metadata
 		final Set< String > standardNames = new HashSet<>( STANDARD_NODE_PROP_NAMES );
 		standardNames.add( trackletProp );
+		standardNames.addAll( Arrays.asList( timeAxisName, xAxisName, yAxisName, colorPropName,
+				radiusPropName, covariance2dPropName, covariance3dPropName ) );
+		if ( zAxisName != null )
+			standardNames.add( zAxisName );
 		final Map< String, Object[] > customPropData = new HashMap<>();
 		if ( metadata.getNodePropsMetadata() != null )
 		{
@@ -871,7 +992,7 @@ public static List< GeffNode > readFromN5( final N5Reader reader, final String g
 						}
 						customPropData.put( propName, nodeVals );
 					}
-					else if ( ndim == 2 )
+					else if ( ndim >= 2 )
 					{
 						if ( isFloat )
 						{
@@ -978,8 +1099,10 @@ else if ( ndim == 2 )
 			final double[] color = colors != null ? colors.rowAt( i ) : DEFAULT_COLOR;
 			final int segmentId = trackIds != null ? trackIds[ i ] : -1;
 			final double r = radius != null ? radius[ i ] : Double.NaN;
-			final double[] covariance2d = covariance2ds != null ? covariance2ds.rowAt( i ) : DEFAULT_COVARIANCE_2D;
-			final double[] covariance3d = covariance3ds != null ? covariance3ds.rowAt( i ) : DEFAULT_COVARIANCE_3D;
+			final double[] covariance2d = covariance2ds != null
+					? permuteCovariance( covariance2ds.rowAt( i ), covariance2dPositions, true ) : DEFAULT_COVARIANCE_2D;
+			final double[] covariance3d = covariance3ds != null
+					? permuteCovariance( asFullCovariance3d( covariance3ds.rowAt( i ) ), covariance3dPositions, true ) : DEFAULT_COVARIANCE_3D;
 			final double[] polygonX = polygonsX != null ? polygonsX[ i ] : null;
 			final double[] polygonY = polygonsY != null ? polygonsY[ i ] : null;
 			final GeffNode node = new GeffNode( id, t, x, y, z, color, segmentId, r, covariance2d, covariance3d, polygonX, polygonY );
@@ -1107,9 +1230,17 @@ public static void writeToN5(
 		if ( zAxisName != null && ( writeAllProps || metadataNodeProps.containsKey( zAxisName ) ) )
 			GeffUtils.writeDoubleArray( nodes, GeffNode::getZ, writer, path + "/nodes/props/" + zAxisName + "/values", chunkSize );
 
+		// Determine the group names of the non-axis standard props from
+		// node_props_metadata, falling back to the standard names for datasets
+		// that do not declare them.
+		final String colorPropName = metadata.getNodePropIdentifier( "color" );
+		final String radiusPropName = metadata.getNodePropIdentifier( "radius" );
+		final String covariance2dPropName = metadata.getNodePropIdentifier( "covariance2d" );
+		final String covariance3dPropName = metadata.getNodePropIdentifier( "covariance3d" );
+
 		// Write color in chunks
-		if ( writeAllProps || metadataNodeProps.containsKey( "color" ) )
-			GeffUtils.writeDoubleMatrix( nodes, 4, GeffNode::getColor, writer, path + "/nodes/props/color/values", chunkSize );
+		if ( writeAllProps || metadata.getNodePropMetadata( "color" ) != null )
+			GeffUtils.writeDoubleMatrix( nodes, 4, GeffNode::getColor, writer, path + "/nodes/props/" + colorPropName + "/values", chunkSize );
 
 		// Write segment IDs in chunks
 		final String trackletProp = metadata.getTrackNodeProps() != null && metadata.getTrackNodeProps().containsKey( "tracklet" )
@@ -1118,16 +1249,25 @@ public static void writeToN5(
 			GeffUtils.writeIntArray( nodes, GeffNode::getSegmentId, writer, path + "/nodes/props/" + trackletProp + "/values", chunkSize );
 
 		// Write radius and covariance attributes if available
-		if ( writeAllProps || metadataNodeProps.containsKey( "radius" ) )
-			GeffUtils.writeDoubleArray( nodes, GeffNode::getRadius, writer, path + "/nodes/props/radius/values", chunkSize );
-
-		// Write covariance2d in chunks
-		if ( writeAllProps || metadataNodeProps.containsKey( "covariance2d" ) )
-			GeffUtils.writeDoubleMatrix( nodes, 4, GeffNode::getCovariance2d, writer, path + "/nodes/props/covariance2d/values", chunkSize );
-
-		// Write covariance3d in chunks
-		if ( writeAllProps || metadataNodeProps.containsKey( "covariance3d" ) )
-			GeffUtils.writeDoubleMatrix( nodes, 6, GeffNode::getCovariance3d, writer, path + "/nodes/props/covariance3d/values", chunkSize );
+		if ( writeAllProps || metadata.getNodePropMetadata( "radius" ) != null )
+			GeffUtils.writeDoubleArray( nodes, GeffNode::getRadius, writer, path + "/nodes/props/" + radiusPropName + "/values", chunkSize );
+
+		// Write the covariances in chunks, with the shapes [N, 2, 2] and
+		// [N, 3, 3] the spec defines, and with their dimensions in the
+		// declaration order of the spatial axes (typically z, y, x) instead of
+		// the x, y, z order this class exposes.
+		final int[] covariance2dPositions = metadata.getSpaceAxisPositions( 2 );
+		final int[] covariance3dPositions = metadata.getSpaceAxisPositions( 3 );
+
+		if ( writeAllProps || metadata.getNodePropMetadata( "covariance2d" ) != null )
+			GeffUtils.writeDoubleMatrix( nodes, new int[] { 2, 2 },
+					node -> permuteCovariance( node.getCovariance2d(), covariance2dPositions, false ),
+					writer, path + "/nodes/props/" + covariance2dPropName + "/values", chunkSize );
+
+		if ( writeAllProps || metadata.getNodePropMetadata( "covariance3d" ) != null )
+			GeffUtils.writeDoubleMatrix( nodes, new int[] { 3, 3 },
+					node -> permuteCovariance( node.getCovariance3d(), covariance3dPositions, false ),
+					writer, path + "/nodes/props/" + covariance3dPropName + "/values", chunkSize );
 
 		// When writeAllProps=true (no nodePropsMetadata provided), populate metadata
 		// with the standard props so the output zarr passes Python structural
@@ -1144,11 +1284,11 @@ public static void writeToN5(
 			{
 				nodePropsMap.put( zAxisName, new PropMetadata( zAxisName, "float64", false, null, null, null ) );
 			}
-			nodePropsMap.put( "color", new PropMetadata( "color", "float64", false, null, null, null ) );
+			nodePropsMap.put( colorPropName, new PropMetadata( colorPropName, "float64", false, null, null, null ) );
 			nodePropsMap.put( trackletProp, new PropMetadata( trackletProp, "int32", false, null, null, null ) );
-			nodePropsMap.put( "radius", new PropMetadata( "radius", "float64", false, null, null, null ) );
-			nodePropsMap.put( "covariance2d", new PropMetadata( "covariance2d", "float64", false, null, null, null ) );
-			nodePropsMap.put( "covariance3d", new PropMetadata( "covariance3d", "float64", false, null, null, null ) );
+			nodePropsMap.put( radiusPropName, new PropMetadata( radiusPropName, "float64", false, null, null, null ) );
+			nodePropsMap.put( covariance2dPropName, new PropMetadata( covariance2dPropName, "float64", false, null, null, null ) );
+			nodePropsMap.put( covariance3dPropName, new PropMetadata( covariance3dPropName, "float64", false, null, null, null ) );
 			metadata.setNodePropsMetadata( nodePropsMap );
 		}
 
diff --git a/src/main/java/org/mastodon/geff/GeffUtils.java b/src/main/java/org/mastodon/geff/GeffUtils.java
index 01903dc..387d5a7 100644
--- a/src/main/java/org/mastodon/geff/GeffUtils.java
+++ b/src/main/java/org/mastodon/geff/GeffUtils.java
@@ -272,6 +272,30 @@ public static < T > void writeDoubleMatrix(
 			final String dataset,
 			final int chunkSize )
 	{
+		writeDoubleMatrix( elements, new int[] { numColumns }, extractor, writer, dataset, chunkSize );
+	}
+
+	/**
+	 * Write one {@code double[]} per element as a dataset of zarr shape
+	 * {@code [numElements, rowShape...]}, e.g. {@code [N, 3, 3]} for
+	 * {@code rowShape = { 3, 3 }}. The values of a row are stored in the order
+	 * the extracted array holds them, i.e. row-major for a matrix.
+	 *
+	 * @param rowShape
+	 *            shape of the values belonging to a single element
+	 */
+	public static < T > void writeDoubleMatrix(
+			final List< T > elements,
+			final int[] rowShape,
+			final Function< T, double[] > extractor,
+			final N5Writer writer,
+			final String dataset,
+			final int chunkSize )
+	{
+		int numColumns = 1;
+		for ( final int d : rowShape )
+			numColumns *= d;
+
 		final int size = elements.size();
 		final double[] data = new double[ numColumns * size ];
 		for ( int i = 0; i < size; ++i )
@@ -281,11 +305,35 @@ public static < T > void writeDoubleMatrix(
 				continue;
 			System.arraycopy( row, 0, data, numColumns * i, numColumns );
 		}
-			final DatasetAttributes attributes = new DatasetAttributes(
-					new long[] { numColumns, size },
-					new int[] { numColumns, chunkSize },
-					DataType.FLOAT64,
-					defaultCompression() );
+
+		// N5 dimensions are in the reverse order of the zarr shape, so the
+		// element dimension comes last.
+		final int n = rowShape.length + 1;
+		final long[] dimensions = new long[ n ];
+		final int[] blockSize = new int[ n ];
+		for ( int d = 0; d < rowShape.length; ++d )
+		{
+			dimensions[ d ] = rowShape[ rowShape.length - 1 - d ];
+			blockSize[ d ] = rowShape[ rowShape.length - 1 - d ];
+		}
+		dimensions[ n - 1 ] = size;
+		blockSize[ n - 1 ] = chunkSize;
+
+		// Overwriting a dataset of a different shape would leave the chunks of
+		// the old shape behind, e.g. when re-writing a covariance3d that was
+		// stored as [N, 6] before.
+		if ( writer.datasetExists( dataset )
+				&& !Arrays.equals( writer.getDatasetAttributes( dataset ).getDimensions(), dimensions ) )
+		{
+			LOG.debug( "removing dataset {} to re-create it with the dimensions {}", dataset, Arrays.toString( dimensions ) );
+			writer.remove( dataset );
+		}
+
+		final DatasetAttributes attributes = new DatasetAttributes(
+				dimensions,
+				blockSize,
+				DataType.FLOAT64,
+				defaultCompression() );
 		writer.createDataset( dataset, attributes );
 		write( data, writer, dataset, attributes );
 	}
@@ -342,18 +390,36 @@ int[] size()
 			return data[ i0 + size[ 0 ] * i1 ];
 		}
 
-		// TODO: remove until needed
 		double at( final int i0, final int i1, final int i2 )
 		{
 			assert size.length == 3;
-			return data[ i0 + size[ 0 ] * ( i1 * i2 * size[ 1 ] ) ];
+			return data[ i0 + size[ 0 ] * ( i1 + size[ 1 ] * i2 ) ];
 		}
 
-		double[] rowAt( final int i1 )
+		/**
+		 * Number of values a single row holds, i.e. the product of all
+		 * dimensions but the last one. For a property stored with the zarr
+		 * shape {@code [N, 3, 3]} this is 9.
+		 */
+		int numValuesPerRow()
 		{
-			assert size.length == 2;
-			final double[] row = new double[ size[ 0 ] ];
-			Arrays.setAll( row, i0 -> at( i0, i1 ) );
+			int numValues = 1;
+			for ( int d = 0; d < size.length - 1; ++d )
+				numValues *= size[ d ];
+			return numValues;
+		}
+
+		/**
+		 * Get the values of one row, i.e. of one node or edge. The values of a
+		 * row are contiguous, so this works for any number of dimensions: a
+		 * property stored with the zarr shape {@code [N, 3, 3]} yields the 9
+		 * values of its 3x3 matrix, in row-major order.
+		 */
+		double[] rowAt( final int i )
+		{
+			final int numValues = numValuesPerRow();
+			final double[] row = new double[ numValues ];
+			System.arraycopy( data, i * numValues, row, 0, numValues );
 			return row;
 		}
 	}
@@ -366,8 +432,8 @@ public static FlattenedDoubles readAsDoubleMatrix( final N5Reader reader, final
 			return null;
 		}
 		final DatasetAttributes attributes = reader.getDatasetAttributes( dataset );
-		if ( attributes.getNumDimensions() != 2 )
-		{ throw new IllegalArgumentException( "Expected 2D array" ); }
+		if ( attributes.getNumDimensions() < 2 )
+		{ throw new IllegalArgumentException( "Expected an array of 2 or more dimensions" ); }
 		return new FlattenedDoubles( convertToDoubleArray( readFully( reader, dataset ), description ), attributes.getDimensions() );
 	}
 
@@ -399,11 +465,27 @@ int at( final int i0, final int i1 )
 			return data[ i0 + size[ 0 ] * i1 ];
 		}
 
-		int[] rowAt( final int i0 )
+		/**
+		 * Number of values a single row holds, i.e. the product of all
+		 * dimensions but the last one.
+		 */
+		int numValuesPerRow()
 		{
-			assert size.length == 2;
-			final int[] row = new int[ size[ 1 ] ];
-			Arrays.setAll( row, i1 -> at( i0, i1 ) );
+			int numValues = 1;
+			for ( int d = 0; d < size.length - 1; ++d )
+				numValues *= size[ d ];
+			return numValues;
+		}
+
+		/**
+		 * Get the values of one row, i.e. of one node or edge. The values of a
+		 * row are contiguous, so this works for any number of dimensions.
+		 */
+		int[] rowAt( final int i )
+		{
+			final int numValues = numValuesPerRow();
+			final int[] row = new int[ numValues ];
+			System.arraycopy( data, i * numValues, row, 0, numValues );
 			return row;
 		}
 	}
@@ -416,8 +498,8 @@ public static FlattenedInts readAsIntMatrix( final N5Reader reader, final String
 			return null;
 		}
 		final DatasetAttributes attributes = reader.getDatasetAttributes( dataset );
-		if ( attributes.getNumDimensions() != 2 )
-		{ throw new IllegalArgumentException( "Expected 2D array" ); }
+		if ( attributes.getNumDimensions() < 2 )
+		{ throw new IllegalArgumentException( "Expected an array of 2 or more dimensions" ); }
 		return new FlattenedInts( convertToIntArray( readFully( reader, dataset ), description ), attributes.getDimensions() );
 	}
 
@@ -507,6 +589,29 @@ public static void verifyLength( final FlattenedInts array, final int expectedLe
 		{ throw new IllegalArgumentException( "property " + name + " does not have expected length (" + array.size()[ array.size().length - 1 ] + " vs " + expectedLength + ")" ); }
 	}
 
+	/**
+	 * Verify that a property holds one of the accepted numbers of values per
+	 * node or edge, regardless of how those values are shaped: a covariance3d
+	 * holds 9 values whether it is stored as {@code [N, 3, 3]} or as
+	 * {@code [N, 9]}.
+	 *
+	 * @param acceptedNumValues
+	 *            the accepted numbers of values per node or edge
+	 */
+	public static void verifyNumValues( final FlattenedDoubles array, final String name, final int... acceptedNumValues )
+	{
+		if ( array == null )
+			return;
+
+		final int numValues = array.numValuesPerRow();
+		for ( final int accepted : acceptedNumValues )
+			if ( numValues == accepted )
+				return;
+
+		throw new IllegalArgumentException( "property " + name + " does not have an expected number of values per element ("
+				+ numValues + " vs " + Arrays.toString( acceptedNumValues ) + ", dimensions " + Arrays.toString( array.size() ) + ")" );
+	}
+
 	// -- write dataset fully --
 
 	public static void write(
diff --git a/src/main/java/org/mastodon/geff/RoundTripGeff.java b/src/main/java/org/mastodon/geff/RoundTripGeff.java
index 9a7f1af..699992f 100644
--- a/src/main/java/org/mastodon/geff/RoundTripGeff.java
+++ b/src/main/java/org/mastodon/geff/RoundTripGeff.java
@@ -68,7 +68,7 @@ public static void main( String[] args )
 
 			// Read edges
 			System.out.println( "\nReading edges..." );
-			List< GeffEdge > edges = GeffEdge.readFromZarr( inputPath, metadata.getGeffVersion() );
+			List< GeffEdge > edges = GeffEdge.readFromZarr( inputPath, metadata );
 			System.out.println( "  Read " + edges.size() + " edges" );
 
 			// Write nodes first so that metadata can be updated (e.g. varlength
diff --git a/src/test/java/org/mastodon/geff/GeffTest.java b/src/test/java/org/mastodon/geff/GeffTest.java
index 9d86745..bbaa2a6 100644
--- a/src/test/java/org/mastodon/geff/GeffTest.java
+++ b/src/test/java/org/mastodon/geff/GeffTest.java
@@ -36,8 +36,15 @@
 import org.junit.jupiter.api.io.TempDir;
 import static org.junit.jupiter.api.Assertions.*;
 
+import org.janelia.saalfeldlab.n5.zarr.N5ZarrReader;
+import org.janelia.saalfeldlab.n5.zarr.N5ZarrWriter;
+
+import java.nio.file.Files;
 import java.nio.file.Path;
+import java.nio.file.Paths;
 import java.util.ArrayList;
+import java.util.HashMap;
+import java.util.Map;
 
 /**
  * Test class for the main Geff functionality
@@ -350,6 +357,148 @@ void testCovariance2dRoundTrip( @TempDir Path tempDir ) throws IOException
 			assertArrayEquals( cov2dValues[ i ], readNodes.get( i ).getCovariance2d(), 1e-9,
 					"covariance2d mismatch at node " + i );
 		}
+
+		// The GEFF spec defines covariance2d as [N, 2, 2]
+		assertArrayEquals( new long[] { cov2dValues.length, 2, 2 },
+				zarrShape( tempPath + "/nodes/props/covariance2d/values" ) );
+	}
+
+	/**
+	 * Read the zarr shape of a written dataset, i.e. the {@code shape} entry of
+	 * its {@code .zarray}.
+	 */
+	private static long[] zarrShape( final String datasetPath ) throws IOException
+	{
+		try (final N5ZarrReader reader = new N5ZarrReader( datasetPath, true ))
+		{
+			// N5 reports the dimensions in the reverse order of the zarr shape
+			final long[] dimensions = reader.getDatasetAttributes( "/" ).getDimensions();
+			final long[] shape = new long[ dimensions.length ];
+			for ( int d = 0; d < dimensions.length; d++ )
+				shape[ d ] = dimensions[ dimensions.length - 1 - d ];
+			return shape;
+		}
+	}
+
+	@Test
+	@DisplayName( "Test covariance dimensions follow the z, y, x axes order on disk" )
+	void testCovarianceAxisOrder( @TempDir Path tempDir ) throws IOException
+	{
+		final String tempPath = tempDir.toString() + "/test-cov-axis-order.zarr/tracks";
+
+		// Axes in the order the spec defines, slowest to fastest
+		final GeffMetadata metadata = new GeffMetadata( Geff.VERSION, true );
+		metadata.setGeffAxes( new GeffAxis[] {
+				new GeffAxis( "t", GeffAxis.TYPE_TIME, GeffAxis.UNIT_SECOND, 0.0, 10.0 ),
+				new GeffAxis( "z", GeffAxis.TYPE_SPACE, GeffAxis.UNIT_MICROMETER, 0.0, 10.0 ),
+				new GeffAxis( "y", GeffAxis.TYPE_SPACE, GeffAxis.UNIT_MICROMETER, 0.0, 10.0 ),
+				new GeffAxis( "x", GeffAxis.TYPE_SPACE, GeffAxis.UNIT_MICROMETER, 0.0, 10.0 )
+		} );
+		assertArrayEquals( new int[] { 2, 1, 0 }, metadata.getSpaceAxisPositions( 3 ) );
+		assertArrayEquals( new int[] { 1, 0 }, metadata.getSpaceAxisPositions( 2 ) );
+
+		// A covariance whose entries are all distinct, so that a wrong
+		// permutation cannot go unnoticed
+		final double[] cov3d = { 1.0, 2.0, 3.0, 2.0, 4.0, 5.0, 3.0, 5.0, 6.0 };
+		final double[] cov2d = { 1.0, 2.0, 2.0, 4.0 };
+
+		final GeffNode node = new GeffNode();
+		node.setId( 0 );
+		node.setT( 0 );
+		node.setX( 1.0 );
+		node.setY( 2.0 );
+		node.setZ( 3.0 );
+		node.setCovariance3d( cov3d );
+		node.setCovariance2d( cov2d );
+		final List< GeffNode > nodes = new ArrayList<>();
+		nodes.add( node );
+
+		GeffNode.writeToZarr( nodes, tempPath, metadata );
+		GeffMetadata.writeToZarr( metadata, tempPath );
+
+		// On disk the dimensions are z, y, x, i.e. the x, y, z matrix with its
+		// rows and columns reversed
+		final double[] storedCov3d = { 6.0, 5.0, 3.0, 5.0, 4.0, 2.0, 3.0, 2.0, 1.0 };
+		final double[] storedCov2d = { 4.0, 2.0, 2.0, 1.0 };
+		try (final N5ZarrReader reader = new N5ZarrReader( tempPath, true ))
+		{
+			assertArrayEquals( storedCov3d,
+					GeffUtils.readAsDoubleMatrix( reader, "/nodes/props/covariance3d/values", "covariance3d" ).rowAt( 0 ),
+					1e-9, "covariance3d is not stored in the z, y, x axes order" );
+			assertArrayEquals( storedCov2d,
+					GeffUtils.readAsDoubleMatrix( reader, "/nodes/props/covariance2d/values", "covariance2d" ).rowAt( 0 ),
+					1e-9, "covariance2d is not stored in the y, x axes order" );
+		}
+
+		// Reading permutes back into the x, y, z order
+		final List< GeffNode > readNodes = GeffNode.readFromZarr( tempPath, GeffMetadata.readFromZarr( tempPath ) );
+		assertEquals( 1, readNodes.size() );
+		assertArrayEquals( cov3d, readNodes.get( 0 ).getCovariance3d(), 1e-9 );
+		assertArrayEquals( cov2d, readNodes.get( 0 ).getCovariance2d(), 1e-9 );
+	}
+
+	@Test
+	@DisplayName( "Test reading covariances written with the legacy [N, 4] and [N, 6] shapes" )
+	void testLegacyCovarianceShapes( @TempDir Path tempDir ) throws IOException
+	{
+		final String tempPath = tempDir.toString() + "/test-legacy-cov.zarr/tracks";
+
+		final List< GeffNode > nodes = new ArrayList<>();
+		for ( int i = 0; i < 3; i++ )
+		{
+			final GeffNode node = new GeffNode();
+			node.setId( i );
+			node.setT( i );
+			node.setX( i * 1.0 );
+			node.setY( i * 2.0 );
+			node.setZ( i * 3.0 );
+			nodes.add( node );
+		}
+
+		final GeffMetadata metadata = new GeffMetadata( Geff.VERSION, true );
+		GeffNode.writeToZarr( nodes, tempPath, metadata );
+		GeffMetadata.writeToZarr( metadata, tempPath );
+
+		// Overwrite the covariances with the flat shapes written by earlier
+		// versions of geff-java: [N, 4] for covariance2d, and the upper
+		// triangle { xx, xy, xz, yy, yz, zz } as [N, 6] for covariance3d
+		final double[][] legacyCov2d = {
+				{ 2.0, 0.5, 0.5, 3.0 },
+				{ 1.5, -0.3, -0.3, 1.8 },
+				{ 4.0, 0.0, 0.0, 4.0 },
+		};
+		final double[][] legacyCov3d = {
+				{ 2.0, 0.5, 0.1, 3.0, 0.2, 1.5 },
+				{ 1.0, 0.0, 0.0, 1.0, 0.0, 1.0 },
+				{ 5.0, -0.2, 0.3, 4.0, -0.1, 3.0 },
+		};
+		final List< Integer > indices = new ArrayList<>();
+		for ( int i = 0; i < nodes.size(); i++ )
+			indices.add( i );
+		try (final N5ZarrWriter writer = new N5ZarrWriter( tempPath, true ))
+		{
+			writer.remove( "/nodes/props/covariance2d/values" );
+			writer.remove( "/nodes/props/covariance3d/values" );
+			GeffUtils.writeDoubleMatrix( indices, 4, i -> legacyCov2d[ i ], writer, "/nodes/props/covariance2d/values", 3 );
+			GeffUtils.writeDoubleMatrix( indices, 6, i -> legacyCov3d[ i ], writer, "/nodes/props/covariance3d/values", 3 );
+		}
+
+		final List< GeffNode > readNodes = GeffNode.readFromZarr( tempPath, GeffMetadata.readFromZarr( tempPath ) );
+		assertEquals( nodes.size(), readNodes.size() );
+		for ( int i = 0; i < nodes.size(); i++ )
+		{
+			// covariance2d is unchanged: [N, 4] and [N, 2, 2] hold the same
+			// 4 values per node
+			assertArrayEquals( legacyCov2d[ i ], readNodes.get( i ).getCovariance2d(), 1e-9,
+					"covariance2d mismatch at node " + i );
+
+			// covariance3d is expanded from the upper triangle to the full
+			// symmetric 3x3 matrix
+			final double[] t = legacyCov3d[ i ];
+			final double[] expected = { t[ 0 ], t[ 1 ], t[ 2 ], t[ 1 ], t[ 3 ], t[ 4 ], t[ 2 ], t[ 4 ], t[ 5 ] };
+			assertArrayEquals( expected, readNodes.get( i ).getCovariance3d(), 1e-9,
+					"covariance3d mismatch at node " + i );
+		}
 	}
 
 	@Test
@@ -359,9 +508,9 @@ void testCovariance3dRoundTrip( @TempDir Path tempDir ) throws IOException
 		final String tempPath = tempDir.toString() + "/test-cov3d.zarr/tracks";
 
 		final double[][] cov3dValues = {
-				{ 2.0, 0.5, 0.1, 3.0, 0.2, 1.5 },
-				{ 1.0, 0.0, 0.0, 1.0, 0.0, 1.0 },
-				{ 5.0, -0.2, 0.3, 4.0, -0.1, 3.0 },
+				{ 2.0, 0.5, 0.1, 0.5, 3.0, 0.2, 0.1, 0.2, 1.5 },
+				{ 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0 },
+				{ 5.0, -0.2, 0.3, -0.2, 4.0, -0.1, 0.3, -0.1, 3.0 },
 		};
 
 		final List< GeffNode > nodes = new ArrayList<>();
@@ -388,6 +537,10 @@ void testCovariance3dRoundTrip( @TempDir Path tempDir ) throws IOException
 			assertArrayEquals( cov3dValues[ i ], readNodes.get( i ).getCovariance3d(), 1e-9,
 					"covariance3d mismatch at node " + i );
 		}
+
+		// The GEFF spec defines covariance3d as [N, 3, 3]
+		assertArrayEquals( new long[] { cov3dValues.length, 3, 3 },
+				zarrShape( tempPath + "/nodes/props/covariance3d/values" ) );
 	}
 
 	@Test
@@ -496,4 +649,91 @@ void testCustomAxisNames( @TempDir Path tempDir ) throws IOException
 			assertEquals( expected.getSegmentId(), actual.getSegmentId(), "segmentId mismatch at node " + i );
 		}
 	}
+
+	@Test
+	@DisplayName( "Test reading/writing radius and covariance3d under custom prop identifiers" )
+	void testCustomPropIdentifiers( @TempDir Path tempDir ) throws IOException
+	{
+		// Declare radius, covariance3d and distance under identifiers that
+		// differ from their standard names.
+		final GeffMetadata metadata = new GeffMetadata( Geff.VERSION, true );
+		metadata.setGeffAxes( new GeffAxis[] {
+				new GeffAxis( "t", GeffAxis.TYPE_TIME, GeffAxis.UNIT_SECOND, 0.0, 100.0 ),
+				new GeffAxis( "y", GeffAxis.TYPE_SPACE, GeffAxis.UNIT_PIXEL, 0.0, 768.0 ),
+				new GeffAxis( "x", GeffAxis.TYPE_SPACE, GeffAxis.UNIT_PIXEL, 0.0, 1024.0 )
+		} );
+
+		final Map< String, PropMetadata > nodeProps = new HashMap<>();
+		nodeProps.put( "t", new PropMetadata( "t", "int32", false, null, null, null ) );
+		nodeProps.put( "x", new PropMetadata( "x", "float64", false, null, null, null ) );
+		nodeProps.put( "y", new PropMetadata( "y", "float64", false, null, null, null ) );
+		nodeProps.put( "radius", new PropMetadata( "cell_radius", "float64", false, null, null, null ) );
+		nodeProps.put( "covariance3d", new PropMetadata( "cell_covariance3d", "float64", false, null, null, null ) );
+		metadata.setNodePropsMetadata( nodeProps );
+
+		final Map< String, PropMetadata > edgeProps = new HashMap<>();
+		edgeProps.put( "distance", new PropMetadata( "edge_distance", "float64", false, null, null, null ) );
+		metadata.setEdgePropsMetadata( edgeProps );
+
+		final List< GeffNode > nodes = new ArrayList<>();
+		for ( int i = 0; i < 3; i++ )
+		{
+			final GeffNode node = new GeffNode();
+			node.setId( i );
+			node.setT( i );
+			node.setX( i * 10.0 );
+			node.setY( i * 20.0 );
+			node.setRadius( 2.0 + i );
+			node.setCovariance3d( new double[] { i, 0.0, 0.0, 0.0, i + 1.0, 0.0, 0.0, 0.0, i + 2.0 } );
+			nodes.add( node );
+		}
+
+		final List< GeffEdge > edges = new ArrayList<>();
+		for ( int i = 0; i < 2; i++ )
+			edges.add( new GeffEdge( i, i, i + 1, 0.5, 3.0 + i ) );
+
+		final String tempPath = tempDir.toString() + "/test-custom-prop-ids.zarr/tracks";
+		GeffNode.writeToZarr( nodes, tempPath, metadata );
+		GeffEdge.writeToZarr( edges, tempPath, metadata );
+		GeffMetadata.writeToZarr( metadata, tempPath );
+
+		// The values are stored under the declared identifiers, not under the
+		// standard names
+		assertTrue( Files.isDirectory( Paths.get( tempPath, "nodes", "props", "cell_radius", "values" ) ) );
+		assertTrue( Files.isDirectory( Paths.get( tempPath, "nodes", "props", "cell_covariance3d", "values" ) ) );
+		assertTrue( Files.isDirectory( Paths.get( tempPath, "edges", "props", "edge_distance", "values" ) ) );
+		assertFalse( Files.exists( Paths.get( tempPath, "nodes", "props", "radius" ) ) );
+		assertFalse( Files.exists( Paths.get( tempPath, "nodes", "props", "covariance3d" ) ) );
+		assertFalse( Files.exists( Paths.get( tempPath, "edges", "props", "distance" ) ) );
+
+		// Reading resolves the identifiers from node_props_metadata and
+		// edge_props_metadata again
+		final GeffMetadata readMetadata = GeffMetadata.readFromZarr( tempPath );
+		assertEquals( "cell_radius", readMetadata.getNodePropIdentifier( "radius" ) );
+		assertEquals( "cell_covariance3d", readMetadata.getNodePropIdentifier( "covariance3d" ) );
+		assertEquals( "edge_distance", readMetadata.getEdgePropIdentifier( "distance" ) );
+		// Props that are not declared keep their standard name
+		assertEquals( "covariance2d", readMetadata.getNodePropIdentifier( "covariance2d" ) );
+
+		final List< GeffNode > readNodes = GeffNode.readFromZarr( tempPath, readMetadata );
+		assertEquals( nodes.size(), readNodes.size() );
+		for ( int i = 0; i < nodes.size(); i++ )
+		{
+			final GeffNode expected = nodes.get( i );
+			final GeffNode actual = readNodes.get( i );
+			assertEquals( expected.getRadius(), actual.getRadius(), 1e-9, "radius mismatch at node " + i );
+			assertArrayEquals( expected.getCovariance3d(), actual.getCovariance3d(), 1e-9,
+					"covariance3d mismatch at node " + i );
+			// A prop read under its own identifier is not also reported as a
+			// custom prop
+			assertNull( actual.getProp( "radius" ) );
+			assertNull( actual.getProp( "cell_radius" ) );
+		}
+
+		final List< GeffEdge > readEdges = GeffEdge.readFromZarr( tempPath, readMetadata );
+		assertEquals( edges.size(), readEdges.size() );
+		for ( int i = 0; i < edges.size(); i++ )
+			assertEquals( edges.get( i ).getDistance(), readEdges.get( i ).getDistance(), 1e-9,
+					"distance mismatch at edge " + i );
+	}
 }