Determining the horizontal and vertical offsets between neighbouring ccd's using satellite tracks

The tracks of two or more satellites can be used to determine the horizontal and vertical offsets between ccd's if the satellites cross these ccd's in a single exposure.

In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from astropy.io import fits

plt.rcParams['image.cmap'] = 'viridis_r'

Let us first retrieve the RawScienceFrames for two neightbouring ccd's, in this case ESO_CCD_#75 and ESO_CCD_#76, for which we know there are two satellite tracks that span both ccd's (Thank's Gert!)

In [2]:
rawname = 'OMEGACAM.2011-08-31T09:35:54.666_{ext}.fits'
#RawScienceFrame(pathname=rawname.format(ext=7)).retrieve()
#RawScienceFrame(pathname=rawname.format(ext=8)).retrieve()
h7 = fits.open(rawname.format(ext=7))[0].header
h8 = fits.open(rawname.format(ext=8))[0].header
print h7['CRPIX1'] - h8['CRPIX1']
print h7['CRPIX2'] - h8['CRPIX2']
c7 = fits.open(rawname.format(ext=7))[0].data[:4100,48:2096]
c8 = fits.open(rawname.format(ext=8))[0].data[:4100,48:2096]

2142.2
-2.37

Now that we've retrieved the RawScienceFrames we can have a look and see where the satellite tracks are.

In [3]:
fig, (fig7,fig8) = plt.subplots(nrows=1,ncols=2,figsize=(14,8))

fig7.set_title(h7['CHIP_ID'])
fig7.imshow(c7,vmin=250,vmax=420)
fig7.annotate('t1',xy=(1070,3410),xycoords='data',xytext=(-30,-50),textcoords='offset points',size=20,arrowprops=dict(arrowstyle='->'))
fig7.annotate('t1a',xy=(0,4000),xycoords='data',xytext=(30,50),textcoords='offset points',size=10,arrowprops=dict(arrowstyle='->'))
fig7.annotate('t1b',xy=(2048,2850),xycoords='data',xytext=(-30,-50),textcoords='offset points',size=10,arrowprops=dict(arrowstyle='->'))
fig7.annotate('t2',xy=(2048,100),xycoords='data',xytext=(-50,15),textcoords='offset points',size=20,arrowprops=dict(arrowstyle='->'))
fig8.set_title(h8['CHIP_ID'])
fig8.imshow(c8,vmin=250,vmax=420)
fig8.annotate('t1',xy=(1070,2160),xycoords='data',xytext=(-30,-40),textcoords='offset points',size=20,arrowprops=dict(arrowstyle='->'))
fig8.annotate('t1c',xy=(0,2790),xycoords='data',xytext=(30,50),textcoords='offset points',size=10,arrowprops=dict(arrowstyle='->'))
fig8.annotate('t1d',xy=(2048,1600),xycoords='data',xytext=(-30,-50),textcoords='offset points',size=10,arrowprops=dict(arrowstyle='->'))
fig8.annotate('t2',xy=(400,1300),xycoords='data',xytext=(50,-50),textcoords='offset points',size=20,arrowprops=dict(arrowstyle='->'))
Out[3]:
<matplotlib.text.Annotation at 0x7f1484d5ee90>

As you can see, there are two satellite tracks:

  • t1, starting a the top left of ESO_CCD_#75 and ending at the lower right half of ESO_CCD_#76.
  • t2, starting a the bottom right of ESO_CCD_#75 and ending at the top of ESO_CCD_#76.

Assuming that the satellite tracks continue from one ccd to the next and they form straight lines, there are a few things we need to know to determine the relative position of the ccd's.

Label the points where each individual satellite track intersects the borders of the ccd's, from left to right: a, b for ESO_CCD_#75 and c, d for ESO_CCD_#76. Then we have

  • The locations (t1ax,t1ay), (t1bx,t1by) that identify the left segment of t1.
  • The locations (t1cx,t1cy), (t1dx,t1dy) that identify the right segment of t1.
  • The locations (t2ax,t2ay), (t2bx,t2by) that identify the left segment of t2.
  • The locations (t2cx,t2cy), (t2dx,t2dy) that identify the right segment of t2.
In [4]:
def find_horizontal_intersection(ccd, xslice, ypos, subplot=None):
    # Note that in the numpy array the y-axis comes first and the x-axis follows.
    crosscut = ccd[ypos,xslice]
    # Find the location of the peak. Assume it's the local maximum in the cross-cut
    peakposx = xslice.start + np.argmax(crosscut)
    if subplot > 0:
        ax = plt.subplot(subplot)
        # Show the cross-cut and mark its local maximum
        ax.plot(np.arange(xslice.start,xslice.stop), crosscut, '.')
        ax.plot([peakposx,peakposx], ax.get_ylim(), 'r')
    return float(peakposx), float(ypos)

def find_vertical_intersection(ccd, xpos, yslice, subplot=None):
    # Note that in the numpy array the y-axis comes first and the x-axis follows.
    crosscut = ccd[yslice,xpos]
    # Find the location of the peak. Assume it's the local maximum in the cross-cut
    peakposy = yslice.start + np.argmax(crosscut)
    if subplot > 0:
        # Show the cross-cut and mark its local maximum
        ax = plt.subplot(subplot)
        ax.plot(np.arange(yslice.start,yslice.stop), crosscut, '.')
        ax.plot([peakposy,peakposy], ax.get_ylim(), 'r')
    return float(xpos), float(peakposy)

f, figs = plt.subplots(nrows=4,ncols=2,figsize=(14,8))

t1ax,t1ay = find_vertical_intersection(c7, 0, slice(3910,4099), subplot=421)
t1bx,t1by = find_vertical_intersection(c7, 2047, slice(2700,2900), subplot=422)
t1cx,t1cy = find_vertical_intersection(c8, 0, slice(2650,2850), subplot=423)
t1dx,t1dy = find_vertical_intersection(c8, 2045, slice(1450,1650), subplot=424)
t2ax,t2ay = find_horizontal_intersection(c7, slice(1900,2048), 0, subplot=425)
t2bx,t2by = find_vertical_intersection(c7, 2047, slice(0,200), subplot=426)
t2cx,t2cy = find_vertical_intersection(c8, 0, slice(300,500), subplot=427)
t2dx,t2dy = find_horizontal_intersection(c8, slice(1300,1500), 4099, subplot=428)
t2ax += 2
t2by -= 3
t2cy += 1
#t1ay+=130
#t1by+=130
#t2ay+=130
#t2by+=130
In [5]:
plt.figure(figsize=(12,8))

dx = 2048 - (-88.5)
dy = - (-3.59)
assert t1bx == t2bx
assert t1cx == t2cx

a1c, b1c = np.polyfit([t1cx,t1dx], [t1cy,t1dy], 1)
a2c, b2c = np.polyfit([t2cx,t2dx], [t2cy,t2dy], 1)

a1a, b1a = np.polyfit([t1ax,t1bx], [t1ay,t1by], 1)
a2a, b2a = np.polyfit([t2ax,t2bx], [t2ay,t2by], 1)

plt.plot([t1ax,t1bx],[t1ay,t1by],'b')
plt.plot([t2ax,t2bx],[t2ay,t2by],'g')
plt.plot([dx+t1cx,dx+t1dx],[dy+t1cy,dy+t1dy],'b')
plt.plot([dx+t2cx,dx+t2dx],[dy+t2cy,dy+t2dy],'g')
plt.plot([dx+t1cx-88.5,dx+t1cx],[dy+a1c*(t1cx-88.5)+b1c,dy+t1cy],'o-r')
plt.plot([dx+t2cx-88.5,dx+t2cx],[dy+a2c*(t2cx-88.5)+b2c,dy+t2cy],'o-r')

print 'The slopes in t1 should be similar in each ccd: {a1a:5.2f}, {a1c:5.2f}'.format(**vars())
print 'The slopes in t2 should be similar in each ccd: {a2a:5.2f}, {a2c:5.2f}'.format(**vars())
x = ((t1by - t2by) - (b1c - b2c)) / (a1c - a2c)
y = a1c*x + b1c - t1by
print a1c*x+b1c-t1by, a2c*x+b2c-t2by
print 'Horizontal and vertical offset of ESO_CCD_#75 w.r.t. ESO_CCD_#76: {xoffset:6.2f}, {yoffset:5.2f}'.format(xoffset=x,yoffset=y)
# t1by, a1c*x + b1c + 3.4, t2by, a2c*x + b2c + 3.4
x = ((t1cy - t2cy) - (b1a - b2a)) / (a1a - a2a)
y = a1a*x + b1a - t1cy
print 'Horizontal and vertical offset of ESO_CCD_#76 w.r.t. ESO_CCD_#75: {offset:6.2f}, {yoffset:5.2f}'.format(offset=x-2048,yoffset=y)
import mpld3
mpld3.display()

The slopes in t1 should be similar in each ccd: -0.58, -0.58
The slopes in t2 should be similar in each ccd:  2.75,  2.74
-3.58928460625 -3.58928460625
Horizontal and vertical offset of ESO_CCD_#75 w.r.t. ESO_CCD_#76: -88.50, -3.59
Horizontal and vertical offset of ESO_CCD_#76 w.r.t. ESO_CCD_#75:  87.23,  3.62
Out[5]:
In [ ]: