coupes_fits.py
8.59 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Sun Mar 15 16:37:04 2020
@author: sacha
"""
def cut_fits(filename,coords,plot_all_cut,average,model_output=[],fig_index=1):
'''Function allowing to make a certain number of cuts between 2 points specified
- can computes error bars by averaging the surronding pixels
- plots the map and the corresponding cuts
input:
filename : str, name of the .fits file in which the cuts are made
coords : list with the element arranged in 3 parts, where 'n' is the number of cuts:
- first list: RA/DEC coordinates of the n starting points, coordinates in deg., each RA/DEC couple is a list,
- second list: RA/DEC coordinates of the n ending points, coordinates in deg.,
- third item: RA/DEC coordinates of the reference point, ex: a star in the case where the cut starts somewhere else
- example with n=3 and the reference being HD200775: coords=[[coord1_1,coord1_2,coord1_3],[coord2_1,coord2_2,coord2_3],coordstar]
plot_all_cut: bool, True if you want to plot the figures,
average: bool, True if you want the average of a 3x3 square and error bars, False if you just want the pixel values
model_output: 2xN array, specific to the photoelectric heating effeciency obtained by the model, allowing to plot these values with the cuts.
leave it as an empty list if you cut something else than the PE efficiency map.
fig_index: index of the figure if you do several plots, default to 1
returns list of the data along the cuts, the corresponding errors and the distance in arcsec.
'''
import numpy as np
import matplotlib.pyplot as plt
from astropy.io import fits
from astropy.wcs import WCS
# from matplotlib.patches import Arrow, Circle
# ==============================
#-- get the data, header and pixel scale in arcsec
hdu=fits.open(filename)[0]
data=hdu.data
hdr=hdu.header
w_ = WCS(hdr)
taille_deg0=hdr['CDELT1'] #taille en degres
taille_sec0=taille_deg0*3600. #bonne taille en arc sec
#-- pixel of reference for the computation of the difference in arcsec
coordstar=[coords[2][0],coords[2][1]]
pixs = w_.wcs_world2pix([coordstar],0)
xs=pixs[0,0].copy()
ys=pixs[0,1].copy()
all_data_average=[]
all_distance_average=[]
all_error=[]
cut_list=[]
#-- copy the data for different purposes
z = data.copy()
# zz=z.copy()
if plot_all_cut:
'''plot of the map to cut'''
fig = plt.figure(fig_index,figsize=(7.48,5.9))
# ax = fig.add_subplot(211)
# pos=ax.imshow(zz[:,::-1],cmap='gray',origin='lower')
# ax.axis('image')
# fig.colorbar(pos,ax=ax,extend='both')
ax2 = fig.add_subplot(111)
ax2.set_xlabel('distance from the star (")')
ax2.set_ylabel('$\epsilon_{PAH}$')
#-- Extract the line...
# Make a line with "num" points...
#-- Loop on the number of cuts wanted
for ind in range(len(coords[0])):
#-- Extract the coordinate in deg from coords list and convert it into pixel index...
#... of both extreme points of the cut. From pix0 to pix1
coord0=[coords[0][ind][0],coords[0][ind][1]]
coord1=[coords[1][ind][0],coords[1][ind][1]]
pix0 = w_.wcs_world2pix([coord0],0) # These are in _pixel_ coordinates!!
pix1 = w_.wcs_world2pix([coord1],0)
x0=pix0[0,0].copy()
y0=pix0[0,1].copy()
x1=pix1[0,0].copy()
y1=pix1[0,1].copy()
#-- Length of the cut in pixels
length = int(round(np.sqrt((x0-x1)**2+(y0-y1)**2)))
# distance=taille_sec0*np.sqrt((pix_list[:,0]-xs)**2+(pix_list[:,1]-ys)**2)
#-- Creation of vectors x and y of a length of "length" between the extreme positions...
#... and make sure that they are integers
x, y = np.linspace(x0, x1, length), np.linspace(y0, y1, length)
x_,y_=x.round(),y.round()
x_int,y_int=x_.astype(np.int), y_.astype(np.int)
#-- coordinates of each pixels
coord=np.zeros([len(x_int),2])
coord[:,0]=x_int
coord[:,1]=y_int
# Extract the values along the line
distance_average=[]
data_average=[]
pixel_value=[]
error=[]
list_of_z=[]
list_of_z2=[]
plouf=[]
g=0
#-- Loop on the pixels inside the cut
for xi,yi in coord:
#-- fix the pixel position and extract the value
xi,yi=xi.astype(np.int), yi.astype(np.int)
#-- distance from the reference point
distance=taille_sec0*np.sqrt((xi-xs)**2+(yi-ys)**2)
distance_average.append(distance)
#-- if you want a value with assymetric errorbars around a mean value
#... average the value of the pixel with all the surronding pixels
if average:
if (yi==z.shape[0]-1):
list_of_z.append([z[yi,xi],z[yi,xi+1],z[yi-1,xi+1],z[yi-1,xi],z[yi-1,xi-1],z[yi,xi-1]])
if np.all(np.isnan(list_of_z)==True):
data_average.append(np.mean(np.array([z[yi,xi],z[yi,xi+1],z[yi-1,xi+1],z[yi-1,xi],z[yi-1,xi-1],z[yi,xi-1]])))
error.append(np.std(np.array([z[yi,xi],z[yi,xi+1],z[yi-1,xi+1],z[yi-1,xi],z[yi-1,xi-1],z[yi,xi-1]])))
else:
for plouf in list_of_z[0]:
g+=1
if np.isnan(plouf)==False:
list_of_z2.append(plouf)
# list_of_z2gamma.append(list_of_zgamma[0][g])
data_average.append(np.mean(list_of_z2))
# gammaverage.append(np.mean(list_of_z2gamma))
error.append(np.std(list_of_z2))
else:
list_of_z.append([z[yi,xi],z[yi+1,xi-1],z[yi+1,xi],z[yi+1,xi+1],z[yi,xi+1],z[yi-1,xi+1],z[yi-1,xi],z[yi-1,xi-1],z[yi,xi-1]])
# list_of_zgamma.append([data_gamma[yi,xi],data_gamma[yi,xi+1],data_gamma[yi-1,xi+1],data_gamma[yi-1,xi],data_gamma[yi-1,xi-1],data_gamma[yi,xi-1]])
if np.all(np.isnan(list_of_z)==True):
data_average.append(np.mean(np.array([z[yi,xi],z[yi+1,xi-1],z[yi+1,xi],z[yi+1,xi+1],z[yi,xi+1],z[yi-1,xi+1],z[yi-1,xi],z[yi-1,xi-1],z[yi,xi-1]])))
error.append(np.std(np.array([z[yi,xi],z[yi+1,xi-1],z[yi+1,xi],z[yi+1,xi+1],z[yi,xi+1],z[yi-1,xi+1],z[yi-1,xi],z[yi-1,xi-1],z[yi,xi-1]])))
# gammaverage.append(np.mean(np.array([data_gamma[yi,xi],data_gamma[yi,xi+1],data_gamma[yi-1,xi+1],data_gamma[yi-1,xi],data_gamma[yi-1,xi-1],data_gamma[yi,xi-1]])))
else:
for plouf in list_of_z[0]:
if np.isnan(plouf)==False:
list_of_z2.append(plouf)
data_average.append(np.mean(list_of_z2))
error.append(np.std(list_of_z2))
list_of_z=[]
list_of_z2=[]
plouf=[]
else:
pixel_value.append(z[yi,xi])
#-- save the vector of distances of the current cut
all_distance_average.append(distance_average)
if average:
#-- If you want just the value of the pixel with no error bar
all_data_average.append(np.array(data_average) )
all_error.append(error)
if plot_all_cut:
ax2.errorbar(np.array(all_distance_average[ind]),all_data_average[ind],yerr=all_error[ind],label='{}'.format(ind))#,label='cut{}'.format(ind))
ax2.legend()
else:
cut_list.append(pixel_value)
if plot_all_cut:
ax2.semilogy(all_distance_average[ind],cut_list[ind])
#-- Plot
if model_output!=[]:
eff=model_output[1,:].copy()
dist=model_output[0,:].copy()
ax2.semilogy(dist,eff,'*',color='red',label='modèle, G=0.5')
ax2.axvline(42.8,lw=0.8,ls='--',color='blue')
ax2.axvline(49.7,lw=0.8,ls='--',color='orange')
ax2.axvline(37.2,lw=0.8,ls='--',color='green')
ax2.axhline(0.0164,lw=0.8,ls='--',color='green')
plt.legend()
return all_data_average,all_error,all_distance_average