<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Multidimentional data - Matrices and Images
Step1:  Let us work with the matrix
Step2:  numpy matrix multiply uses the dot() function
Step3:  Caution the * will just multiply the matricies on an element-by-element basis
Step4:  Solving system of linear equations
$$
\begin{array}{c}
x + 2y = 4 \
x + y = 3 \
\end{array}
\hspace{2cm}
\left[
\begin{array}{cc}
1 & 2 \
1 & 1 \
\end{array}
\right]
\left[
\begin{array}{c}
x\
y
\end{array}
\right]
=
\left[
\begin{array}{c}
4\
3\ 
\end{array}
\right]
\hspace{2cm}
{\bf A}x = {\bf b}
\hspace{2cm}
\left[
\begin{array}{c}
x\
y
\end{array}
\right]
=
\left[
\begin{array}{cc}
1 & 2 \
1 & 1 \
\end{array}
\right]^{-1}
\left[
\begin{array}{c}
4\
3\ 
\end{array}
\right]
=
\left[
\begin{array}{c}
2\
1
\end{array}
\right]
$$
Step5:  System of 3 equations example (Numpy)
Step6:  System of 3 equations example (SymPy) - Python's Symbolic Math Package
Step7:  SymPy is slower than NumPy
Step8:  Images are just 2-d arrays - imshow will display 2-d arrays as images
Step9:  Read in some data
Step10:  Math on images applies to every value (pixel)
Step11:  Show the image represenation of I with a colorbar
Step12:  Colormap reference
Step13:  WARNING! Common image formats DO NOT preserve dynamic range of original data!!
Common image formats
Step14:  Creating images from math
Step15:  Fancy Image Display
Step16:  Reading in images (imread) - Common Formats
Step17:  Images are just arrays that can be sliced.
For common image formats the origin is the upper left hand corner
Step18:  Simple image manipulation
Step19:  ndimage can do much more
Step20:  You can use masks on images
Step21:  You can add and subtract images
Step22:  The two images above may look the same but they are not! Subtracting the two images reveals the truth.
Step23:  FITS Tables - An astronomical example
Stellar spectra data from the ESO Library of Stellar Spectra
Step24:  Pseudocolor - All color astronomy images are fake.
Color images are composed of three 2-d images
Step25:  We just want to read in one of the three channels
Step26:  Need to create a blank 3-d array to hold all of the images
Step27:  Fill the array with the filtered images
<ASSISTANT_TASK:>
Python Code:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy import linalg
plt.style.use('ggplot')
plt.rc('axes', grid=False)   # turn off the background grid for images
Explanation: Multidimentional data - Matrices and Images
End of explanation
my_matrix = np.array([[1,2],[1,1]])
print(my_matrix.shape)
print(my_matrix)
my_matrix_transposed = np.transpose(my_matrix)
print(my_matrix_transposed)
my_matrix_inverse = linalg.inv(my_matrix)
print(my_matrix_inverse)
Explanation: Let us work with the matrix:
$
\left[
\begin{array}{cc}
1 & 2 \
1 & 1
\end{array}
\right]
$
End of explanation
my_matrix_inverse.dot(my_matrix)
Explanation: numpy matrix multiply uses the dot() function:
End of explanation
my_matrix_inverse * my_matrix_inverse
Explanation: Caution the * will just multiply the matricies on an element-by-element basis:
End of explanation
A = np.array([[1,2],[1,1]])
print(A)
b = np.array([[4],[3]])
print(b)
# Solve by inverting A and then mulitply by b
linalg.inv(A).dot(b) 
# Cleaner looking
linalg.solve(A,b)
Explanation: Solving system of linear equations
$$
\begin{array}{c}
x + 2y = 4 \
x + y = 3 \
\end{array}
\hspace{2cm}
\left[
\begin{array}{cc}
1 & 2 \
1 & 1 \
\end{array}
\right]
\left[
\begin{array}{c}
x\
y
\end{array}
\right]
=
\left[
\begin{array}{c}
4\
3\ 
\end{array}
\right]
\hspace{2cm}
{\bf A}x = {\bf b}
\hspace{2cm}
\left[
\begin{array}{c}
x\
y
\end{array}
\right]
=
\left[
\begin{array}{cc}
1 & 2 \
1 & 1 \
\end{array}
\right]^{-1}
\left[
\begin{array}{c}
4\
3\ 
\end{array}
\right]
=
\left[
\begin{array}{c}
2\
1
\end{array}
\right]
$$
End of explanation
A = np.array([[1,3,5],[2,5,1],[2,3,8]])
b = np.array([[10],[8],[3]])
print(linalg.inv(A))
print(linalg.solve(A,b))
Explanation: System of 3 equations example (Numpy):
$$
\begin{array}{c}
x + 3y + 5z = 10 \
2x + 5y + z = 8 \
2x + 3y + 8z = 3 \
\end{array}
\hspace{3cm}
\left[
\begin{array}{ccc}
1 & 3 & 5 \
2 & 5 & 1 \
2 & 3 & 8 
\end{array}
\right]
\left[
\begin{array}{c}
x\
y\
z 
\end{array}
\right]
=
\left[
\begin{array}{c}
10\
8\
3 
\end{array}
\right]
$$
End of explanation
import sympy as sym
AA = sym.Matrix([[1,3,5],[2,5,1],[2,3,8]])
bb = sym.Matrix([[10],[8],[3]])
print(AA**-1)
print(AA**-1 * bb)
Explanation: System of 3 equations example (SymPy) - Python's Symbolic Math Package
End of explanation
%timeit AA**-1 * bb
%timeit linalg.solve(A,b)
Explanation: SymPy is slower than NumPy
End of explanation
print(A)
plt.imshow(A, interpolation='nearest', cmap=plt.cm.Blues);
Explanation: Images are just 2-d arrays - imshow will display 2-d arrays as images
End of explanation
I = np.load("test_data.npy")    # load in a saved numpy array
I.ndim, I.shape, I.dtype
print("The minimum value of the array I is {0:.2f}".format(I.min()))
print("The maximum value of the array I is {0:.2f}".format(I.max()))
print("The mean value of the array I is {0:.2f}".format(I.mean()))
print("The standard deviation of the array I is {0:.2f}".format(I.std()))
#flatten() collapses n-dimentional data into 1-d
plt.hist(I.flatten(),bins=30);
Explanation: Read in some data
End of explanation
II = I + 8
print("The minimum value of the array II is {0:.2f}".format(II.min()))
print("The maximum value of the array II is {0:.2f}".format(II.max()))
print("The mean value of the array II is {0:.2f}".format(II.mean()))
print("The standard deviation of the array II is {0:.2f}".format(II.std()))
Explanation: Math on images applies to every value (pixel)
End of explanation
plt.imshow(I, cmap=plt.cm.gray)
plt.colorbar();
Explanation: Show the image represenation of I with a colorbar
End of explanation
fig, ax = plt.subplots(1,5,sharey=True)
fig.set_size_inches(12,6)
fig.tight_layout()
ax[0].imshow(I, cmap=plt.cm.viridis)
ax[0].set_xlabel('viridis')
ax[1].imshow(I, cmap=plt.cm.hot)
ax[1].set_xlabel('hot')
ax[2].imshow(I, cmap=plt.cm.magma)
ax[2].set_xlabel('magma')
ax[3].imshow(I, cmap=plt.cm.spectral)
ax[3].set_xlabel('spectral')
ax[4].imshow(I, cmap=plt.cm.gray)
ax[4].set_xlabel('gray')
Explanation: Colormap reference: http://matplotlib.org/examples/color/colormaps_reference.html
End of explanation
plt.imsave('Splash.png', I, cmap=plt.cm.gray)     # Write the array I to a PNG file
Ipng = plt.imread('Splash.png')                   # Read in the PNG file
print("The original data has a min = {0:.2f} and a max = {1:.2f}".format(I.min(), I.max()))
print("The PNG file has a min = {0:.2f} and a max = {1:.2f}".format(Ipng.min(), Ipng.max()))
Explanation: WARNING! Common image formats DO NOT preserve dynamic range of original data!!
Common image formats: jpg, gif, png, tiff
Common image formats will re-scale your data values to [0:1]
Common image formats are NOT suitable for scientific data!
End of explanation
X = np.linspace(-5, 5, 500)
Y = np.linspace(-5, 5, 500)
X, Y = np.meshgrid(X, Y)     # turns two 1-d arrays (X, Y) into one 2-d grid
Z = np.sqrt(X**2+Y**2)+np.sin(X**2+Y**2)
Z.min(), Z.max(), Z.mean()
Explanation: Creating images from math
End of explanation
from matplotlib.colors import LightSource
ls = LightSource(azdeg=0,altdeg=40)
shadedfig = ls.shade(Z,plt.cm.copper)
fig, ax = plt.subplots(1,3)
fig.set_size_inches(12,6)
fig.tight_layout()
ax[0].imshow(shadedfig)
contlevels = [1,2,Z.mean()]
ax[1].axis('equal')
ax[1].contour(Z,contlevels)
ax[2].imshow(shadedfig)
ax[2].contour(Z,contlevels);
Explanation: Fancy Image Display
End of explanation
I2 = plt.imread('doctor5.png')
print("The image I2 has a shape [height,width] of {0}".format(I2.shape))
print("The image I2 is made up of data of type {0}".format(I2.dtype))
print("The image I2 has a maximum value of {0}".format(I2.max()))
print("The image I2 has a minimum value of {0}".format(I2.min()))
plt.imshow(I2,cmap=plt.cm.gray);
Explanation: Reading in images (imread) - Common Formats
End of explanation
fig, ax = plt.subplots(1,4)
fig.set_size_inches(12,6)
fig.tight_layout()
# You can show just slices of the image - Rememeber: The origin is the upper left corner
ax[0].imshow(I2, cmap=plt.cm.gray)
ax[0].set_xlabel('Original')
ax[1].imshow(I2[0:300,0:100], cmap=plt.cm.gray)
ax[1].set_xlabel('[0:300,0:100]')                 # 300 rows, 100 columns
ax[2].imshow(I2[:,0:100], cmap=plt.cm.gray)       # ":" = whole range
ax[2].set_xlabel('[:,0:100]')                     # all rows, 100 columns
ax[3].imshow(I2[:,::-1], cmap=plt.cm.gray);
ax[3].set_xlabel('[:,::-1]')                      # reverse the columns
fig, ax = plt.subplots(1,2)
fig.set_size_inches(12,6)
fig.tight_layout()
CutLine = 300
ax[0].imshow(I2, cmap=plt.cm.gray)
ax[0].hlines(CutLine, 0, 194, color='b', linewidth=3)
ax[1].plot(I2[CutLine,:], color='b', linewidth=3)
ax[1].set_xlabel("X Value")
ax[1].set_ylabel("Pixel Value")
Explanation: Images are just arrays that can be sliced.
For common image formats the origin is the upper left hand corner
End of explanation
from scipy import ndimage
fig, ax = plt.subplots(1,5)
fig.set_size_inches(14,6)
fig.tight_layout()
ax[0].imshow(I2, cmap=plt.cm.gray)
I3 = ndimage.rotate(I2,45,cval=0.75)               # cval is the value to set pixels outside of image
ax[1].imshow(I3, cmap=plt.cm.gray)                 # Rotate and reshape
I4 = ndimage.rotate(I2,45,reshape=False,cval=0.75) # Rotate and do not reshape
ax[2].imshow(I4, cmap=plt.cm.gray)
I5 = ndimage.shift(I2,(10,30),cval=0.75)           # Shift image      
ax[3].imshow(I5, cmap=plt.cm.gray)
I6 = ndimage.gaussian_filter(I2,5)                # Blur image
ax[4].imshow(I6, cmap=plt.cm.gray);
Explanation: Simple image manipulation
End of explanation
import astropy.io.fits as fits
x = fits.open('bsg01.fits')
x.info()
x[0].header
xd = x[0].data
print("The image x has a shape [height,width] of {0}".format(xd.shape))
print("The image x is made up of data of type {0}".format(xd.dtype))
print("The image x has a maximum value of {0}".format(xd.max()))
print("The image x has a minimum value of {0}".format(xd.min()))
fig, ax = plt.subplots(1,2)
fig.set_size_inches(12,6)
fig.tight_layout()
ax[0].imshow(xd,cmap=plt.cm.gray)
ax[1].hist(xd.flatten(),bins=20);
Explanation: ndimage can do much more: http://scipy-lectures.github.io/advanced/image_processing/
FITS file (Flexible Image Transport System) - Standard Astro File Format
FITS format preserves dynamic range of data
FITS format can include lists, tables, images, and combunations of different types of data
End of explanation
CopyData = np.copy(xd)
CutOff = 40
mask = np.where(CopyData > CutOff)
CopyData[mask] = 50                   # You can not just throw data away, you have to set it to something.
fig, ax = plt.subplots(1,2)
fig.set_size_inches(12,6)
fig.tight_layout()
ax[0].imshow(CopyData,cmap=plt.cm.gray)
ax[1].hist(CopyData.flatten(),bins=20);
Explanation: You can use masks on images
End of explanation
fig, ax = plt.subplots(1,2)
fig.set_size_inches(12,6)
fig.tight_layout()
ax[0].imshow(xd, cmap=plt.cm.gray)
# Open another file 'bsg02.fits'
y = fits.open('bsg02.fits')
yd = y[0].data
ax[1].imshow(yd, cmap=plt.cm.gray);
Explanation: You can add and subtract images
End of explanation
fig, ax = plt.subplots(1,3)
fig.set_size_inches(12,6)
fig.tight_layout()
ax[0].imshow(xd, cmap=plt.cm.gray)
ax[1].imshow(yd, cmap=plt.cm.gray)
z = xd - yd                          # Subtract the images pixel by pixel
ax[2].imshow(z, cmap=plt.cm.gray);
Explanation: The two images above may look the same but they are not! Subtracting the two images reveals the truth.
End of explanation
S = fits.open('SolarSpectra.fits')
S.info()
Data = S[0].data
Head = S[0].header
Head
# The FITS header has the information to make an array of wavelengths
Start = Head['CRVAL1']
Number = Head['NAXIS1']
Delta  = Head['CDELT1']
End = Start + (Number * Delta)
Wavelength = np.arange(Start,End,Delta)
fig, ax = plt.subplots(2,1)
fig.set_size_inches(11,8.5)
fig.tight_layout()
# Full spectra
ax[0].plot(Wavelength, Data, color='b')
ax[0].set_ylabel("Flux")
ax[0].set_xlabel("Wavelength [angstroms]")
# Just the visible range with the hydrogen Balmer lines
ax[1].set_xlim(4000,7000)
ax[1].set_ylim(0.6,1.2)
ax[1].plot(Wavelength, Data, color='b')
ax[1].set_ylabel("Flux")
ax[1].set_xlabel("Wavelength [angstroms]")
H_Balmer = [6563,4861,4341,4102,3970,3889,3835,3646]
ax[1].vlines(H_Balmer,0,2, color='r', linewidth=3, alpha = 0.25)
Explanation: FITS Tables - An astronomical example
Stellar spectra data from the ESO Library of Stellar Spectra
End of explanation
redfilter = plt.imread('sphereR.jpg')
redfilter.shape,redfilter.dtype
Explanation: Pseudocolor - All color astronomy images are fake.
Color images are composed of three 2-d images: <img src="images/Layers.png" width="150">
JPG images are 3-d, even grayscale images
End of explanation
redfilter = plt.imread('sphereR.jpg')[:,:,0]
redfilter.shape,redfilter.dtype
plt.imshow(redfilter,cmap=plt.cm.gray);
greenfilter = plt.imread('sphereG.jpg')[:,:,0]
bluefilter = plt.imread('sphereB.jpg')[:,:,0]
fig, ax = plt.subplots(1,3)
fig.set_size_inches(12,3)
fig.tight_layout()
ax[0].set_title("Red Filter")
ax[1].set_title("Green Filter")
ax[2].set_title("Blue Filter")
ax[0].imshow(redfilter,cmap=plt.cm.gray)
ax[1].imshow(greenfilter,cmap=plt.cm.gray)
ax[2].imshow(bluefilter,cmap=plt.cm.gray);
Explanation: We just want to read in one of the three channels
End of explanation
rgb = np.zeros((480,640,3),dtype='uint8')
print(rgb.shape, rgb.dtype)
plt.imshow(rgb,cmap=plt.cm.gray);
Explanation: Need to create a blank 3-d array to hold all of the images
End of explanation
rgb[:,:,0] = redfilter
rgb[:,:,1] = greenfilter
rgb[:,:,2] = bluefilter
fig, ax = plt.subplots(1,4)
fig.set_size_inches(14,3)
fig.tight_layout()
ax[0].set_title("Red Filter")
ax[1].set_title("Green Filter")
ax[2].set_title("Blue Filter")
ax[3].set_title("All Filters Stacked")
ax[0].imshow(redfilter,cmap=plt.cm.gray)
ax[1].imshow(greenfilter,cmap=plt.cm.gray)
ax[2].imshow(bluefilter,cmap=plt.cm.gray)
ax[3].imshow(rgb,cmap=plt.cm.gray);
print("The image rgb has a shape [height,width] of {0}".format(rgb.shape))
print("The image rgb is made up of data of type {0}".format(rgb.dtype))
print("The image rgb has a maximum value of {0}".format(rgb.max()))
print("The image rgb has a minimum value of {0}".format(rgb.min()))
rgb[:,:,0] = redfilter * 1.5
plt.imshow(rgb)
Explanation: Fill the array with the filtered images
End of explanation
<END_TASK>