docs: def _read_coord_h5(files, shapes, header, twod):
      """Read all coord hdf5 files of a snapshot.
  
      Args:
          files (list of pathlib.Path): list of NodeCoordinates files of
              a snapshot.
          shapes (list of (int,int)): shape of mesh grids.
          header (dict): geometry info.
          twod (str): 'XZ', 'YZ' or None depending on what is relevant.
      """
      meshes = []
      for h5file, shape in zip(files, shapes):
          meshes.append({})
          with h5py.File(h5file, 'r') as h5f:
              for coord, mesh in h5f.items():
                  # for some reason, the array is transposed!
                  meshes[-1][coord] = mesh[()].reshape(shape).T
                  meshes[-1][coord] = _make_3d(meshes[-1][coord], twod)
  
      header['ncs'] = _ncores(meshes, twod)
      header['nts'] = list((meshes[0]['X'].shape[i] - 1) * header['ncs'][i]
                           for i in range(3))
      header['nts'] = np.array([max(1, val) for val in header['nts']])
      # meshes could also be defined in legacy parser, so that these can be used
      # in geometry setup
      meshes = _conglomerate_meshes(meshes, header)
      if np.any(meshes['Z'][:, :, 0] != 0):
          # spherical
          header['x_mesh'] = np.copy(meshes['Y'])  # annulus geometry...
          header['y_mesh'] = np.copy(meshes['Z'])
          header['z_mesh'] = np.copy(meshes['X'])
          header['r_mesh'] = np.sqrt(header['x_mesh']**2 + header['y_mesh']**2 +
                                     header['z_mesh']**2)
          header['t_mesh'] = np.arccos(header['z_mesh'] / header['r_mesh'])
          header['p_mesh'] = np.roll(
              np.arctan2(header['y_mesh'], -header['x_mesh']) + np.pi, -1, 1)
          header['e1_coord'] = header['t_mesh'][:, 0, 0]
          header['e2_coord'] = header['p_mesh'][0, :, 0]
          header['e3_coord'] = header['r_mesh'][0, 0, :]
      else:
          header['e1_coord'] = meshes['X'][:, 0, 0]
          header['e2_coord'] = meshes['Y'][0, :, 0]
          header['e3_coord'] = meshes['Z'][0, 0, :]
      header['aspect'] = (header['e1_coord'][-1] - header['e2_coord'][0],
                          header['e1_coord'][-1] - header['e2_coord'][0])
      header['rcmb'] = header['e3_coord'][0]
      if header['rcmb'] == 0:
          header['rcmb'] = -1
      else:
          # could make the difference between r_coord and z_coord
          header['e3_coord'] = header['e3_coord'] - header['rcmb']
      if twod is None or 'X' in twod:
          header['e1_coord'] = header['e1_coord'][:-1]
      if twod is None or 'Y' in twod:
          header['e2_coord'] = header['e2_coord'][:-1]
      header['e3_coord'] = header['e3_coord'][:-1]