-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathpca2.R
233 lines (201 loc) · 10.3 KB
/
pca2.R
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
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
library(RStoolbox)
library(raster)
library(doSNOW)
library(doParallel)
library(parallel)
# img.dir1 = paste0(data.dir,"/S2A_MSIL1C_20170322T081611_N0204_R121_T34JEP_20170322T084728.SAFE/GRANULE/L1C_T34JEP_A009127_20170322T084728/IMG_DATA")
#
# cat("Mapping PCA 1.........")
# #### Load Raster
# #### 23-03-2017
# band1 = raster::raster(paste0(img.dir1,'/T34JEP_20170322T081611_B01.jp2'))
# band2 = raster::raster(paste0(img.dir1,'/T34JEP_20170322T081611_B02.jp2'))
# band3 = raster::raster(paste0(img.dir1,'/T34JEP_20170322T081611_B03.jp2'))
# band4 = raster::raster(paste0(img.dir1,'/T34JEP_20170322T081611_B04.jp2'))
# band5 = raster::raster(paste0(img.dir1,'/T34JEP_20170322T081611_B05.jp2'))
# band6 = raster::raster(paste0(img.dir1,'/T34JEP_20170322T081611_B06.jp2'))
# band7 = raster::raster(paste0(img.dir1,'/T34JEP_20170322T081611_B07.jp2'))
# band8 = raster::raster(paste0(img.dir1,'/T34JEP_20170322T081611_B08.jp2'))
#band8A = raster::raster(paste0(img.dir1,'/T34JEP_20170322T081611_B8A.jp2'))
# band9 = raster::raster(paste0(img.dir1,'/T34JEP_20170322T081611_B09.jp2'))
# band10 = raster::raster(paste0(img.dir1,'/T34JEP_20170322T081611_B10.jp2'))
# band11 = raster::raster(paste0(img.dir1,'/T34JEP_20170322T081611_B11.jp2'))
# bandTC = raster::raster(paste0(img.dir1,'/T34JEP_20170322T081611_TCI.jp2'))
##### RASTER PCA
# Time the code execution
# start.time <- Sys.time()
# # Create a cluster to work on 10 logical cores.
# cl <- makeCluster(8, type = "PSOCK")
# registerDoParallel(cl)
# #stack1 = rasterPCA(raster::stack(band1,band9,band10))
# stack2 = rasterPCA(raster::stack(band2,band4))
# stack3 = rasterPCA(raster::stack(band5,band6))
#
# # pc = extract(stack1$map,train_points)
# # colnames(pc) = c(paste0("stack1_PC",1:3))
# pc2 = extract(stack2$map,train_points)
# colnames(pc2) = c(paste0("stack2_PC",1:2))
# pc3 = extract(stack3$map,train_points)
# colnames(pc3) = c(paste0("stack3_PC",1:2))
#
# trainpca3 = cbind(pc2,pc3)
#
# ### test
# # pc = extract(stack1$map,test_points)
# # colnames(pc) = c(paste0("stack1_5_PC",1:3))
# pc2 = extract(stack2$map,test_points)
# colnames(pc2) = c(paste0("stack2_PC",1:2))
# pc3 = extract(stack3$map,test_points)
# colnames(pc3) = c(paste0("stack3_PC",1:2))
# testpca3 = cbind(pc2,pc3)
#
#
# img.dir1 = paste0(data.dir,"/S2A_MSIL1C_20170322T081611_N0204_R121_T34JFP_20170322T084728.SAFE/GRANULE/L1C_T34JFP_A009127_20170322T084728/IMG_DATA")
#
# cat("Mapping PCA 2.........")
# ###### 2017-03-22
# band1 = raster::raster(paste0(img.dir1,'/T34JFP_20170322T081611_B01.jp2'))
# band2 = raster::raster(paste0(img.dir1,'/T34JFP_20170322T081611_B02.jp2'))
# band3 = raster::raster(paste0(img.dir1,'/T34JFP_20170322T081611_B03.jp2'))
# band4 = raster::raster(paste0(img.dir1,'/T34JFP_20170322T081611_B04.jp2'))
# band5 = raster::raster(paste0(img.dir1,'/T34JFP_20170322T081611_B05.jp2'))
# band6 = raster::raster(paste0(img.dir1,'/T34JFP_20170322T081611_B06.jp2'))
# band7 = raster::raster(paste0(img.dir1,'/T34JFP_20170322T081611_B07.jp2'))
# band8 = raster::raster(paste0(img.dir1,'/T34JFP_20170322T081611_B08.jp2'))
# #band8A = raster::raster(paste0(img.dir1,'/T34JFP_20170322T081611_B8A.jp2'))
# # band9 = raster::raster(paste0(img.dir1,'/T34JFP_20170322T081611_B09.jp2'))
# # band10 = raster::raster(paste0(img.dir1,'/T34JFP_20170322T081611_B10.jp2'))
# # band11 = raster::raster(paste0(img.dir1,'/T34JFP_20170322T081611_B11.jp2'))
#
# ##### RASTER PCA
#
# #stack1 = rasterPCA(raster::stack(band1,band9,band10))
# stack2 = rasterPCA(raster::stack(band2,band4))
# stack3 = rasterPCA(raster::stack(band5,band6))
#
# # pc = extract(stack1$map,train_points)
# # colnames(pc) = c(paste0("stack1_PC",1:3))
# pc2 = extract(stack2$map,train_points)
# colnames(pc2) = c(paste0("stack2_PC",1:2))
# pc3 = extract(stack3$map,train_points)
# colnames(pc3) = c(paste0("stack3_PC",1:2))
#
# trainpca4 = cbind(pc2,pc3)
#
# # pc = extract(stack1$map,test_points)
# # colnames(pc) = c(paste0("stack1_5_PC",1:3))
# pc2 = extract(stack2$map,test_points)
# colnames(pc2) = c(paste0("stack2_PC",1:2))
# pc3 = extract(stack3$map,test_points)
# colnames(pc3) = c(paste0("stack3_PC",1:2))
# testpca4 = cbind(pc2,pc3)
img.dir6 = paste0(data.dir,"/S2A_MSIL1C_20170131T082151_N0204_R121_T34JFP_20170131T084118.SAFE/GRANULE/L1C_T34JFP_A008412_20170131T084118/IMG_DATA")
###### 2017-01-31
band1 = raster::raster(paste0(img.dir6,'/T34JFP_20170131T082151_B01.jp2'))
band2 = raster::raster(paste0(img.dir6,'/T34JFP_20170131T082151_B02.jp2'))
band3 = raster::raster(paste0(img.dir6,'/T34JFP_20170131T082151_B03.jp2'))
band4 = raster::raster(paste0(img.dir6,'/T34JFP_20170131T082151_B04.jp2'))
band5 = raster::raster(paste0(img.dir6,'/T34JFP_20170131T082151_B05.jp2'))
band6 = raster::raster(paste0(img.dir6,'/T34JFP_20170131T082151_B06.jp2'))
band7 = raster::raster(paste0(img.dir6,'/T34JFP_20170131T082151_B07.jp2'))
band8 = raster::raster(paste0(img.dir6,'/T34JFP_20170131T082151_B08.jp2'))
band8A = raster::raster(paste0(img.dir6,'/T34JFP_20170131T082151_B8A.jp2'))
band9 = raster::raster(paste0(img.dir6,'/T34JFP_20170131T082151_B09.jp2'))
band10 = raster::raster(paste0(img.dir6,'/T34JFP_20170131T082151_B10.jp2'))
band11 = raster::raster(paste0(img.dir6,'/T34JFP_20170131T082151_B11.jp2'))
bandTC = raster::raster(paste0(img.dir6,'/T34JFP_20170131T082151_TCI.jp2'))
##### RASTER PCA
# Time the code execution
start.time <- Sys.time()
# Create a cluster to work on 10 logical cores.
cl <- makeCluster(8, type = "PSOCK")
registerDoParallel(cl)
cat("Mapping PCA 1.........")
#stack1 = rasterPCA(raster::stack(band1,band9,band10))
stack2 = rasterPCA(raster::stack(band2,band3,band4,band8),nComp = 2)
stack3 = rasterPCA(raster::stack(band5,band6,band7,band11),nComp = 2)
# pc = extract(stack1$map,train_points)
# colnames(pc) = c(paste0("stack1_PC",1:3))
pc2 = extract(stack2$map,train_points)
colnames(pc2) = c(paste0("stack2_1_PC",1:2))
pc3 = extract(stack3$map,train_points)
colnames(pc3) = c(paste0("stack3_1_PC",1:2))
trainpca6 = cbind(pc2,pc3)
### test
# pc = extract(stack1$map,test_points)
# colnames(pc) = c(paste0("stack1_5_PC",1:3))
pc2 = extract(stack2$map,test_points)
colnames(pc2) = c(paste0("stack2_1_PC",1:2))
pc3 = extract(stack3$map,test_points)
colnames(pc3) = c(paste0("stack3_1_PC",1:2))
testpca6 = cbind(pc2,pc3)
img.dir2 = paste0(data.dir,"/S2A_MSIL1C_20170210T082051_N0204_R121_T34JFP_20170210T083752.SAFE/GRANULE/L1C_T34JFP_A008555_20170210T083752/IMG_DATA")
#### 10-02-2017
band1 = raster::raster(paste0(img.dir2,'/T34JFP_20170210T082051_B01.jp2'))
band2 = raster::raster(paste0(img.dir2,'/T34JFP_20170210T082051_B02.jp2'))
band3 = raster::raster(paste0(img.dir2,'/T34JFP_20170210T082051_B03.jp2'))
band4 = raster::raster(paste0(img.dir2,'/T34JFP_20170210T082051_B04.jp2'))
band5 = raster::raster(paste0(img.dir2,'/T34JFP_20170210T082051_B05.jp2'))
band6 = raster::raster(paste0(img.dir2,'/T34JFP_20170210T082051_B06.jp2'))
band7 = raster::raster(paste0(img.dir2,'/T34JFP_20170210T082051_B07.jp2'))
band8 = raster::raster(paste0(img.dir2,'/T34JFP_20170210T082051_B08.jp2'))
band8A = raster::raster(paste0(img.dir2,'/T34JFP_20170210T082051_B8A.jp2'))
band9 = raster::raster(paste0(img.dir2,'/T34JFP_20170210T082051_B09.jp2'))
band10 = raster::raster(paste0(img.dir2,'/T34JFP_20170210T082051_B10.jp2'))
band11 = raster::raster(paste0(img.dir2,'/T34JFP_20170210T082051_B11.jp2'))
bandTC = raster::stack(paste0(img.dir2,'/T34JFP_20170210T082051_TCI.jp2'))
cat("Mapping PCA 2.........")
#stack1 = rasterPCA(raster::stack(band1,band9,band10))
stack2 = rasterPCA(raster::stack(band2,band3,band4,band8),nComp = 2)
stack3 = rasterPCA(raster::stack(band5,band6,band7,band11),nComp = 2)
# pc = extract(stack1$map,train_points)
# colnames(pc) = c(paste0("stack1_PC",1:3))
pc2 = extract(stack2$map,train_points)
colnames(pc2) = c(paste0("stack2_2_PC",1:2))
pc3 = extract(stack3$map,train_points)
colnames(pc3) = c(paste0("stack3_2_PC",1:2))
trainpca6 = cbind(trainpca6,pc2,pc3)
### test
# pc = extract(stack1$map,test_points)
# colnames(pc) = c(paste0("stack1_5_PC",1:3))
pc2 = extract(stack2$map,test_points)
colnames(pc2) = c(paste0("stack2_2_PC",1:2))
pc3 = extract(stack3$map,test_points)
colnames(pc3) = c(paste0("stack3_2_PC",1:2))
testpca6 = cbind(testpca6,pc2,pc3)
img.dir66 = paste0(data.dir,"/S2A_MSIL1C_20170620T082011_N0205_R121_T34JFP_20170620T084200.SAFE/GRANULE/L1C_T34JFP_A010414_20170620T084200/IMG_DATA")
#### 20-06-2017
band1 = raster::raster(paste0(img.dir66,'/T34JFP_20170620T082011_B01.jp2'))
band2 = raster::raster(paste0(img.dir66,'/T34JFP_20170620T082011_B02.jp2'))
band3 = raster::raster(paste0(img.dir66,'/T34JFP_20170620T082011_B03.jp2'))
band4 = raster::raster(paste0(img.dir66,'/T34JFP_20170620T082011_B04.jp2'))
band5 = raster::raster(paste0(img.dir66,'/T34JFP_20170620T082011_B05.jp2'))
band6 = raster::raster(paste0(img.dir66,'/T34JFP_20170620T082011_B06.jp2'))
band7 = raster::raster(paste0(img.dir66,'/T34JFP_20170620T082011_B07.jp2'))
band8 = raster::raster(paste0(img.dir66,'/T34JFP_20170620T082011_B08.jp2'))
band8A = raster::raster(paste0(img.dir66,'/T34JFP_20170620T082011_B8A.jp2'))
band9 = raster::raster(paste0(img.dir66,'/T34JFP_20170620T082011_B09.jp2'))
band10 = raster::raster(paste0(img.dir66,'/T34JFP_20170620T082011_B10.jp2'))
band11 = raster::raster(paste0(img.dir66,'/T34JFP_20170620T082011_B11.jp2'))
bandTC = raster::raster(paste0(img.dir66,'/T34JFP_20170620T082011_TCI.jp2'))
cat("Mapping PCA 3.........")
#stack1 = rasterPCA(raster::stack(band1,band9,band10))
stack2 = rasterPCA(raster::stack(band2,band3,band4,band8),nComp = 2)
stack3 = rasterPCA(raster::stack(band5,band6,band7,band11),nComp = 2)
# pc = extract(stack1$map,train_points)
# colnames(pc) = c(paste0("stack1_PC",1:3))
pc2 = extract(stack2$map,train_points)
colnames(pc2) = c(paste0("stack2_6_PC",1:2))
pc3 = extract(stack3$map,train_points)
colnames(pc3) = c(paste0("stack3_6_PC",1:2))
trainpca6 = cbind(trainpca6,pc2,pc3)
### test
# pc = extract(stack1$map,test_points)
# colnames(pc) = c(paste0("stack1_5_PC",1:3))
pc2 = extract(stack2$map,test_points)
colnames(pc2) = c(paste0("stack2_6_PC",1:2))
pc3 = extract(stack3$map,test_points)
colnames(pc3) = c(paste0("stack3_6_PC",1:2))
testpca6 = cbind(testpca6,pc2,pc3)
stopCluster(cl)
Sys.time() - start.time