added exception for bad headers

Metadata_extractor/0.1/Dockerfile

@@ -2,9 +2,8 @@ From ubuntu:14.04.3
 MAINTAINER HUANIAN ZHANG <fantasyzhn@email.arizona.edu>
 LABEL "This Dockerfile is for Metadata_extractor:ME_v01"
 
-RUN apt-get update && apt-get install -y git make gcc wget python-dev python2.7 python-numpy python-pip
-RUN pip install --no-deps astropy 
-RUN pip install fitsio==0.9.7
+RUN apt-get update && apt-get install -y git make gcc wget python-dev python2.7 python-numpy python-pip matplotlib
+RUN pip install --no-deps astropy
 ENV BINPATH /usr/bin
 ADD extract_metadata.py $BINPATH
 RUN chmod +x $BINPATH/extract_metadata.py

Metadata_extractor/0.1/extract_metadata.py

@@ -1,35 +1,133 @@
 #!/usr/bin/env python
-
-import fitsio
-import sys
-from astropy.wcs import WCS
-import csv
-from astropy.io import fits
 import os
+import sys
+
+import astropy.units as u
+import matplotlib.pylab as plt
 import numpy as np
 
+from astropy.coordinates import SkyCoord
+from astropy.io import fits
+from astropy.wcs import WCS
+from astropy.wcs.utils import proj_plane_pixel_scales
+from matplotlib.colors import LogNorm
+
+# FROM ASTROPY PYREGION, WILL ADD IMPORT WHEN ITS OUT OF DEV VERSION!!!
+def _calculate_rotation_angle(reg_coordinate_frame, header):
+    """Calculates the rotation angle from the region to the header's frame
+
+    This attempts to be compatible with the implementation used by SAOImage
+    DS9. In particular, this measures the rotation of the north axis as
+    measured at the center of the image, and therefore requires a
+    `~astropy.io.fits.Header` object with defined 'NAXIS1' and 'NAXIS2'
+    keywords.
+
+    Parameters
+    ----------
+    reg_coordinate_frame : str
+        Coordinate frame used by the region file
+
+    header : `~astropy.io.fits.Header` instance
+        Header describing the image
+
+    Returns
+    -------
+    y_axis_rot : float
+        Degrees by which the north axis in the region's frame is rotated when
+        transformed to pixel coordinates
+    """
+    new_wcs = WCS(header)
+    region_frame = SkyCoord(
+        '0d 0d',
+        frame=reg_coordinate_frame,
+        obstime='J2000')
+    region_frame = SkyCoord(
+        '0d 0d',
+        frame=reg_coordinate_frame,
+        obstime='J2000',
+        equinox=region_frame.equinox)
+
+    origin = SkyCoord.from_pixel(
+        header['NAXIS1']/2,
+        header['NAXIS2']/2,
+        wcs=new_wcs,
+        origin=1).transform_to(region_frame)
+
+    offset = proj_plane_pixel_scales(new_wcs)[1]
+
+    origin_x, origin_y = origin.to_pixel(new_wcs, origin=1)
+    origin_lon = origin.data.lon.degree
+    origin_lat = origin.data.lat.degree
+
+    offset_point = SkyCoord(
+        origin_lon, origin_lat+offset, unit='degree',
+        frame=origin.frame.name, obstime='J2000')
+    offset_x, offset_y = offset_point.to_pixel(new_wcs, origin=1)
+
+    north_rot = np.arctan2(
+        offset_y-origin_y,
+        offset_x-origin_x) / np.pi*180.
+
+    cdelt = new_wcs.wcs.get_cdelt()
+    if (cdelt > 0).all() or (cdelt < 0).all():
+        return north_rot - 90
+    else:
+        return -(north_rot - 90)
+
+public_folder = '/iplant/home/shared/Astrolabe/ViewinWWT/'
+base_url = 'http://www.worldwidetelescope.org/wwtweb/ShowImage.aspx?'
+
+def viewinwwt(filename, imageurl):
+    # dirpath = sys.argv[1]
+    # filenames = [f for f in os.listdir(filepath) if f.endswith('fits')]
+    fname = '{:s}.txt'.format(os.path.splitext(os.path.split(filename)[1])[0])
+    outfile = os.path.join(public_folder, fname)
+    header = fits.getheader(filename)
+    reqd = {'imageurl': imageurl}
+    try:
+        reqd['x'] = header['CRPIX1']
+        reqd['y'] = header['CRPIX2']
+        ra_str = header['RA']
+        dec_str = header['DEC']
+    except:
+        print('{:s} does not have the needed header keys'.format(filename))
+        return
+    coord = SkyCoord('{} {}'.format(ra_str, dec_str), unit=(u.hourangle, u.deg))
+    reqd['ra'] = coord.ra.value
+    reqd['dec'] = coord.dec.value
+
+    reqd['rotation'] = _calculate_rotation_angle('icrs', header)
+
+    wcs = WCS(header)
+    reqd['scale'] = proj_plane_pixel_scales(wcs)[0]
+    reqd['name'] = os.path.split(filename)[0]
+
+    request = 'name={name}&ra={ra}&dec={dec}&x={x}&y={y}&rotation={rotation}&imageurl={imageurl}'.format(**reqd)
+    with open(outfile, 'w') as outp:
+        outp.write('{0:s}{1:s}'.format(base_url, request))
+    return
+
+
+def makepng(filename):
+    header = fits.getheader(filename)
+    data = fits.getdata(filename)
+    try:
+        plt.imshow(data, cmap=plt.cm.Greys, norm=LogNorm(),
+                   vmin=header['IRAF-MIN'], vmax=header['IRAF-MAX'])
+    except KeyError:
+        plt.imshow(data, cmap=plt.cm.Greys, norm=LogNorm())
+
+    plt.axis('off')
+    plt.tight_layout()
+    fname = '{0:s}.png'.format(os.path.splitext(os.path.split(filename)[1])[0])
+    outfig = os.path.join(public_folder, fname)
+    plt.savefig(outfig)
+    return outfig
+
 
-#filepath = '/data2/fantasyzhn/Globular_cluster/Galex_image/'
-filepath = sys.argv[1]
-filename = []
-for file in os.listdir(filepath):
-    if file.endswith(".fits"):
-        filename.append(os.path.join(filepath, file))
-        #filename.append(file)
-
-keys = ['NAXIS','NAXIS1','NAXIS2','CDELT1','CDELT2','CTYPE1','CTYPE2','CRPIX1','CRPIX2','CRVAL1','CRVAL2']
-fieldsname = ['filepath','Number of axis','Dimension of axis1','Dimension of axis2','Interval of axis1','Interval of axis2','Projection type of axis1','Projection type of axis2','Central pixel of axis1','Central pixel of axis2','Value of central pixel1','Value of central pixel2']
-#writer = csv.DictWriter('metadata.csv',fieldnames=fieldsname,extrasaction='ignore',delimiter = ',')
-#writer.writeheader()
-writer=csv.writer(open('metadata.csv','wa'))
-writer.writerow(fieldsname)
-#for word in fieldsname:
-#    writer.writerow([word])
-for i, name in enumerate(filename):
-    fits_index = fitsio.FITS(name)
-    hdu = fits_index[0].read_header()
-    keys_value = [name]
-    for j in keys:
-        keys_value.append(hdu[j])
-    keys_value = np.array(keys_value)	
-    writer.writerow(keys_value)
+if __name__ == "__main__":
+    directory = sys.argv[1]
+    filenames = [f for f in os.path.listdir(directory) if f.endswith('fits')]
+    for filename in filenames:
+        imageurl = makepng(filename)
+        viewinwwt(filename, imageurl)