"""
This is a recipe to show how to open a dataset and extract it to a file at a
fixed resolution with no interpolation or smoothing.  Additionally, this
recipe shows how to insert a dataset into an external HDF5 file using h5py.
For more information see :class:`covering_grid`.
"""
from yt.mods import *

# For this example we will use h5py to write to our output file.
import h5py

fn = "RedshiftOutput0005" # parameter file to load
pf = load(fn) # load data

# This is the resolution we will extract at
DIMS = 128

# Now, we construct an object that describes the data region and structure we
# want
cube = pf.h.covering_grid(2, # The level we are willing to extract to; higher
                             # levels than this will not contribute to the data!
                          # Now we set our spatial extent...
                          left_edge=[0.0, 0.0, 0.0], 
                          right_edge=[1.0, 1.0, 1.0],
                          # How many dimensions along each axis
                          dims=[DIMS,DIMS,DIMS],
                          # And any fields to preload (this is optional!)
                          fields=["Density"])

# Now we open our output file using h5py
# Note that we open with 'w' which will overwrite existing files!
f = h5py.File("my_data.h5", "w") 

# We create a dataset at the root note, calling it density...
f.create_dataset("/density", data=cube["Density"])

# We close our file
f.close()