-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathraster_pca.R
219 lines (191 loc) · 10.1 KB
/
raster_pca.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
library(raster)
library(RStoolbox)
library(readr)
library(base)
library(doSNOW)
library(doParallel)
library(parallel)
img.dir = paste0(data.dir,"/S2A_MSIL1C_20170101T082332_N0204_R121_T34JFP_20170101T084543.SAFE/GRANULE/L1C_T34JFP_A007983_20170101T084543/IMG_DATA")
img.dir3 = paste0(data.dir,"/S2A_MSIL1C_20170312T082001_N0204_R121_T34JFP_20170312T084235.SAFE/GRANULE/L1C_T34JFP_A008984_20170312T084235/IMG_DATA")
# img.dir6 = paste0(data.dir,"/S2A_MSIL1C_20170131T082151_N0204_R121_T34JEP_20170131T084118.SAFE/GRANULE/L1C_T34JFP_A008412_20170131T084118/IMG_DATA")
#
# img.dir4 = paste0(data.dir,"/S2B_MSIL1C_20170715T081609_N0205_R121_T34JEP_20170715T084650.SAFE/GRANULE/L1C_T34JEP_A001863_20170715T084650/IMG_DATA")
# img.dir7 = paste0(data.dir,"/S2B_MSIL1C_20170804T081559_N0205_R121_T34JEP_20170804T084631.SAFE/GRANULE/L1C_T34JEP_A002149_20170804T084631/IMG_DATA")
img.dir8 = paste0(data.dir,"/S2A_MSIL1C_20170819T082011_N0205_R121_T34JFP_20170819T084427.SAFE/GRANULE/L1C_T34JFP_A011272_20170819T084427/IMG_DATA")
img.dir5 = paste0(data.dir,"/S2A_MSIL1C_20170531T082011_N0205_R121_T34JFP_20170531T084246.SAFE/GRANULE/L1C_T34JFP_A010128_20170531T084246/IMG_DATA")
# img.dir2 = paste0(data.dir,"/S2A_MSIL1C_20170210T082051_N0204_R121_T34JEP_20170210T083752.SAFE/GRANULE/L1C_T34JEP_A008555_20170210T083752/IMG_DATA")
#
cat("Mapping PCA 1.........")
#### Load Raster
band1 = raster::raster(paste0(img.dir,'/T34JFP_20170101T082332_B01.jp2'))
band2 = raster::raster(paste0(img.dir,'/T34JFP_20170101T082332_B02.jp2'))
band3 = raster::raster(paste0(img.dir,'/T34JFP_20170101T082332_B03.jp2'))
band4 = raster::raster(paste0(img.dir,'/T34JFP_20170101T082332_B04.jp2'))
band5 = raster::raster(paste0(img.dir,'/T34JFP_20170101T082332_B05.jp2'))
band6 = raster::raster(paste0(img.dir,'/T34JFP_20170101T082332_B06.jp2'))
band7 = raster::raster(paste0(img.dir,'/T34JFP_20170101T082332_B07.jp2'))
band8 = raster::raster(paste0(img.dir,'/T34JFP_20170101T082332_B08.jp2'))
band9 = raster::raster(paste0(img.dir,'/T34JFP_20170101T082332_B09.jp2'))
band10 = raster::raster(paste0(img.dir,'/T34JFP_20170101T082332_B10.jp2'))
band11 = raster::raster(paste0(img.dir,'/T34JFP_20170101T082332_B11.jp2'))
#bandTC = raster::raster(paste0(img.dir,'/T34JFP_20170101T082332_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,band3,band4,band8))
stack3 = rasterPCA(raster::stack(band5,band6,band7,band11))
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:4))
pc3 = extract(stack3$map,train_points)
colnames(pc3) = c(paste0("stack3_PC",1:4))
trainpca2 = cbind(pc,pc2,pc3)
### test
pc = extract(stack1$map,test_points)
colnames(pc) = c(paste0("stack1_PC",1:3))
pc2 = extract(stack2$map,test_points)
colnames(pc2) = c(paste0("stack2_PC",1:4))
pc3 = extract(stack3$map,test_points)
colnames(pc3) = c(paste0("stack3_PC",1:4))
testpca2 = cbind(pc,pc2,pc3)
cat("Mapping PCA 2.........")
###### 2017-05-31
band1 = raster::raster(paste0(img.dir5,'/T34JFP_20170531T082011_B01.jp2'))
band2 = raster::raster(paste0(img.dir5,'/T34JFP_20170531T082011_B02.jp2'))
band3 = raster::raster(paste0(img.dir5,'/T34JFP_20170531T082011_B03.jp2'))
band4 = raster::raster(paste0(img.dir5,'/T34JFP_20170531T082011_B04.jp2'))
band5 = raster::raster(paste0(img.dir5,'/T34JFP_20170531T082011_B05.jp2'))
band6 = raster::raster(paste0(img.dir5,'/T34JFP_20170531T082011_B06.jp2'))
band7 = raster::raster(paste0(img.dir5,'/T34JFP_20170531T082011_B07.jp2'))
band8 = raster::raster(paste0(img.dir5,'/T34JFP_20170531T082011_B08.jp2'))
band9 = raster::raster(paste0(img.dir5,'/T34JFP_20170531T082011_B09.jp2'))
band10 = raster::raster(paste0(img.dir5,'/T34JFP_20170531T082011_B10.jp2'))
band11 = raster::raster(paste0(img.dir5,'/T34JFP_20170531T082011_B11.jp2'))
##### RASTER PCA
stack1_5 = rasterPCA(raster::stack(band1,band9,band10))
stack2_5 = rasterPCA(raster::stack(band2,band3,band4))
stack3_5 = rasterPCA(raster::stack(band5,band6,band7))
pc = extract(stack1_5$map,train_points)
colnames(pc) = c(paste0("stack1_5_PC",1:3))
pc2 = extract(stack2_5$map,train_points)
colnames(pc2) = c(paste0("stack2_5_PC",1:3))
pc3 = extract(stack3_5$map,train_points)
colnames(pc3) = c(paste0("stack3_5_PC",1:3))
trainpca2 = cbind(trainpca2,pc,pc2,pc3)
### test
pc = extract(stack1_5$map,test_points)
colnames(pc) = c(paste0("stack1_5_PC",1:3))
pc2 = extract(stack2_5$map,test_points)
colnames(pc2) = c(paste0("stack2_5_PC",1:3))
pc3 = extract(stack3_5$map,test_points)
colnames(pc3) = c(paste0("stack3_5_PC",1:3))
testpca2 = cbind(testpca2,pc,pc2,pc3)
# cat("Mapping PCA 3...........")
# ###### 2017-08-19
# band1 = raster::raster(paste0(img.dir8,'/T34JEP_20170819T082011_B01.jp2'))
# band2 = raster::raster(paste0(img.dir8,'/T34JEP_20170819T082011_B02.jp2'))
# band3 = raster::raster(paste0(img.dir8,'/T34JEP_20170819T082011_B03.jp2'))
# band4 = raster::raster(paste0(img.dir8,'/T34JEP_20170819T082011_B04.jp2'))
# band5 = raster::raster(paste0(img.dir8,'/T34JEP_20170819T082011_B05.jp2'))
# band6 = raster::raster(paste0(img.dir8,'/T34JEP_20170819T082011_B06.jp2'))
# band7 = raster::raster(paste0(img.dir8,'/T34JEP_20170819T082011_B07.jp2'))
# band8 = raster::raster(paste0(img.dir8,'/T34JEP_20170819T082011_B08.jp2'))
# band9 = raster::raster(paste0(img.dir8,'/T34JEP_20170819T082011_B09.jp2'))
# band10 = raster::raster(paste0(img.dir8,'/T34JEP_20170819T082011_B10.jp2'))
# band11 = raster::raster(paste0(img.dir8,'/T34JEP_20170819T082011_B11.jp2'))
#
# ##### RASTER PCA
# stack1_8 = rasterPCA(raster::stack(band1,band9,band10))
# stack2_8 = rasterPCA(raster::stack(band2,band3,band4))
# stack3_8 = rasterPCA(raster::stack(band5,band6,band7))
#
# pc = extract(stack1_8$map,train_points)
# colnames(pc) = c(paste0("stack1_8_PC",1:3))
# pc2 = extract(stack2_8$map,train_points)
# colnames(pc2) = c(paste0("stack2_8_PC",1:3))
# pc3 = extract(stack3_8$map,train_points)
# colnames(pc3) = c(paste0("stack3_8_PC",1:3))
# trainpca = cbind(trainpca,pc,pc2,pc3)
# ### test
# pc = extract(stack1_8$map,test_points)
# colnames(pc) = c(paste0("stack1_8_PC",1:3))
# pc2 = extract(stack2_8$map,test_points)
# colnames(pc2) = c(paste0("stack2_8_PC",1:3))
# pc3 = extract(stack3_8$map,test_points)
# colnames(pc3) = c(paste0("stack3_8_PC",1:3))
# testpca = cbind(testpca,pc,pc2,pc3)
cat("Mapping PCA 4.........")
###### 2017-03-12
band1 = raster::raster(paste0(img.dir3,'/T34JFP_20170312T082001_B01.jp2'))
band2 = raster::raster(paste0(img.dir3,'/T34JFP_20170312T082001_B02.jp2'))
band3 = raster::raster(paste0(img.dir3,'/T34JFP_20170312T082001_B03.jp2'))
band4 = raster::raster(paste0(img.dir3,'/T34JFP_20170312T082001_B04.jp2'))
band5 = raster::raster(paste0(img.dir3,'/T34JFP_20170312T082001_B05.jp2'))
band6 = raster::raster(paste0(img.dir3,'/T34JFP_20170312T082001_B06.jp2'))
band7 = raster::raster(paste0(img.dir3,'/T34JFP_20170312T082001_B07.jp2'))
band8 = raster::raster(paste0(img.dir3,'/T34JFP_20170312T082001_B08.jp2'))
band9 = raster::raster(paste0(img.dir3,'/T34JFP_20170312T082001_B09.jp2'))
band10 = raster::raster(paste0(img.dir3,'/T34JFP_20170312T082001_B10.jp2'))
band11 = raster::raster(paste0(img.dir3,'/T34JFP_20170312T082001_B11.jp2'))
##### RASTER PCA
stack1_12 = rasterPCA(raster::stack(band1,band9,band10))
stack2_12 = rasterPCA(raster::stack(band2,band3,band4))
stack3_12 = rasterPCA(raster::stack(band5,band6,band7))
pc = extract(stack1_12$map,train_points)
colnames(pc) = c(paste0("stack1_12_PC",1:3))
pc2 = extract(stack2_12$map,train_points)
colnames(pc2) = c(paste0("stack2_12_PC",1:3))
pc3 = extract(stack3_12$map,train_points)
colnames(pc3) = c(paste0("stack3_12_PC",1:3))
trainpca2 = cbind(trainpca2,pc,pc2,pc3)
### test
pc = extract(stack1_12$map,test_points)
colnames(pc) = c(paste0("stack1_12_PC",1:3))
pc2 = extract(stack2_12$map,test_points)
colnames(pc2) = c(paste0("stack2_12_PC",1:3))
pc3 = extract(stack3_12$map,test_points)
colnames(pc3) = c(paste0("stack3_12_PC",1:3))
testpca2 = cbind(testpca2,pc,pc2,pc3)
# cat("Mapping PCA 5...........")
# ###### 2017-08-04
# band1 = raster::raster(paste0(img.dir7,'/T34JEP_20170804T081559_B01.jp2'))
# band2 = raster::raster(paste0(img.dir7,'/T34JEP_20170804T081559_B02.jp2'))
# band3 = raster::raster(paste0(img.dir7,'/T34JEP_20170804T081559_B03.jp2'))
# band4 = raster::raster(paste0(img.dir7,'/T34JEP_20170804T081559_B04.jp2'))
# band5 = raster::raster(paste0(img.dir7,'/T34JEP_20170804T081559_B05.jp2'))
# band6 = raster::raster(paste0(img.dir7,'/T34JEP_20170804T081559_B06.jp2'))
# band7 = raster::raster(paste0(img.dir7,'/T34JEP_20170804T081559_B07.jp2'))
# band8 = raster::raster(paste0(img.dir7,'/T34JEP_20170804T081559_B08.jp2'))
# band9 = raster::raster(paste0(img.dir7,'/T34JEP_20170804T081559_B09.jp2'))
# band10 = raster::raster(paste0(img.dir7,'/T34JEP_20170804T081559_B10.jp2'))
# band11 = raster::raster(paste0(img.dir7,'/T34JEP_20170804T081559_B11.jp2'))
#
# ##### RASTER PCA
# stack1_7 = rasterPCA(raster::stack(band1,band9,band10))
# stack2_7 = rasterPCA(raster::stack(band2,band3,band4))
# stack3_7 = rasterPCA(raster::stack(band5,band6,band7))
#
# pc = extract(stack1_7$map,train_points)
# colnames(pc) = c(paste0("stack1_7_PC",1:3))
# pc2 = extract(stack2_7$map,train_points)
# colnames(pc2) = c(paste0("stack2_7_PC",1:3))
# pc3 = extract(stack3_7$map,train_points)
# colnames(pc3) = c(paste0("stack3_7_PC",1:3))
# trainpca = cbind(trainpca,pc,pc2,pc3)
# ### test
# pc = extract(stack1_7$map,test_points)
# colnames(pc) = c(paste0("stack1_7_PC",1:3))
# pc2 = extract(stack2_7$map,test_points)
# colnames(pc2) = c(paste0("stack2_7_PC",1:3))
# pc3 = extract(stack3_7$map,test_points)
# colnames(pc3) = c(paste0("stack3_7_PC",1:3))
# testpca = cbind(testpca,pc,pc2,pc3)
stopCluster(cl)
Sys.time() - start.time
# save(trainpca, file = paste0(save.files.dir,"/trainpca_data.RData"))
# save(testpca, file = paste0(save.files.dir,"/testpca_data.RData"))