From 29c549a1b1ec7de2694d8f8781479b2ca331c967 Mon Sep 17 00:00:00 2001 From: Phil Rosenfield Date: Mon, 25 Sep 2017 17:25:42 -0700 Subject: [PATCH 1/2] fixed up code for hunter data --- Metadata_extractor/0.1/Dockerfile | 5 +- Metadata_extractor/0.1/extract_metadata.py | 155 +++++++++++++++++---- 2 files changed, 127 insertions(+), 33 deletions(-) diff --git a/Metadata_extractor/0.1/Dockerfile b/Metadata_extractor/0.1/Dockerfile index 1abc5b3..927ff0e 100644 --- a/Metadata_extractor/0.1/Dockerfile +++ b/Metadata_extractor/0.1/Dockerfile @@ -2,9 +2,8 @@ From ubuntu:14.04.3 MAINTAINER HUANIAN ZHANG 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 diff --git a/Metadata_extractor/0.1/extract_metadata.py b/Metadata_extractor/0.1/extract_metadata.py index 4e8a3bc..22cc99b 100644 --- a/Metadata_extractor/0.1/extract_metadata.py +++ b/Metadata_extractor/0.1/extract_metadata.py @@ -1,35 +1,130 @@ #!/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} + reqd['x'] = header['CRPIX1'] + reqd['y'] = header['CRPIX2'] + ra_str = header['RA'] + dec_str = header['DEC'] + + 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) From 05433bf7397050ec187a0116a93e71db74406175 Mon Sep 17 00:00:00 2001 From: Phil Rosenfield Date: Fri, 29 Sep 2017 08:52:08 -0400 Subject: [PATCH 2/2] added exception for bad headers --- Metadata_extractor/0.1/extract_metadata.py | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/Metadata_extractor/0.1/extract_metadata.py b/Metadata_extractor/0.1/extract_metadata.py index 22cc99b..3a31375 100644 --- a/Metadata_extractor/0.1/extract_metadata.py +++ b/Metadata_extractor/0.1/extract_metadata.py @@ -84,11 +84,14 @@ def viewinwwt(filename, imageurl): outfile = os.path.join(public_folder, fname) header = fits.getheader(filename) reqd = {'imageurl': imageurl} - reqd['x'] = header['CRPIX1'] - reqd['y'] = header['CRPIX2'] - ra_str = header['RA'] - dec_str = header['DEC'] - + 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