docs: def read_data(header, fh=None, filename=None, index_order='F'):
      """Read data from file into :class:`numpy.ndarray`
  
      The two parameters :obj:`fh` and :obj:`filename` are optional depending on the parameters but it never hurts to
      specify both. The file handle (:obj:`fh`) is necessary if the header is attached with the NRRD data. However, if
      the NRRD data is detached from the header, then the :obj:`filename` parameter is required to obtain the absolute
      path to the data file.
  
      See :ref:`user-guide:Reading NRRD files` for more information on reading NRRD files.
  
      Parameters
      ----------
      header : :class:`dict` (:class:`str`, :obj:`Object`)
          Parsed fields/values obtained from :meth:`read_header` function
      fh : file-object, optional
          File object pointing to first byte of data. Only necessary if data is attached to header.
      filename : :class:`str`, optional
          Filename of the header file. Only necessary if data is detached from the header. This is used to get the
          absolute data path.
      index_order : {'C', 'F'}, optional
          Specifies the index order of the resulting data array. Either 'C' (C-order) where the dimensions are ordered from
          slowest-varying to fastest-varying (e.g. (z, y, x)), or 'F' (Fortran-order) where the dimensions are ordered
          from fastest-varying to slowest-varying (e.g. (x, y, z)).
  
      Returns
      -------
      data : :class:`numpy.ndarray`
          Data read from NRRD file
  
      See Also
      --------
      :meth:`read`, :meth:`read_header`
      """
  
      if index_order not in ['F', 'C']:
          raise NRRDError('Invalid index order')
  
      # Check that the required fields are in the header
      for field in _NRRD_REQUIRED_FIELDS:
          if field not in header:
              raise NRRDError('Header is missing required field: "%s".' % field)
  
      if header['dimension'] != len(header['sizes']):
          raise NRRDError('Number of elements in sizes does not match dimension. Dimension: %i, len(sizes): %i' % (
              header['dimension'], len(header['sizes'])))
  
      # Determine the data type from the header
      dtype = _determine_datatype(header)
  
      # Determine the byte skip, line skip and the data file
      # These all can be written with or without the space according to the NRRD spec, so we check them both
      line_skip = header.get('lineskip', header.get('line skip', 0))
      byte_skip = header.get('byteskip', header.get('byte skip', 0))
      data_filename = header.get('datafile', header.get('data file', None))
  
      # If the data file is separate from the header file, then open the data file to read from that instead
      if data_filename is not None:
          # If the pathname is relative, then append the current directory from the filename
          if not os.path.isabs(data_filename):
              if filename is None:
                  raise NRRDError('Filename parameter must be specified when a relative data file path is given')
  
              data_filename = os.path.join(os.path.dirname(filename), data_filename)
  
          # Override the fh parameter with the data filename
          # Note that this is opened without a "with" block, thus it must be closed manually in all circumstances
          fh = open(data_filename, 'rb')
  
      # Get the total number of data points by multiplying the size of each dimension together
      total_data_points = header['sizes'].prod()
  
      # Skip the number of lines requested when line_skip >= 0
      # Irrespective of the NRRD file having attached/detached header
      # Lines are skipped before getting to the beginning of the data
      if line_skip >= 0:
          for _ in range(line_skip):
              fh.readline()
      else:
          # Must close the file because if the file was opened above from detached filename, there is no "with" block to
          # close it for us
          fh.close()
  
          raise NRRDError('Invalid lineskip, allowed values are greater than or equal to 0')
  
      # Skip the requested number of bytes or seek backward, and then parse the data using NumPy
      if byte_skip < -1:
          # Must close the file because if the file was opened above from detached filename, there is no "with" block to
          # close it for us
          fh.close()
  
          raise NRRDError('Invalid byteskip, allowed values are greater than or equal to -1')
      elif byte_skip >= 0:
          fh.seek(byte_skip, os.SEEK_CUR)
      elif byte_skip == -1 and header['encoding'] not in ['gzip', 'gz', 'bzip2', 'bz2']:
          fh.seek(-dtype.itemsize * total_data_points, os.SEEK_END)
      else:
          # The only case left should be: byte_skip == -1 and header['encoding'] == 'gzip'
          byte_skip = -dtype.itemsize * total_data_points
  
      # If a compression encoding is used, then byte skip AFTER decompressing
      if header['encoding'] == 'raw':
          data = np.fromfile(fh, dtype)
      elif header['encoding'] in ['ASCII', 'ascii', 'text', 'txt']:
          data = np.fromfile(fh, dtype, sep=' ')
      else:
          # Handle compressed data now
          # Construct the decompression object based on encoding
          if header['encoding'] in ['gzip', 'gz']:
              decompobj = zlib.decompressobj(zlib.MAX_WBITS | 16)
          elif header['encoding'] in ['bzip2', 'bz2']:
              decompobj = bz2.BZ2Decompressor()
          else:
              # Must close the file because if the file was opened above from detached filename, there is no "with" block
              # to close it for us
              fh.close()
  
              raise NRRDError('Unsupported encoding: "%s"' % header['encoding'])
  
          # Loop through the file and read a chunk at a time (see _READ_CHUNKSIZE why it is read in chunks)
          decompressed_data = bytearray()
  
          # Read all of the remaining data from the file
          # Obtain the length of the compressed data since we will be using it repeatedly, more efficient
          compressed_data = fh.read()
          compressed_data_len = len(compressed_data)
          start_index = 0
  
          # Loop through data and decompress it chunk by chunk
          while start_index < compressed_data_len:
              # Calculate the end index = start index plus chunk size
              # Set to the string length to read the remaining chunk at the end
              end_index = min(start_index + _READ_CHUNKSIZE, compressed_data_len)
  
              # Decompress and append data
              decompressed_data += decompobj.decompress(compressed_data[start_index:end_index])
  
              # Update start index
              start_index = end_index
  
          # Delete the compressed data since we do not need it anymore
          # This could potentially be using a lot of memory
          del compressed_data
  
          # Byte skip is applied AFTER the decompression. Skip first x bytes of the decompressed data and parse it using
          # NumPy
          data = np.frombuffer(decompressed_data[byte_skip:], dtype)
  
      # Close the file, even if opened using "with" block, closing it manually does not hurt
      fh.close()
  
      if total_data_points != data.size:
          raise NRRDError('Size of the data does not equal the product of all the dimensions: {0}-{1}={2}'
                          .format(total_data_points, data.size, total_data_points - data.size))
  
      # In the NRRD header, the fields are specified in Fortran order, i.e, the first index is the one that changes
      # fastest and last index changes slowest. This needs to be taken into consideration since numpy uses C-order
      # indexing.
      
      # The array shape from NRRD (x,y,z) needs to be reversed as numpy expects (z,y,x).
      data = np.reshape(data, tuple(header['sizes'][::-1]))
  
      # Transpose data to enable Fortran indexing if requested.
      if index_order == 'F':
          data = data.T
  
      return data