extracted SEDs, images, animations

Author: koepferl

getanimation.py

@@ -0,0 +1,49 @@
+import os
+
+import numpy as np
+
+import matplotlib as mpl
+mpl.use('Agg')
+import matplotlib.pyplot as plt
+
+from hyperion.model import ModelOutput
+from hyperion.util.constants import pc
+
+# Create output directory if it does not already exist
+if not os.path.exists('frames'):
+    os.mkdir('frames')
+
+# Open model
+m = ModelOutput('tutorial_model.rtout')
+
+# Read image from model
+wav, nufnu = m.get_image(group=2, distance=300 * pc)
+
+# nufnu is now an array with four dimensions (n_view, n_wav, n_y, n_x)
+
+# Fix the wavelength to the first one and cycle through viewing angles
+iwav = 0
+print "Wavelength is %g microns" % wav[iwav]
+
+for iview in range(nufnu.shape[0]):
+
+    # Open figure and create axes
+    fig = plt.figure()
+    ax = fig.add_subplot(1, 1, 1)
+
+    # This is the command to show the image. The parameters vmin and vmax are
+    # the min and max levels for the grayscale (remove for default values).
+    # The colormap is set here to be a heat map. Other possible heat maps
+    # include plt.cm.gray (grayscale), plt.cm.gist_yarg (inverted grayscale),
+    # plt.cm.jet (default, colorful). The np.sqrt() is used to plot the
+    # images on a sqrt stretch.
+    ax.imshow(np.sqrt(nufnu[iview, :, :, iwav]), vmin=0, vmax=np.sqrt(1.e-11), \
+              cmap=plt.cm.gist_heat, origin='lower')
+
+    # Save figure. The facecolor='black' and edgecolor='black' are for
+    # esthetics, and hide the axes
+    fig.savefig('frames/frame_%05i.png' % iview, \
+                facecolor='black', edgecolor='black')
+
+    # Close figure
+    plt.close(fig)

getfits.py

@@ -0,0 +1,20 @@
+import pyfits
+
+from hyperion.model import ModelOutput
+from hyperion.util.constants import pc
+
+# Open the model - we specify the name without the .rtout extension
+m = ModelOutput('tutorial_model.rtout')
+
+# Extract the image for the first inclination, and scale to 300pc. We
+# have to specify group=1 as there is no image in group 0
+wav, nufnu = m.get_image(group=1, inclination=0, distance=300 * pc)
+
+# The image extracted above is a 3D array. We can write it out to FITS.
+# We need to swap some of the directions around so as to be able to use
+# the ds9 slider to change the wavelength of the image.
+pyfits.writeto('image_cube.fits', nufnu.swapaxes(0, 2).swapaxes(1, 2), \
+               clobber=True)
+
+# We can also just output one of the wavelengths
+pyfits.writeto('image_slice.fits', nufnu[:, :, 0], clobber=True)

getfitsWCS.py

@@ -0,0 +1,31 @@
+import pywcs
+import pyfits
+
+from hyperion.model import ModelOutput
+from hyperion.util.constants import pc
+
+m = ModelOutput('tutorial_model.rtout')
+wav, nufnu = m.get_image(group=1, inclination=0, distance=300 * pc)
+
+# Initialize WCS information
+wcs = pywcs.WCS(naxis=2)
+
+# Use the center of the image as projection center
+wcs.wcs.crpix = [nufnu.shape[2] / 2. + 0.5,
+                 nufnu.shape[1] / 2. + 0.5]
+
+# Set the coordinates of the image center
+wcs.wcs.crval = [233.4452, 1.2233]
+
+# Set the pixel scale (in deg/pix)
+wcs.wcs.cdelt = [1./3600., 1./3600.]
+
+# Set the coordinate system
+wcs.wcs.ctype = ['GLON-CAR', 'GLAT-CAR']
+
+# And produce a FITS header
+header = wcs.to_header()
+
+# Write out to a file including the new header
+pyfits.writeto('image_slice_wcs.fits', nufnu[:, :, 0], header,
+               clobber=True)

getsed.py

@@ -0,0 +1,16 @@
+import numpy as np
+
+from hyperion.model import ModelOutput
+from hyperion.util.constants import pc
+
+# Open the model - we specify the name without the .rtout extension
+m = ModelOutput('tutorial_model.rtout')
+
+# Extract the SED for the smallest inclination and largest aperture, and
+# scale to 300pc. In Python, negative indices can be used for lists and
+# arrays, and indicate the position from the end. So to get the SED in the
+# largest aperture, we set aperture=-1.
+wav, nufnu = m.get_sed(inclination=1, aperture=-1, distance=300 * pc)
+
+# Write out the SED to file
+np.savetxt('sed.txt', zip(wav, nufnu), fmt="%11.4e %11.4e")