3d_print_the_universe.ipynb

{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import h5py\n",
      "import numpy as np\n",
      "import yt\n",
      "from scipy.ndimage.filters import gaussian_filter\n",
      "from scipy.misc import imsave"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# load volumetric density data\n",
      "# we used hdf5, but any 3d numpy array works\n",
      "f = h5py.File('dset.hdf5','r')\n",
      "rhof = f.get(\"RHO\")\n",
      "rho = np.array(rhof)\n",
      "f.close()\n",
      "\n",
      "# log-scale the data, and create a smoothed version for \n",
      "# better-looking 3d models\n",
      "rholog = np.log10(rho)\n",
      "rhosmooth = gaussian_filter(rholog, 1.0)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# make a log-scaled heightmap of projected density,\n",
      "# scaled to the range [0, 1]\n",
      "hmap = np.log10(np.sum(rho[40:90,:,:], axis=0))\n",
      "hmin = np.min(hmap)\n",
      "hmax = np.max(hmap)\n",
      "hmap = (hmap - hmin)/(hmax - hmin)\n",
      "\n",
      "# Save as a grayscale image \n",
      "imsave('hmap.png', (255.0*hmap).astype(np.uint8))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# load data into yt\n",
      "data = dict(density = (rhosmooth, \"1/Mpc**3\"))\n",
      "bbox = np.array([[0, 40], [0, 40], [0, 40]])\n",
      "ds = yt.load_uniform_grid(data, rhosmooth.shape, length_unit=\"Mpc\", bbox=bbox)\n",
      "\n",
      "# plot it\n",
      "# slc = yt.SlicePlot(ds, \"z\", [\"density\"])\n",
      "# slc.set_cmap(\"density\", \"Blues\")\n",
      "# slc.annotate_grids(cmap=None)\n",
      "# slc.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# extract isocontours from the density field\n",
      "# drange = ds.sphere(\"max\", (10.0, \"Mpc\")) # selects a 10 Mpc sphere around the max density\n",
      "drange = ds.all_data() # selects the entire data cube\n",
      "surface = ds.surface(drange, \"density\", 0.0)\n",
      "\n",
      "# export a .obj file \n",
      "surface.export_obj(\"contours\")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}