In this exercise, we will:
Keys: NumPy, array indexing, slicing, view and reduction, astropy.stats and masked arrays
#Open the file
from astropy.io import fits
with fits.open('data/mods1r.20171125.0026.fits') as hdulist:
#Retrieve image data copy
data = hdulist[0].data.copy()
data
array([[ 0, 1373, 488, ..., 981, 2860, 0], [ 645, 317, 322, ..., 320, 302, 636], [ 644, 321, 324, ..., 320, 305, 639], ..., [ 645, 304, 322, ..., 335, 324, 670], [ 639, 302, 323, ..., 334, 323, 665], [ 0, 1332, 499, ..., 756, 1942, 0]], dtype=uint16)
np.array
is a grid of valuesshape
of an array is a tuple of integers giving the size of the array along each dimensionimport numpy as np
# 1 dimensional array
a = np.array([10., 20., 30., 40., 50.])
# 2 dimensional array
b = np.array([[ 10., 20., 30., 40., 50.], \
[ 60., 70., 80., 90., 100.]])
# 3 dimensional array
c = np.array([[[ 10., 20., 30., 40., 50.], \
[ 60., 70., 80., 90., 100.]],\
[[ 110., 120., 130., 140., 150.], \
[ 160., 170., 180., 190., 200.]],\
[[ 210., 220., 230., 240., 250.], \
[ 260., 270., 280., 290., 300.]]])
# and so on......
#1D array
a_shape = np.arange(0, 15)
#Create 2D array by reshape
b_shape = a_shape.reshape((3, 5))
BE CAREFUL: (we will see later) this is not a new array, but a view
of the the original array
b_shape
array([[ 0, 1, 2, 3, 4], [ 5, 6, 7, 8, 9], [10, 11, 12, 13, 14]])
We can also reshape in place
a_shape.shape=(5,3)
a_shape
array([[ 0, 1, 2], [ 3, 4, 5], [ 6, 7, 8], [ 9, 10, 11], [12, 13, 14]])
There are several general mechanisms for creating arrays
import numpy as np
#Array full of 1.0 of length 2
np.ones((2,))
#Empty array of length 4
np.empty((4,))
#Array full of -infinite of length
np.full((5,), -np.inf)
z=np.zeros((3,))
#Array with the same shape of z full of 3.14
np.full_like(z, 3.14)
#Values are generated between the half-open interval [-2.3,1.5) 0.4 is the step
np.arange(-2.3, 1.5, 0.4)
#Returns 5 evenly spaced numbers over the closed interval [-2.3, 1.5]
np.linspace(-2.3, 1.5, 5)
np.linspace(-2.3, 1.5, 5)
array([-2.3 , -1.35, -0.4 , 0.55, 1.5 ])
import numpy as np
a = np.array([-1, 0, 3])
.dtype
attribute# Data type
data.dtype.type ### All items have the same type: np.array are homogenus!
# Data type name
data.dtype.name
#Memory size (byte) of each elements
data.dtype.itemsize
a
array([-1, 0, 3])
#a = np.array([-1, 0, 3])
b = np.array([-1.,0.,3.])
a.dtype, b.dtype
(dtype('int64'), dtype('float64'))
c=np.array(a, dtype=bool)
c
array([ True, False, True])
a = 0.8*np.arange(0, 5) #Define the array
a[0] # First element
a[4] # Last element
3.2
a[4]
3.2
a[-1] # Last element
a[-5] # First element
0.0
a[-5]
0.0
#Get 2nd and 3th element [1,3)
a[1:3]
array([0.8, 1.6])
a[1:3]
array([0.8, 1.6])
Basic slicing
extends Python’s basic concept of slicing to N dimensions. Basic slicing occurs when obj is a slice
object (constructed by start:stop:step
notation inside of brackets)
b=a
b is a # True
True
c = a[3:5]
c+=1
a #changed
array([0. , 0.8, 1.6, 3.4, 4.2])
Just define the array for examples
a=10.*np.arange(0, 40).reshape((5,8))
a
array([[ 0., 10., 20., 30., 40., 50., 60., 70.], [ 80., 90., 100., 110., 120., 130., 140., 150.], [160., 170., 180., 190., 200., 210., 220., 230.], [240., 250., 260., 270., 280., 290., 300., 310.], [320., 330., 340., 350., 360., 370., 380., 390.]])
row,col
indexing a[0][1] #Element of 1st row, 2nd column
a[0,1] #The same
10.0
NOTE: arrays are used to map image data, but arrays indexing is different from image (x,y) coordinates
a[0][1]
10.0
a[0] #The 1st row
a[0,:]
array([ 0., 10., 20., 30., 40., 50., 60., 70.])
a[0]
array([ 0., 10., 20., 30., 40., 50., 60., 70.])
a[:,1] #The 2nd column
array([ 10., 90., 170., 250., 330.])
a[:,1]
array([ 10., 90., 170., 250., 330.])
a[2:4,3:6]
a[::2,0::3]
array([[ 0., 30., 60.], [160., 190., 220.], [320., 350., 380.]])
There are two types of advanced indexing
: integer and boolean. Advanced indexing always returns a copy of the data (contrast with basic slicing that returns a view).
i = np.array([0, 2, 3])
j = np.array([1, 1, 5])
b = a[i,j]
b
array([ 10., 170., 290.])
i = np.array([True, False, False, True, True])
c = a[:,3][i] # b is not a view of a
c
array([ 30., 270., 350.])
a=np.arange(1, 5)
a-0.5
a*2.5
a/2
(returns float)a//2
(rounded away from +infinity)a%3
a**2
np.log(a)
-5//2.
-3.0
a=np.array([1., 2., 3.])
b=np.array([2., -1., 0.])
a-b
a*b
a/b
a//b
a%b
a**2
a.dot(b)
or a@b
np.tensordot(a,b,axes=0)
np.tensordot(a,b,axes=0)
array([[ 2., -1., 0.], [ 4., -2., 0.], [ 6., -3., 0.]])
Operations like sum, mean, max, ... are array reduction operations.
We call these reduction operations because operations like sum have the effect of slicing out n-D arrays from the input array, and reducing these n-D arrays to (n-m)-D arrays or scalars.
a=10.*np.arange(0, 40).reshape((5,8))
a.mean()
a.std()
a.sum()
data.min()
xprofile = data.sum(axis=0)
xprofile
array([1998260, 974663, 999642, ..., 1013854, 982238, 2025375], dtype=uint64)
yprofile = data.sum(axis=1)
np.allclose(xprofile.sum(), yprofile.sum())
True
xprofile.sum()
8939461635
Due to rounding errors, most floating-point numbers end up being slightly imprecise. As long as this imprecision stays small, it can usually be ignored. However, it also means that numbers expected to be equal (e.g. when calculating the same result through different correct methods) often differ slightly, and a simple equality test fails
a=0.15+0.15
b=0.1+0.2
a==b
False
numpy.allclose
: use this function to compare arrays (or scalars)
If the following equation is element-wise True
then allclose
returns True
.
a=np.array([.5, 1., 1.05+0.45])
b=np.array([.5, 1., 1.+0.5])
a==b
np.allclose(a,b)
True
np.allclose?
with fits.open("data/mods1r.20171125.0026.fits") as hdulist:
image = hdulist[0].data.astype(float)
NOTE: LUCI and MODS files have been downloaded from i2a archives in Trieste.
quadrants = [image[:3088//2,:8288//2],
image[:3088//2,8288//2:],
image[3088//2:,:8288//2],
image[3088//2:,8288//2:]]
prscans = [image[:3088//2,:48],
image[:3088//2,-48:],
image[3088//2:,:48],
image[3088//2:,-48:]]
For each prescan region, we compute the bias level following this recipe:
Here we compute the bias level for each channel and remove these levels from image
from astropy import stats
for i in range(4):
med = np.median(prscans[i], axis=0)
bias_even = stats.sigma_clip(med[::2], sigma=2.0, maxiters=4)
bias_odd = stats.sigma_clip(med[1::2], sigma=2.0, maxiters=4)
quadrants[i][:, ::2]-=bias_even.mean()
quadrants[i][:, 1::2]-=bias_odd.mean()
fits.PrimaryHDU(data=image).writeto("tmp/bias_sub.fits", overwrite=True)
The astropy.stats
package holds statistical functions or algorithms used in astronomy. While the scipy.stats
and statsmodel packages contains a wide range of statistical tools, they are general-purpose packages and are missing some tools that are particularly useful or specific to astronomy. This package is intended to provide such functionality, but not to replace scipy.stats
if its implementation satisfies astronomers’ needs.
numpy.ma
a.data
array([1, 2, 3])
a.mask
array([False, True, False])
data
operations¶sum
no masked data (default)a.sum()
4
sum
all dataa.data.sum()
6
np.log(a)
masked_array(data=[0.0, --, 1.0986122886681098], mask=[False, True, False], fill_value=1e+20)
mask
operations¶a[~a.mask]
masked_array(data=[1, 3], mask=[False, False], fill_value=999999)
a.mask=[True, False, False]
a[0]=2
a
masked_array(data=[2, 2, 3], mask=[False, False, False], fill_value=999999)