In this tutorial, we will learn to classify spectral data using the Principal Components Analysis (PCA) method.
Objectives
After completing this tutorial, you will be able to:
- Classify spectral remote sensing data using Principal Components Analysis.
Install Python Packages
- numpy
- gdal
- matplotlib
- matplotlib.pyplot
Download Data
Additional Materials
This tutorial was prepared in conjunction with a presentation on spectral classification that can be downloaded.
Set up
First, we'll start by setting up the necessary environment.
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from scipy import linalg
from scipy import io
from mpl_toolkits.mplot3d import Axes3D
def PlotSpectraAndMean(Spectra, Wv, fignum):
### Spectra is NBands x NSamps
mu = np.mean(Spectra, axis=1)
print(np.shape(mu))
plt.figure(fignum)
plt.plot(Wv, Spectra, 'c')
plt.plot(Wv, mu, 'r')
plt.show()
return mu
Now we can load the spectra.
filename = 'OSBSTinyIm.mat'
ImDict = io.loadmat(filename)
OSBSTinyIm = ImDict['OSBSTinyIm']
TinySize = np.shape(OSBSTinyIm)
NRows = TinySize[0]
NCols = TinySize[1]
NBands = TinySize[2]
print('{0:4d} {1:4d} {2:4d}'.format(NRows, NCols, NBands))
62 194 346
Now we can extract wavelengths.
### LOAD WAVELENGTHS WITH WATER BANDS ###
### AND BAD BEGINNING AND ENDING BANDS REMOVED ###
Wv = io.loadmat("NEONWvsNBB")
Wv = Wv['NEONWvsNBB']
print(np.shape(Wv))
plt.figure(1)
plt.plot(range(346), Wv)
plt.show()
(346, 1)
Let's load indices for Red, Green, and Blue for NEON hyperspectral data.
### HAVE TO SUBTRACT AN OFFSET BECAUSE OF BAD BAND ###
### REMOVAL AND 0-BASED Python vs 1-Based MATLAB ###
Offset = 7
### LOAD & PRINT THE INDICES FOR THE COLORS ###
### AND DIG THEM OUT OF MANY LAYERS OF ARRAYS ###
NEONColors = io.loadmat('NEONColors.mat')
NEONRed = NEONColors['NEONRed']
NEONGreen = NEONColors['NEONGreen']
NEONBlue = NEONColors['NEONBlue']
NEONNir = NEONColors['NEONNir']
NEONRed = NEONRed[0][0]-Offset
NEONGreen = NEONGreen[0][0]-Offset
NEONBlue = NEONBlue[0][0]-Offset
NEONNir = NEONNir[0][0]-Offset
print('Indices: {0:4d} {1:4d} {2:4d} {3:4d}'.format(NEONRed, NEONGreen, NEONBlue, NEONNir))
### CONVERT THE INDICES TO WAVELENGTHS ###
NEONRedWv = Wv[NEONRed][0]
NEONGreenWv = Wv[NEONGreen][0]
NEONBlueWv = Wv[NEONBlue][0]
NEONNirWv = Wv[NEONNir][0]
print('Wavelengths: {0:4d} {1:4d} {2:4d} {3:4d}'.format(NEONRedWv, NEONGreenWv, NEONBlueWv, NEONNirWv))
Indices: 47 25 6 119
Wavelengths: 645 535 440 1005
Now we can make a color image and display it.
RGBIm = OSBSTinyIm[:, :, [NEONRed, NEONGreen, NEONBlue]]
RGBIm = np.sqrt(RGBIm)
plt.figure(2)
plt.imshow(RGBIm)
plt.show()
Now let's turn the image into a sequence of vectors so we can use matrix algebra.
### HAVE TO TAKE INTO ACCOUNT DIFFERENCES BETWEEN Python AND Matlab ###
### Python USES THE C PROGRAMMING LANGUAGE ORDERING ###
### MATLAB USERS THE FORTRAN PROGRAMMING LANGUAGE ORDERING ###
### Python WOULD RESHAPE BY REFERENCE AND MATLAB BY VALUE ###
### THEREFORE, WE NEED TO COPY THE VALUES EXPLICITLY ###
TinyVecs = OSBSTinyIm.reshape(NRows*NCols, NBands, order='F').copy()
### MATLAB TREATS THE ROWS AS DATA SAMPLES ###
### np TREATS THE COLS AS DATA SAMPLES ###
TinyVecs = np.transpose(TinyVecs)
NSamps = np.shape(TinyVecs)[1]
np.shape(TinyVecs)
(346, 12028)
### EXERCISE
SpecIndices = range(1000, 2000, 100)
SomeSpectra = TinyVecs[:, range(1000, 2000, 100)]
mymu = PlotSpectraAndMean(SomeSpectra, Wv, 3)
(346,)
np.shape(mymu)
(346,)
Let's plot some spectra.
### Indices of Spectra to Try ###
### SpecIndices = range(0, 1000, 100) ###
SpecIndices = range(1000, 2000, 100)
SomeSpectra = TinyVecs[:, range(1000, 2000, 100)]
plt.figure(3)
plt.plot(Wv, SomeSpectra)
plt.xlabel('Wavelengths in nm')
plt.ylabel('Reflectance')
plt.show()
Compute the average spectrum and plot it.
mu = np.mean(TinyVecs, axis=1)
plt.figure(4)
plt.plot(Wv, SomeSpectra, 'c')
plt.plot(Wv, mu, 'k')
plt.xlabel('Wavelengths in nm')
plt.ylabel('Reflectance')
plt.show()
Now we want to subtract the mean from every sample.
np.shape(mu)
(346,)
TinyVecsZ = np.zeros((NBands, NSamps))
for n in range(NSamps):
TinyVecsZ[range(NBands),n]= TinyVecs[(range(NBands), n)]-mu
muz = np.mean(TinyVecsZ, axis=1)
plt.figure(5)
plt.plot(Wv, muz, 'k')
#plt.ylim(-1,1)
plt.show()
Let's calculate the covariance.
C = np.cov(TinyVecs)
np.shape(C)
(346, 346)
We can look at some of the values but there are too many to look at them all. We can also view C as an image.
plt.figure(6)
plt.imshow(C)
plt.show()
# PRINT OUT SOME "AMPLIFIED" COVARIANCE VALUES %%%
for cn in range(0, 50,5):
w = int(Wv[cn])
if cn==0:
print(" ", end=" ")
else:
print('{0:5d}'.format(w), end=" ")
print('\n')
for rn in range(5, 50, 5):
w = int(Wv[rn])
print('{0:5d}'.format(w), end=" ")
for cn in range(5,50,5):
CovVal = int(100000*C[rn, rn])
print('{0:5d}'.format(CovVal), end=" ")
print('\n')
#print(round(100000*C[NEONBlue, NEONNir]))
#print(round(100000*C[NEONGreen, NEONNir]))
#print(round(100000*C[NEONRed, NEONNir]))
#print(round(100000*C[NEONGreen, NEONRed]))
435 460 485 510 535 560 585 610 635
435 166 166 166 166 166 166 166 166 166
460 198 198 198 198 198 198 198 198 198
485 226 226 226 226 226 226 226 226 226
510 245 245 245 245 245 245 245 245 245
535 289 289 289 289 289 289 289 289 289
560 360 360 360 360 360 360 360 360 360
585 428 428 428 428 428 428 428 428 428
610 479 479 479 479 479 479 479 479 479
635 544 544 544 544 544 544 544 544 544
Notice that there are no negative values. Why?
What if we normalize the vectors to have Norm 1 (a common strategy).
Norms = np.sqrt(np.sum(TinyVecs*TinyVecs, axis=0))
plt.figure(7)
plt.plot(Norms)
### Too many Norms. How do we fix?
plt.show()
np.shape(Norms)
(12028,)
np.shape(TinyVecs)
(346, 12028)
### Allocate Memory
TinyVecsNorm = np.zeros((NBands, NSamps))
for samp in range(NSamps):
NormSamp = Norms[samp]
for band in range(NBands):
TinyVecsNorm[band, samp] = TinyVecs[band,samp]/NormSamp
Norms1 = np.sqrt(np.sum(TinyVecsNorm*TinyVecsNorm, axis=0))
plt.figure(7)
plt.plot(Norms1)
plt.show()
BigNorm = np.max(Norms1)
LitNorm = np.min(Norms1)
print('{0:4f} {1:4f}'.format(BigNorm, LitNorm))
### Too many Norms. How do we fix?
1.000000 1.000000
Challenge: Plotting Spectra with Mean Function
Turn the script for plotting spectra and their mean above into a function.
### EXERCISE
SpecIndices = range(1000, 2000, 100)
SomeSpectraNorm = TinyVecsNorm[:, range(1000, 2000, 100)]
MuNorm = PlotSpectraAndMean(SomeSpectraNorm, Wv, 3)
(346,)
CNorm = np.cov(TinyVecsNorm)
plt.figure()
plt.imshow(CNorm)
plt.show()
# PRINT OUT SOME "AMPLIFIED" COVARIANCE VALUES %%%
for cn in range(0, 50,5):
w = int(Wv[cn])
if cn==0:
print(" ", end=" ")
else:
print('{0:5d}'.format(w), end=" ")
print('\n')
for rn in range(5, 50, 5):
w = int(Wv[rn])
print('{0:5d}'.format(w), end=" ")
for cn in range(5,50,5):
CovVal = int(10000000*CNorm[rn, rn])
print('{0:5d}'.format(CovVal), end=" ")
print('\n')
435 460 485 510 535 560 585 610 635
435 253 253 253 253 253 253 253 253 253
460 269 269 269 269 269 269 269 269 269
485 289 289 289 289 289 289 289 289 289
510 284 284 284 284 284 284 284 284 284
535 263 263 263 263 263 263 263 263 263
560 303 303 303 303 303 303 303 303 303
585 397 397 397 397 397 397 397 397 397
610 455 455 455 455 455 455 455 455 455
635 523 523 523 523 523 523 523 523 523
print(np.shape(TinyVecs))
print(NEONNir)
print(NEONRed)
NIRVals = TinyVecs[NEONNir, range(NSamps)]
RedVals = TinyVecs[NEONRed, range(NSamps)]
NDVIVals = (NIRVals-RedVals)/(NIRVals+RedVals)
np.shape(NDVIVals)
NDVIIm = np.reshape(NDVIVals,(NRows, NCols),order='F')
(346, 12028)
119
47
print(np.shape(NDVIIm))
(62, 194)
plt.figure()
plt.hist(NDVIVals)
plt.show()
HiNDVI = NDVIIm*(NDVIIm>0.8)
plt.figure()
plt.imshow(HiNDVI)
plt.show()
# plt.figure()
# plt.plot(nonzero(NDVIVals>0.8))
# plt.show()
VegIndices = np.nonzero(NDVIVals>0.8)
# print(VegIndices[0])
print(np.shape(VegIndices))
# print(np.shape(TinyVecs))
VegSpectra = TinyVecs[:, VegIndices[0]]
print(np.shape(VegSpectra))
(1, 3136)
(346, 3136)
CVeg = np.cov(VegSpectra)
plt.figure(9)
plt.imshow?
plt.imshow(CVeg,extent=(np.amin(Wv), np.amax(Wv),np.amax(Wv), np.amin(Wv)))
plt.colorbar()
plt.show()
PCA
OK, let's do PCA
Recall that TinyVecs is the mean-subtracted version of the original spectra.
C = np.cov(TinyVecs)
D,V = linalg.eig(C)
D = D.real
print(np.shape(D))
print(np.shape(V))
print(TinyVecs.shape)
print(V[0,0])
(346,)
(346, 346)
(346, 12028)
-0.0124822914573
plt.figure(10)
print(D.shape)
DiagD = np.diag(D)
print(D.shape)
plt.plot(DiagD)
#Exercise
#plt.plot(D[range(10)])
#plt.plot(D[range(10, 30, 10)])
plt.show()
(346,)
(346,)
TinyVecsPCA = np.dot(V.T, TinyVecs)
PCACovar = np.cov(TinyVecsPCA)
D,V = linalg.eig(C)
D = D.real
print(D.shape)
print(PCACovar.shape)
for r in range(10):
print('{0:5f} {1:5f}'.format(D[r], PCACovar[r,r]))
print()
for r in range(10):
for c in range(10):
NextVal = int(10000*PCACovar[r,c])
print('{0:5d}'.format(NextVal), end=" ")
print('\n')
# #Delta = np.sum(np.sum((PCACovar-D), axis=0), axis=0)
# print(Delta)
# plt.figure(11)
# plt.plot(np.diag(PCACovar))
# plt.show()
(346,)
(346, 346)
4.911250 4.911250
0.653442 0.653442
0.083168 0.083168
0.013260 0.013260
0.004038 0.004038
0.003201 0.003201
0.001390 0.001390
0.000931 0.000931
0.000705 0.000705
0.000561 0.000561
49112 0 0 0 0 0 0 0 0 0
0 6534 0 0 0 0 0 0 0 0
0 0 831 0 0 0 0 0 0 0
0 0 0 132 0 0 0 0 0 0
0 0 0 0 40 0 0 0 0 0
0 0 0 0 0 32 0 0 0 0
0 0 0 0 0 0 13 0 0 0
0 0 0 0 0 0 0 9 0 0
0 0 0 0 0 0 0 0 7 0
0 0 0 0 0 0 0 0 0 5
Notice that the values on the diagonal are the variances of each coordinate in the PCA transformed data. They drop off rapidly which is why one can reduce dimensionality by discarding components that have low variance. Also, notice that the diagonal matrix D produce by diagonalizing the covariance of x is the covariance of y = PCA(x).
If the data are Gaussian, then the coordinates of y are uncorrelated and independent. If not, then only uncorrelated.
Let's pull out the first 3 dimensions and plot them.
%matplotlib notebook
fig = plt.figure(13)
ax = fig.add_subplot(111, projection='3d')
ax.scatter(TinyVecsPCA[0,range(NSamps)],TinyVecsPCA[1,range(NSamps)],TinyVecsPCA[2,range(NSamps)], marker='o')
plt.show()
We can also display principal components as images
for coord in range(3):
P1 = TinyVecsPCA[coord, :]
PCAIm = np.reshape(P1, (NRows, NCols), order='F')
plt.figure(14+coord)
plt.imshow(np.abs(PCAIm))
plt.colorbar()
plt.show()
পল গ্যাডার
এই টিউটোরিয়ালে, আমরা প্রিন্সিপাল কম্পোনেন্ট অ্যানালাইসিস (পিসিএ) পদ্ধতি ব্যবহার করে বর্ণালী ডেটা শ্রেণিবদ্ধ করতে শিখব।
উদ্দেশ্য
এই টিউটোরিয়ালটি শেষ করার পরে, আপনি সক্ষম হবেন:
- অধ্যক্ষ উপাদান বিশ্লেষণ ব্যবহার করে বর্ণালী দূরবর্তী সেন্সিং ডেটা শ্রেণিবদ্ধ করুন।
পাইথন প্যাকেজ ইনস্টল করুন
- numpy
- gdal
- matplotlib
- matplotlib.pyplot
ডাউনলোড ডেটা
অতিরিক্ত উপকরণ
এই টিউটোরিয়ালটি বর্ণালী শ্রেণিবিন্যাসের ডাউনলোডের জন্য উপস্থাপনার সাথে তৈরি করা হয়েছিল।
সেট আপ করুন
প্রথমত, আমরা প্রয়োজনীয় পরিবেশ স্থাপন করে শুরু করব।
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from scipy import linalg
from scipy import io
from mpl_toolkits.mplot3d import Axes3D
def PlotSpectraAndMean(Spectra, Wv, fignum):
### Spectra is NBands x NSamps
mu = np.mean(Spectra, axis=1)
print(np.shape(mu))
plt.figure(fignum)
plt.plot(Wv, Spectra, 'c')
plt.plot(Wv, mu, 'r')
plt.show()
return mu
এখন আমরা স্পেকট্রা লোড করতে পারি।
filename = 'OSBSTinyIm.mat'
ImDict = io.loadmat(filename)
OSBSTinyIm = ImDict['OSBSTinyIm']
TinySize = np.shape(OSBSTinyIm)
NRows = TinySize[0]
NCols = TinySize[1]
NBands = TinySize[2]
print('{0:4d} {1:4d} {2:4d}'.format(NRows, NCols, NBands))
62 194 346
এখন আমরা তরঙ্গদৈর্ঘ্য বের করতে পারি।
### LOAD WAVELENGTHS WITH WATER BANDS ###
### AND BAD BEGINNING AND ENDING BANDS REMOVED ###
Wv = io.loadmat("NEONWvsNBB")
Wv = Wv['NEONWvsNBB']
print(np.shape(Wv))
plt.figure(1)
plt.plot(range(346), Wv)
plt.show()
(346, 1)
নিওন হাইপারস্পেকট্রাল ডেটার জন্য লাল, সবুজ এবং নীলের সূচকগুলি লোড করা যাক।
### HAVE TO SUBTRACT AN OFFSET BECAUSE OF BAD BAND ###
### REMOVAL AND 0-BASED Python vs 1-Based MATLAB ###
Offset = 7
### LOAD & PRINT THE INDICES FOR THE COLORS ###
### AND DIG THEM OUT OF MANY LAYERS OF ARRAYS ###
NEONColors = io.loadmat('NEONColors.mat')
NEONRed = NEONColors['NEONRed']
NEONGreen = NEONColors['NEONGreen']
NEONBlue = NEONColors['NEONBlue']
NEONNir = NEONColors['NEONNir']
NEONRed = NEONRed[0][0]-Offset
NEONGreen = NEONGreen[0][0]-Offset
NEONBlue = NEONBlue[0][0]-Offset
NEONNir = NEONNir[0][0]-Offset
print('Indices: {0:4d} {1:4d} {2:4d} {3:4d}'.format(NEONRed, NEONGreen, NEONBlue, NEONNir))
### CONVERT THE INDICES TO WAVELENGTHS ###
NEONRedWv = Wv[NEONRed][0]
NEONGreenWv = Wv[NEONGreen][0]
NEONBlueWv = Wv[NEONBlue][0]
NEONNirWv = Wv[NEONNir][0]
print('Wavelengths: {0:4d} {1:4d} {2:4d} {3:4d}'.format(NEONRedWv, NEONGreenWv, NEONBlueWv, NEONNirWv))
Indices: 47 25 6 119
Wavelengths: 645 535 440 1005
এখন আমরা একটি রঙিন চিত্র তৈরি করতে পারি এবং এটি প্রদর্শন করতে পারি।
RGBIm = OSBSTinyIm[:, :, [NEONRed, NEONGreen, NEONBlue]]
RGBIm = np.sqrt(RGBIm)
plt.figure(2)
plt.imshow(RGBIm)
plt.show()
এখন আসুন চিত্রটিকে ভেক্টরের অনুক্রমে পরিণত করি যাতে আমরা ম্যাট্রিক্স বীজগণিতটি ব্যবহার করতে পারি।
### HAVE TO TAKE INTO ACCOUNT DIFFERENCES BETWEEN Python AND Matlab ###
### Python USES THE C PROGRAMMING LANGUAGE ORDERING ###
### MATLAB USERS THE FORTRAN PROGRAMMING LANGUAGE ORDERING ###
### Python WOULD RESHAPE BY REFERENCE AND MATLAB BY VALUE ###
### THEREFORE, WE NEED TO COPY THE VALUES EXPLICITLY ###
TinyVecs = OSBSTinyIm.reshape(NRows*NCols, NBands, order='F').copy()
### MATLAB TREATS THE ROWS AS DATA SAMPLES ###
### np TREATS THE COLS AS DATA SAMPLES ###
TinyVecs = np.transpose(TinyVecs)
NSamps = np.shape(TinyVecs)[1]
np.shape(TinyVecs)
(346, 12028)
### EXERCISE
SpecIndices = range(1000, 2000, 100)
SomeSpectra = TinyVecs[:, range(1000, 2000, 100)]
mymu = PlotSpectraAndMean(SomeSpectra, Wv, 3)
(346,)
np.shape(mymu)
(346,)
আসুন কিছু বর্ণালী প্লট করা যাক।
### Indices of Spectra to Try ###
### SpecIndices = range(0, 1000, 100) ###
SpecIndices = range(1000, 2000, 100)
SomeSpectra = TinyVecs[:, range(1000, 2000, 100)]
plt.figure(3)
plt.plot(Wv, SomeSpectra)
plt.xlabel('Wavelengths in nm')
plt.ylabel('Reflectance')
plt.show()
গড় বর্ণালী গণনা করুন এবং এটি প্লট করুন।
mu = np.mean(TinyVecs, axis=1)
plt.figure(4)
plt.plot(Wv, SomeSpectra, 'c')
plt.plot(Wv, mu, 'k')
plt.xlabel('Wavelengths in nm')
plt.ylabel('Reflectance')
plt.show()
এখন আমরা প্রতিটি নমুনা থেকে গড় বিয়োগ করতে চাই।
np.shape(mu)
(346,)
TinyVecsZ = np.zeros((NBands, NSamps))
for n in range(NSamps):
TinyVecsZ[range(NBands),n]= TinyVecs[(range(NBands), n)]-mu
muz = np.mean(TinyVecsZ, axis=1)
plt.figure(5)
plt.plot(Wv, muz, 'k')
#plt.ylim(-1,1)
plt.show()
আসল গণনা করা যাক।
C = np.cov(TinyVecs)
np.shape(C)
(346, 346)
আমরা কয়েকটি মান দেখতে পারি তবে সেগুলি দেখতে অনেকগুলি রয়েছে। আমরা সি হিসাবে একটি চিত্র হিসাবে দেখতে পারেন।
plt.figure(6)
plt.imshow(C)
plt.show()
# PRINT OUT SOME "AMPLIFIED" COVARIANCE VALUES %%%
for cn in range(0, 50,5):
w = int(Wv[cn])
if cn==0:
print(" ", end=" ")
else:
print('{0:5d}'.format(w), end=" ")
print('\n')
for rn in range(5, 50, 5):
w = int(Wv[rn])
print('{0:5d}'.format(w), end=" ")
for cn in range(5,50,5):
CovVal = int(100000*C[rn, rn])
print('{0:5d}'.format(CovVal), end=" ")
print('\n')
#print(round(100000*C[NEONBlue, NEONNir]))
#print(round(100000*C[NEONGreen, NEONNir]))
#print(round(100000*C[NEONRed, NEONNir]))
#print(round(100000*C[NEONGreen, NEONRed]))
435 460 485 510 535 560 585 610 635
435 166 166 166 166 166 166 166 166 166
460 198 198 198 198 198 198 198 198 198
485 226 226 226 226 226 226 226 226 226
510 245 245 245 245 245 245 245 245 245
535 289 289 289 289 289 289 289 289 289
560 360 360 360 360 360 360 360 360 360
585 428 428 428 428 428 428 428 428 428
610 479 479 479 479 479 479 479 479 479
635 544 544 544 544 544 544 544 544 544
লক্ষ্য করুন যে কোনও নেতিবাচক মান নেই। কেন?
যদি আমরা ভেক্টরগুলিকে নরম 1 (সাধারণ কৌশল) রাখার জন্য স্বাভাবিক করি।
Norms = np.sqrt(np.sum(TinyVecs*TinyVecs, axis=0))
plt.figure(7)
plt.plot(Norms)
### Too many Norms. How do we fix?
plt.show()
np.shape(Norms)
(12028,)
np.shape(TinyVecs)
(346, 12028)
### Allocate Memory
TinyVecsNorm = np.zeros((NBands, NSamps))
for samp in range(NSamps):
NormSamp = Norms[samp]
for band in range(NBands):
TinyVecsNorm[band, samp] = TinyVecs[band,samp]/NormSamp
Norms1 = np.sqrt(np.sum(TinyVecsNorm*TinyVecsNorm, axis=0))
plt.figure(7)
plt.plot(Norms1)
plt.show()
BigNorm = np.max(Norms1)
LitNorm = np.min(Norms1)
print('{0:4f} {1:4f}'.format(BigNorm, LitNorm))
### Too many Norms. How do we fix?
1.000000 1.000000
প্রতিদ্বন্দ্বিতা: স্পেনট্রা প্লেনিং মানে মিডিয়া ফাংশন
স্পেকট্রা প্লট করার জন্য স্ক্রিপ্টটি এবং উপরে তাদের ফাংশন রূপান্তর করুন।
### EXERCISE
SpecIndices = range(1000, 2000, 100)
SomeSpectraNorm = TinyVecsNorm[:, range(1000, 2000, 100)]
MuNorm = PlotSpectraAndMean(SomeSpectraNorm, Wv, 3)
(346,)
CNorm = np.cov(TinyVecsNorm)
plt.figure()
plt.imshow(CNorm)
plt.show()
# PRINT OUT SOME "AMPLIFIED" COVARIANCE VALUES %%%
for cn in range(0, 50,5):
w = int(Wv[cn])
if cn==0:
print(" ", end=" ")
else:
print('{0:5d}'.format(w), end=" ")
print('\n')
for rn in range(5, 50, 5):
w = int(Wv[rn])
print('{0:5d}'.format(w), end=" ")
for cn in range(5,50,5):
CovVal = int(10000000*CNorm[rn, rn])
print('{0:5d}'.format(CovVal), end=" ")
print('\n')
435 460 485 510 535 560 585 610 635
435 253 253 253 253 253 253 253 253 253
460 269 269 269 269 269 269 269 269 269
485 289 289 289 289 289 289 289 289 289
510 284 284 284 284 284 284 284 284 284
535 263 263 263 263 263 263 263 263 263
560 303 303 303 303 303 303 303 303 303
585 397 397 397 397 397 397 397 397 397
610 455 455 455 455 455 455 455 455 455
635 523 523 523 523 523 523 523 523 523
print(np.shape(TinyVecs))
print(NEONNir)
print(NEONRed)
NIRVals = TinyVecs[NEONNir, range(NSamps)]
RedVals = TinyVecs[NEONRed, range(NSamps)]
NDVIVals = (NIRVals-RedVals)/(NIRVals+RedVals)
np.shape(NDVIVals)
NDVIIm = np.reshape(NDVIVals,(NRows, NCols),order='F')
(346, 12028)
119
47
print(np.shape(NDVIIm))
(62, 194)
plt.figure()
plt.hist(NDVIVals)
plt.show()
HiNDVI = NDVIIm*(NDVIIm>0.8)
plt.figure()
plt.imshow(HiNDVI)
plt.show()
# plt.figure()
# plt.plot(nonzero(NDVIVals>0.8))
# plt.show()
VegIndices = np.nonzero(NDVIVals>0.8)
# print(VegIndices[0])
print(np.shape(VegIndices))
# print(np.shape(TinyVecs))
VegSpectra = TinyVecs[:, VegIndices[0]]
print(np.shape(VegSpectra))
(1, 3136)
(346, 3136)
CVeg = np.cov(VegSpectra)
plt.figure(9)
plt.imshow?
plt.imshow(CVeg,extent=(np.amin(Wv), np.amax(Wv),np.amax(Wv), np.amin(Wv)))
plt.colorbar()
plt.show()
পিসিএ
ঠিক আছে, পিসিএ করি
মনে করুন যে টিনিভেকস হ'ল মূল বর্ণালির গড়-বিয়োগ সংস্করণ।
C = np.cov(TinyVecs)
D,V = linalg.eig(C)
D = D.real
print(np.shape(D))
print(np.shape(V))
print(TinyVecs.shape)
print(V[0,0])
(346,)
(346, 346)
(346, 12028)
-0.0124822914573
plt.figure(10)
print(D.shape)
DiagD = np.diag(D)
print(D.shape)
plt.plot(DiagD)
#Exercise
#plt.plot(D[range(10)])
#plt.plot(D[range(10, 30, 10)])
plt.show()
(346,)
(346,)
TinyVecsPCA = np.dot(V.T, TinyVecs)
PCACovar = np.cov(TinyVecsPCA)
D,V = linalg.eig(C)
D = D.real
print(D.shape)
print(PCACovar.shape)
for r in range(10):
print('{0:5f} {1:5f}'.format(D[r], PCACovar[r,r]))
print()
for r in range(10):
for c in range(10):
NextVal = int(10000*PCACovar[r,c])
print('{0:5d}'.format(NextVal), end=" ")
print('\n')
# #Delta = np.sum(np.sum((PCACovar-D), axis=0), axis=0)
# print(Delta)
# plt.figure(11)
# plt.plot(np.diag(PCACovar))
# plt.show()
(346,)
(346, 346)
4.911250 4.911250
0.653442 0.653442
0.083168 0.083168
0.013260 0.013260
0.004038 0.004038
0.003201 0.003201
0.001390 0.001390
0.000931 0.000931
0.000705 0.000705
0.000561 0.000561
49112 0 0 0 0 0 0 0 0 0
0 6534 0 0 0 0 0 0 0 0
0 0 831 0 0 0 0 0 0 0
0 0 0 132 0 0 0 0 0 0
0 0 0 0 40 0 0 0 0 0
0 0 0 0 0 32 0 0 0 0
0 0 0 0 0 0 13 0 0 0
0 0 0 0 0 0 0 9 0 0
0 0 0 0 0 0 0 0 7 0
0 0 0 0 0 0 0 0 0 5
লক্ষ্য করুন যে তির্যকটির মানগুলি পিসিএ রূপান্তরিত তথ্যের প্রতিটি সমন্বয়ের বৈকল্পিক। তারা দ্রুত বিসর্জন দেয় যে কারণে কম বৈকল্পিক রয়েছে এমন উপাদানগুলি ত্যাগ করে কেউ মাত্রিকতা হ্রাস করতে পারে। এছাড়াও, লক্ষ করুন যে এক্স এর সমবায়ু দ্বারা তির্যক ম্যাট্রিক্স ডি উত্পাদন করে y = PCA (x) এর সমবায়তা।
যদি ডেটা গাউসিয়ান হয় তবে y এর স্থানাঙ্কগুলি অসংলগ্ন এবং স্বতন্ত্র। যদি তা না হয় তবে কেবল অনিয়ন্ত্রিত।
আসুন প্রথম 3 টি মাত্রা বের করুন এবং সেগুলি প্লট করুন।
%matplotlib notebook
fig = plt.figure(13)
ax = fig.add_subplot(111, projection='3d')
ax.scatter(TinyVecsPCA[0,range(NSamps)],TinyVecsPCA[1,range(NSamps)],TinyVecsPCA[2,range(NSamps)], marker='o')
plt.show()
আমরা চিত্র হিসাবে মূল উপাদানগুলি প্রদর্শন করতে পারি
for coord in range(3):
P1 = TinyVecsPCA[coord, :]
PCAIm = np.reshape(P1, (NRows, NCols), order='F')
plt.figure(14+coord)
plt.imshow(np.abs(PCAIm))
plt.colorbar()
plt.show()
চ্যালেঞ্জ: পিসিএ শ্রেণিবদ্ধকরণ ফাংশন
স্পেকট্রার পিসিএ গণনা করে প্রথম ফাংশন প্লট করে এমন একটি ফাংশন লিখুন।
0 comments:
Post a Comment