from __future__ import print_function
import astropy.io.fits as pyfits
import numpy as np
import os
import urllib
import astropy.visualization as viz
import matplotlib.pyplot as plt
%matplotlib inline
plt.rcParams['figure.figsize'] = (10.0, 10.0)

"""
Explanation: A First Look at an X-ray Image Dataset
Images are data. They can be 2D, from cameras, or 1D, from spectrographs, or 3D, from IFUs (integral field units). In each case, the data come packaged as an array of numbers, which we can visualize, and do calculations with.
Let's suppose we are interested in clusters of galaxies. We choose one, Abell 1835, and propose to observe it with the XMM-Newton space telescope. We are successful, we design the observations, and they are taken for us. Next: we download the data, and take a look at it.
Getting the Data
We will download our images from HEASARC, the online archive where XMM data are stored.
End of explanation
"""


targdir = 'a1835_xmm'
if not os.path.isdir(targdir):
    os.mkdir()

filenames = ('P0098010101M2U009IMAGE_3000.FTZ', 
             'P0098010101M2U009EXPMAP3000.FTZ',
             'P0098010101M2X000BKGMAP3000.FTZ')

remotedir = 'http://heasarc.gsfc.nasa.gov/FTP/xmm/data/rev0/0098010101/PPS/'

for filename in filenames:
    path = os.path.join(targdir, filename)
    url = os.path.join(remotedir, filename)
    if not os.path.isfile(path):
        urllib.urlretrieve(url, path)

imagefile, expmapfile, bkgmapfile = [os.path.join(targdir, filename) for filename in filenames]
        
for filename in os.listdir(targdir):
    print('{0:>10.2f} KB   {1}'.format(os.path.getsize(os.path.join(targdir, filename))/1024.0, filename))

"""
Explanation: Download the example data files if we don't already have them.
End of explanation
"""


imfits = pyfits.open(imagefile)
imfits.info()

"""
Explanation: The XMM MOS2 image
Let's find the "science" image taken with the MOS2 camera, and display it.
End of explanation
"""


im = imfits[0].data

"""
Explanation: imfits is a FITS object, containing multiple data structures. The image itself is an array of integer type, and size 648x648 pixels, stored in the primary "header data unit" or HDU. 

If we need it to be floating point for some reason, we need to cast it:
im = imfits[0].data.astype('np.float32')
Note that this (probably?) prevents us from using the pyfits "writeto" method to save any changes. Assuming the integer type is ok, just get a pointer to the image data.

Accessing the .data member of the FITS object returns the image data as a numpy ndarray.
End of explanation
"""


!ds9 -log "$imagefile"

"""
Explanation: Let's look at this with ds9.
End of explanation
"""


plt.imshow(viz.scale_image(im, scale='log', max_cut=40), cmap='gray', origin='lower');
plt.savefig("figures/cluster_image.png")

"""
Explanation: If you don't have the image viewing tool ds9, you should install it - it's very useful astronomical software. You can download it (later!) from this webpage.

We can also display the image in the notebook:
End of explanation
"""


im[350:359,350:359]

index = np.unravel_index(im.argmax(), im.shape)
print("image dimensions:",im.shape)
print("location of maximum pixel value:",index)
print("maximum pixel value: ",im[index])

"""
Explanation: Exercise
What is going on in this image? 
Make a list of everything that is interesting about this image with your neighbor, and we'll discuss the features you identify in about 5 minutes time.

Just to prove that images really are arrays of numbers:
End of explanation
"""


