Python库 | testreverseautoslitcode-0.1.13.tar.gz

#Create the quick look plot of the mask results and save as a pdf import subprocess import os from astropy.io import fits from astropy.wcs import WCS,utils from astropy.table import Table from matplotlib import pyplot as plt from astropy.coordinates import SkyCoord, ICRS, Galactic, FK4, FK5 import numpy as np import astropy.units as u from matplotlib.patches import Polygon as Pgon from astropy.visualization import PercentileInterval, ImageNormalize from glob import glob from matplotlib.gridspec import GridSpec #Function to create the quick look plots of teh mask results def create_err_plot(data,catalog_obj_ra,catalog_obj_dec,objects_ra,objects_dec,output_file): #Calculate the center of the mask to provide PanSTAMPS to download the image for the mask field (should I just send over the actual center coords) x_center = np.average(data['Calc_RA']) y_center = np.average(data['Calc_Dec']) #Download the PanSTARRS field image for the mask location panstampscommand = 'panstamps --width=15 stack '+str(x_center)+' '+str(y_center) process2 = subprocess.Popen(panstampscommand.split(), stdout=subprocess.PIPE) output, error = process2.communicate() process2.wait() #Open the field image. Just picks the first image on the list generated by glob field_file_path = str(x_center)+'*/''*.fits' all_field_files = glob(field_file_path) hdu = fits.open(all_field_files[0])[0] #================================================================================================================================= #Set up the grid for the plot panels using GridSpec wcs = WCS(hdu.header) fig = plt.figure(constrained_layout=True,figsize=(7,4)) gs = GridSpec(3, 3, figure=fig,width_ratios = [3,1,1],height_ratios=[1,1,1]) ax1 = fig.add_subplot(gs[:, :-2],projection=wcs) ax2 = fig.add_subplot(gs[0, 1],projection=wcs) ax3 = fig.add_subplot(gs[1, 1],projection=wcs) ax4 = fig.add_subplot(gs[2, 1],projection=wcs) ax5 = fig.add_subplot(gs[0, 2],projection=wcs) ax6 = fig.add_subplot(gs[1, 2],projection=wcs) ax7 = fig.add_subplot(gs[2, 2],projection=wcs) #For clarity of the plot, set the ticks of the zoom in panels on the alignment boxes to be invisible lon2 = ax2.coords[0] lat2 = ax2.coords[1] lon2.set_ticks_visible(False) lat2.set_ticklabel_visible(False) lat2.set_ticks_visible(False) lon2.set_ticklabel_visible(False) lon3 = ax3.coords[0] lat3 = ax3.coords[1] lon3.set_ticks_visible(False) lat3.set_ticklabel_visible(False) lat3.set_ticks_visible(False) lon3.set_ticklabel_visible(False) lon4 = ax4.coords[0] lat4 = ax4.coords[1] lon4.set_ticks_visible(False) lat4.set_ticklabel_visible(False) lat4.set_ticks_visible(False) lon4.set_ticklabel_visible(False) lon5 = ax5.coords[0] lat5 = ax5.coords[1] lon5.set_ticks_visible(False) lat5.set_ticklabel_visible(False) lat5.set_ticks_visible(False) lon5.set_ticklabel_visible(False) lon6 = ax6.coords[0] lat6 = ax6.coords[1] lon6.set_ticks_visible(False) lat6.set_ticklabel_visible(False) lat6.set_ticks_visible(False) lon6.set_ticklabel_visible(False) lon7 = ax7.coords[0] lat7 = ax7.coords[1] lon7.set_ticks_visible(False) lat7.set_ticklabel_visible(False) lat7.set_ticks_visible(False) lon7.set_ticklabel_visible(False) #================================================================================================================================= #Create the main plot. This is the full mask view. norm = ImageNormalize(hdu.data, interval=PercentileInterval(99)) norm2 = ImageNormalize(hdu.data, interval=PercentileInterval(99)) i=0 avg_y_pos = [] box_numbers = [] #order the slits by y positions while i in range(len(data)): avg_y_pos.append(np.average(data['Y'][i:i+4])) i = i+4 avg_y_pos = np.sort(avg_y_pos) i=0 while i in range(int(len(data))): vertices = np.zeros(shape=(4,2)) averages = np.zeros(shape=(2,2)) vertices[0:4,0] = data['Calc_RA'][i:i+4] vertices[0:4,1] = data['Calc_Dec'][i:i+4] averages[0,0]=np.average((vertices[1,0],vertices[2,0])) averages[1,0] = np.average((vertices[0,0],vertices[3,0])) averages[0,1]=np.average((vertices[1,1],vertices[2,1])) averages[1,1]=np.average((vertices[0,1],vertices[3,1])) ra_center = np.average(vertices[:,0]) dec_center = np.average(vertices[:,1]) y_avg = np.average(data['Y'][i:i+4]) for k in range(len(avg_y_pos)): if y_avg == avg_y_pos[k]: ind = k+1 ra1,ra2,ra3,ra4 = data['Calc_RA'][i],data['Calc_RA'][i+1],data['Calc_RA'][i+2],data['Calc_RA'][i+3] dec1,dec2,dec3,dec4 = data['Calc_Dec'][i],data['Calc_Dec'][i+1],data['Calc_Dec'][i+2],data['Calc_Dec'][i+3] side1 = 3600*np.sqrt(((ra2-ra1)*np.cos(dec1*np.pi/180))**2+(dec1-dec2)**2) side2 = 3600*np.sqrt(((ra2-ra3)*np.cos(dec2*np.pi/180))**2+(dec3-dec2)**2) ratio = side1/side2 if 0.9 < ratio < 1.1: box_numbers.append(ind) k=k+1 pixelcoords = [] for l in range(0,4): temp = SkyCoord(vertices[l,0],vertices[l,1], unit=(u.deg,u.deg),frame='fk5') temp2 = utils.skycoord_to_pixel(temp, wcs=wcs, origin=1, mode='all') pixelcoords.append(temp2) r = Pgon(pixelcoords, closed = True, edgecolor='green', facecolor='none',alpha=0.6,zorder=4) ax1.add_patch(r) ax1.text(ra_center,dec_center,str(ind),color='blue',transform=ax1.get_transform('fk5')) i = i+4 im = ax1.imshow(hdu.data,origin='lower',zorder=0,cmap='gray_r',norm=norm) temp1a = SkyCoord(np.min(data['Calc_RA'])-50/3600.,np.min(data['Calc_Dec'])-50/3600.,unit=(u.deg,u.deg)) temp1b = SkyCoord(np.max(data['Calc_RA'])+50/3600.,np.max(data['Calc_Dec'])+50/3600.,unit=(u.deg,u.deg)) pixels1a=utils.skycoord_to_pixel(temp1a, wcs=wcs, origin=0, mode='all') pixels1b=utils.skycoord_to_pixel(temp1b, wcs=wcs, origin=0, mode='all') ax1.set_xlim([pixels1a[0],pixels1b[0]]) ax1.set_ylim([pixels1a[1],pixels1b[1]]) ax1.invert_xaxis() #================================================================================================================================= #Find the first six boxes and plot the cutout panels in the quick look plot. i=0 boxes=0 while i <= len(data)-4 and boxes < 6: ra1,ra2,ra3,ra4 = data['Calc_RA'][i],data['Calc_RA'][i+1],data['Calc_RA'][i+2],data['Calc_RA'][i+3] dec1,dec2,dec3,dec4 = data['Calc_Dec'][i],data['Calc_Dec'][i+1],data['Calc_Dec'][i+2],data['Calc_Dec'][i+3] side1 = 3600*np.sqrt(((ra2-ra1)*np.cos(dec1*np.pi/180))**2+(dec1-dec2)**2) side2 = 3600*np.sqrt(((ra2-ra3)*np.cos(dec2*np.pi/180))**2+(dec3-dec2)**2) ratio = side1/side2 if 0.9 < ratio < 1.1: boxes = boxes+1 vertices = np.zeros(shape=(4,2)) vertices[0:4,0] = data['Calc_RA'][i:i+4] vertices[0:4,1] = data['Calc_Dec'][i:i+4] ra_buffer = (5/(3600.*np.cos(vertices[0,1]*np.pi/180.))) dec_buffer = 5/3600. if boxes == 1: pixelcoords = [] for l in range(len(vertices)): temp = SkyCoord(vertices[l,0],vertices[l,1], unit=(u.deg,u.deg),frame='fk5') temp2 = utils.skycoord_to_pixel(temp, wcs=wcs, origin=1, mode='all') pixelcoords.append(temp2) r = Pgon(pixelcoords, closed = True, edgecolor='green', facecolor='none',alpha=0.6,zorder=4) ax2.add_patch(r) im2 = ax2.imshow(hdu.data,origin='lower',zorder=0,cmap='gray_r',norm=norm2) temp1a = SkyCoord(np.min(vertices[0:4,0])-ra_buffer,np.min(vertices[0:4,1])-dec_buffer,unit=(u.deg,u.