02-regions.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lesson 02: Exploring data with yt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-info\">\n",
    " \n",
    "## Goals:\n",
    "\n",
    "* Exploring fields (min/max), the quantities interface\n",
    "* Selecting regions of data\n",
    "* slices and rays\n",
    "* visualizing our selections\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we've learned a bit about how to load data with yt and what available objects there are in a field, we can use yt to perform analysis on those fields. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exploring fields"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, let's explore some values that occur in our dat. For this tutorial we'll be using `enzo tiny cosmology`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import yt\n",
    "ds = yt.load_sample('enzo_tiny_cosmology')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This particular dataset is a *time series*. yt automatically loads the last file in a time series dataset; in this case it is `DD0046`.\n",
    "\n",
    "Now let's explore a bit of the data contained in this file. First, let's see what fields we can explore:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds.field_list"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Great! That's a lot of fields. You might know that yt has a lot of tools in visualization, but it can also be used to inspect data. Here's an example where we find the extrema (the min and max) values of the \"density\" field. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dd = ds.all_data()\n",
    "dd.quantities.extrema(\"density\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ok, a few things have happened here. First, we've done this operation `ds.all_data()` which returns a *region* in yt based on the domain boundaries detected when the data was read in. We've returned all of the data in this dataset. \n",
    "\n",
    "We can check this by inspecting the `dd` object, which is a `YTRegion`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Region objects (and in fact any selection of the data) have **quantites** that we can use on them, but they need to be a type data selection. In the previous example, we accessed the `extrema` quantity, whcih returns a unyt array of length two of the minimum and maximum values in the density field of this dataset. \n",
    "\n",
    "Let's see what other quantities are available: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print (list(dd.quantities.keys()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "ok, now let's examine what arguments we need to pass to get something we want. How about `WeightedAverageQuantity`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dd.quantities.weighted_average_quantity?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and now let's calculate the temperature weighted density of this dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dd.quantities.weighted_average_quantity(\"density\", weight=\"temperature\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "we can also pass a list of fields with which to calculate a weighted average quantity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dd.quantities.weighted_average_quantity([\"density\", \"temperature\"], weight=\"temperature\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "\n",
    "## Interactive Exercise 01\n",
    "    \n",
    "Load in the dataset `enzo_cosmology_plus` and find the  following quantities:\n",
    "    \n",
    "* min and maximum values of each `enzo` field in the entire dataset\n",
    "* the total gas mass in the dataset\n",
    "* the `cell_volume` weighted average density\n",
    "* the location of the min and max values \n",
    "    \n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This is starting cell to do exercise 01"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">\n",
    " \n",
    "### Tips:\n",
    "    \n",
    "* instead of doing `dd = ds.all_data()` and then selecting data with dens = dd[\"density\"], yt has the option to do an automatioc `region selector` with `dens = ds.r[\"density\"]`. `ds.r()` is a `RegionExpression` object and doesn't have any quantities associated with it. When used with a field argument, this function returns a flattened array of the data. \n",
    "* for this tutorial we've only loaded in the last file `DD0046`. However, because this is a time-series dataset there are a number of other files. yt can recognize this by loading with wildcard operators. e.g. `ts = yt.load(\"enzo_tiny_cosmology/DD????/DD????\")`. \n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Volumetric Region Selectors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So far we've loaded in the entire dataset with `all_data()`. We learned that this is a special function that returns a YTRegion that spans the entire domain of the data. However, yt has other data selectors that are available, like `sphere`. Let's do a region selection of the entire dataset manually and check that the values are the same as what we saw in the first section. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bx = ds.box([0.,0.,0.], [1.,1.,1.]) \n",
    "# we could also specify a center and use ds.region(). ds.box() assumes the center of this prism is the \n",
    "# centerpoint of the left and right edges. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bx"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "bx is a YTRegion object, which is what would be returned had we used the `ds.region()` selector as well. Let's compare it to the dd object we used before. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dd "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "They look the same at first pass! Let's see what happens if we look at the min, the max, and a random weighted field! "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(bx.quantities.extrema('temperature'))\n",
    "print(dd.quantities.extrema('temperature'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(bx.quantities.weighted_average_quantity('temperature', weight='ones'))\n",
    "print(dd.quantities.weighted_average_quantity('temperature', weight='ones'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ok! So now we can be reasonably sure that the way we selected data with `ds.all_data()` is a nice shorthand to select all data in our dataset. However, we can create the same object with `box()`\n",
    "\n",
    "Now let's do a region selector with a sphere. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds.sphere?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ok, so let's define a centerpoint at 0.5, 0.5, 0.5 in code units and extend the sphere outwards to 0.25 code units. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "center = [0.5, 0.5, 0.5]\n",
    "sp = ds.sphere(center , 0.25) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Is the sphere a YTRegion object like the selector we saw before? Nope! It's a YTSphere. Both are subclasses of the same type of yt selector object and so have similar operations available to them though. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's take a look at the extrema and compare them to what we saw in the larger dataset! "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(sp.quantities.extrema('density'))\n",
    "print(dd.quantities.extrema('density'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(sp.quantities.weighted_average_quantity('density', weight='cell_volume'))\n",
    "print(dd.quantities.weighted_average_quantity('density', weight='cell_volume'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(sp.mean('density', weight='cell_volume'))\n",
    "print(dd.mean('density', weight='cell_volume'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "\n",
    "## Interactive Exercise 02\n",
    "    \n",
    "With your previously loaded dataset `enzo_cosmology_plus` and select a region (from the available options [in the docs](https://yt-project.org/doc/analyzing/objects.html#region-reference)) in the domain center with random dimensions. \n",
    "    \n",
    "Now compare the total gas mass, the min, and the max values of the density field for the region you've selected to that of the total dataset that you did in exercise 01. Do they differ? Do they look the same? \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This is a starting cell to do exercise 02 "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">\n",
    " \n",
    "### Tips:\n",
    "    \n",
    "* `obj.mean()`, `.min()`, and `.max()` are easy accessors for `.quantiy.weighted_average_quantity()`, and `.quantity.extrema()`\n",
    "* The way this tutorial has specified the left, right, and center of the object is in code units, which generally span from 0 to 1. However, you can use any unit quantity you'd like thanks to yt's unit interface. Maybe you'd like to use MPc instead? No problem! Try `sp = ds.sphere(center, (10, 'Mpc'))`\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Rays and Slices "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So far we've seen some volumetric selector objects and how we can calculate different things on those object with the quantities interface. These are very useful to get intuitions into our data! We can find out a lot about what's going on in our fields by using these selectors. \n",
    "\n",
    "However, yt has other selector objects that can give us insight into our data. Let's start with `ray` objects. We can define a path through our data with a starting point and an ending point. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds.ray?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ra = ds.ray([0.1, 0.2, 0.3], [0.9, 0.8, 0.7])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ra"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can plot this with matplotlib if we'd like. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t = ra['t']\n",
    "dens = ra['density']\n",
    "temp = ra['temperature']\n",
    "fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, figsize=(12, 4))\n",
    "ax1.plot(t, dens)\n",
    "ax1.set_yscale('log')\n",
    "ax1.set(xlabel='distance along ray', ylabel='density (g/cm**3)',\n",
    "       title='hit me with that laser beam!')\n",
    "\n",
    "ax2.plot(t, temp)\n",
    "ax2.set_yscale('log')\n",
    "ax2.set(xlabel='distance along ray', ylabel='temperature (K)',\n",
    "       title='hit me with that  thermometer!')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also use the same `quantity` accessors that we used before! "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(ra.quantities.extrema('density'))\n",
    "print(ra.quantities.extrema('temperature'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Rays are a good way to gain intuition to what our data looks like along a specific path. This might be hard to immediately glean by inspecting values or visualizing the whole dataset. \n",
    "\n",
    "We can also create 2d data objects, such as a `slice`! "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds.slice?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sl1 = ds.slice('z', 0.4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sl1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(sl1.min('density'))\n",
    "print(sl1.max('density'))\n",
    "print(sl1.mean('density'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also do things like [off-axis slices](https://yt-project.org/docs/dev/quickstart/data_objects_and_time_series.html#Off-Axis-Slices)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds.cutting?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sl2 = ds.cutting([0.2, 0.3, 0.5], \"min\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sl2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(sl2.min('density'))\n",
    "print(sl2.max('density'))\n",
    "print(sl2.mean('density'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "\n",
    "## Interactive Exercise 03\n",
    "    \n",
    "* Find the locations of the maximum and minumum values in the temperature field of `enzo cosmology plus` and create a ray with it. How do the quantity values of the ray compare to the shape you created in exercise 02? \n",
    "* Create a slice object at an arbitrary location with this dataset. \n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This is a cell to start for exercise 03 "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">\n",
    " \n",
    "### Tips:\n",
    "    \n",
    "* `ds.r()` gives us shorthand notation to slice our data without using the `slice()` method. For example `ds.r[:,:, 0.5]` will return a YTSlice object at the z midplane of the dataset. \n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualizing our Selections and Advanced Selections"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `.plot()` method is an easily accessible way for us to visualize the selections we've been doing with our data. Let's try them out on the objects we've already created in the previous selection. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p1 = sl1.plot('density')\n",
    "p1.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p2 = sl1.plot('temperature')\n",
    "p2.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also use volume selectors in slices to do some more advanced data selection. Let's use a sphere and a slice together and see what that looks like! "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sp2 = ds.sphere(center, 0.4) \n",
    "sl3 = ds.slice('z', 0.4, data_source=sp2)\n",
    "p3 = sl3.plot('density')\n",
    "p3.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can even chain together volumes with an union and plot that! "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sp3 = ds.sphere([0.4, 0.5, 0.5], 0.18)\n",
    "sp4 = ds.sphere([0.7, 0.65, 0.5], 0.28)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "isp = ds.union( [sp3, sp4] )\n",
    "sl4 = ds.slice('z', 0.5, data_source=isp)\n",
    "p4 = sl4.plot('temperature')\n",
    "p4.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "\n",
    "## Interactive Exercise 04\n",
    "    \n",
    "* create a a slice object object in the `temperature` field and visualize it with the .plot() object. \n",
    "* bonus: create an intersection object with a slice object and the volumetric object you made in exercise 03 and visualize that! \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This is a cell to start exercise 04"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">\n",
    " \n",
    "### Tips:\n",
    "    \n",
    "* .plot() works yt slice and projection objects\n",
    "* you can visualize a volumetric object by passing it through the `data_source` arg of a slice. However, you need to be careful that the regions actually intersect, or else your returned object will be all zeroes! \n",
    "* Another type of selection that wasn't covered in this tutorial is a [Profile](https://yt-project.org/docs/dev/quickstart/derived_fields_and_profiles.html) object. If you have time, try to create one! \n",
    "* Another type of object we didn't create was a projection! Try it out. \n",
    "* more plotting will happen in the next lesson! \n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-danger\">\n",
    "\n",
    "## Bringing It All Together: Challenge Exercise\n",
    "    \n",
    "With `enzo_cosmology_plus`, create an intersection object composed of a disk and a sphere at arbitrary locations in the data (make sure they overlap). Then make a slice of this object at its midpoint and visualize it with `.plot()`. Calculate the total mass of each object you created. \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This is a starting cell to do the challenge exercise"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Takeaways"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">\n",
    "\n",
    "### There are many ways we can do the same operation with yt\n",
    "### Derived quantities can be calculated on YT selection objects. These objects can be:\n",
    "    \n",
    "* [geometric](https://yt-project.org/doc/analyzing/objects.html#geometric-objects) (object is based on geometry)\n",
    "* [filtering objects](https://yt-project.org/doc/analyzing/objects.html#filtering-and-collection-objects) (object is based on field criteria)\n",
    "* [construction objects](https://yt-project.org/doc/analyzing/objects.html#construction-objects) (data is based on additional analysis)\n",
    "* [collections](https://yt-project.org/doc/analyzing/objects.html#filtering-and-collection-objects) (object is a collection of other objects).\n",
    "    \n",
    "### Slice objects can be visualized with the `.plot()` method \n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}