def read_field_h5(xdmf_file, fieldname, snapshot, header=None):
    """Extract field data from hdf5 files.

    Args:
        xdmf_file (:class:`pathlib.Path`): path of the xdmf file.
        fieldname (str): name of field to extract.
        snapshot (int): snapshot number.
        header (dict): geometry information.
    Returns:
        (dict, numpy.array): geometry information and field data. None
            is returned if data is unavailable.
    """
    if header is None:
        header, xdmf_root = read_geom_h5(xdmf_file, snapshot)
    else:
        xdmf_root = xmlET.parse(str(xdmf_file)).getroot()

    npc = header['nts'] // header['ncs']  # number of grid point per node
    flds = np.zeros(_flds_shape(fieldname, header))
    data_found = False

    for elt_subdomain in xdmf_root[0][0][snapshot].findall('Grid'):
        ibk = int(elt_subdomain.get('Name').startswith('meshYang'))
        for data_attr in elt_subdomain.findall('Attribute'):
            if data_attr.get('Name') != fieldname:
                continue
            icore, fld = _get_field(xdmf_file, data_attr.find('DataItem'))
            # for some reason, the field is transposed
            fld = fld.T
            shp = fld.shape
            if shp[-1] == 1 and header['nts'][0] == 1:  # YZ
                fld = fld.reshape((shp[0], 1, shp[1], shp[2]))
                if header['rcmb'] < 0:
                    fld = fld[(2, 0, 1), ...]
            elif shp[-1] == 1:  # XZ
                fld = fld.reshape((shp[0], shp[1], 1, shp[2]))
                if header['rcmb'] < 0:
                    fld = fld[(0, 2, 1), ...]
            elif header['nts'][1] == 1:  # cart XZ
                fld = fld.reshape((1, shp[0], 1, shp[1]))
            ifs = [icore // np.prod(header['ncs'][:i]) % header['ncs'][i] *
                   npc[i] for i in range(3)]
            if header['zp']:  # remove top row
                fld = fld[:, :, :, :-1]
            flds[:,
                 ifs[0]:ifs[0] + npc[0] + header['xp'],
                 ifs[1]:ifs[1] + npc[1] + header['yp'],
                 ifs[2]:ifs[2] + npc[2],
                 ibk] = fld
            data_found = True

    flds = _post_read_flds(flds, header)

    return (header, flds) if data_found else None
Extract field data from hdf5 files.

    Args:
        xdmf_file (:class:`pathlib.Path`): path of the xdmf file.
        fieldname (str): name of field to extract.
        snapshot (int): snapshot number.
        header (dict): geometry information.
    Returns:
        (dict, numpy.array): geometry information and field data. None
            is returned if data is unavailable.