docs: def read_geom_h5(xdmf_file, snapshot):
      """Extract geometry information from hdf5 files.
  
      Args:
          xdmf_file (:class:`pathlib.Path`): path of the xdmf file.
          snapshot (int): snapshot number.
      Returns:
          (dict, root): geometry information and root of xdmf document.
      """
      header = {}
      xdmf_root = xmlET.parse(str(xdmf_file)).getroot()
      if snapshot is None:
          return None, xdmf_root
  
      # Domain, Temporal Collection, Snapshot
      # should check that this is indeed the required snapshot
      elt_snap = xdmf_root[0][0][snapshot]
      header['ti_ad'] = float(elt_snap.find('Time').get('Value'))
      header['mo_lambda'] = _maybe_get(elt_snap, 'mo_lambda', 'Value', float)
      header['mo_thick_sol'] = _maybe_get(elt_snap, 'mo_thick_sol', 'Value',
                                          float)
      header['ntb'] = 1
      coord_h5 = []  # all the coordinate files
      coord_shape = []  # shape of meshes
      twod = None
      for elt_subdomain in elt_snap.findall('Grid'):
          if elt_subdomain.get('Name').startswith('meshYang'):
              header['ntb'] = 2
              break  # iterate only through meshYin
          elt_geom = elt_subdomain.find('Geometry')
          if elt_geom.get('Type') == 'X_Y' and twod is None:
              twod = ''
              for data_item in elt_geom.findall('DataItem'):
                  coord = data_item.text.strip()[-1]
                  if coord in 'XYZ':
                      twod += coord
          data_item = elt_geom.find('DataItem')
          coord_shape.append(_get_dim(data_item))
          coord_h5.append(
              xdmf_file.parent / data_item.text.strip().split(':/', 1)[0])
      _read_coord_h5(coord_h5, coord_shape, header, twod)
      return header, xdmf_root