-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathdata details.R
1038 lines (853 loc) · 45.5 KB
/
data details.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
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
options(warn = -1)
library(plyr)
library(tidyverse)
library(raster)
library(sp)
library(sf)
library(rgdal)
library(rgeos)
library(velox)
library(celestial)
library(caret)
library(fastICA)
library(SOAR)
library(RStoolbox)
library(spdep)
source("utils.R")
#####
path.dir = getwd()
data.dir = paste0(path.dir,"/Data")
img.dir = paste0(data.dir,"/S2A_MSIL1C_20170101T082332_N0204_R121_T34JEP_20170101T084543.SAFE/GRANULE/L1C_T34JEP_A007983_20170101T084543/IMG_DATA")
img.dir3 = paste0(data.dir,"/S2A_MSIL1C_20170312T082001_N0204_R121_T34JEP_20170312T084235.SAFE/GRANULE/L1C_T34JEP_A008984_20170312T084235/IMG_DATA")
img.dir6 = paste0(data.dir,"/S2A_MSIL1C_20170131T082151_N0204_R121_T34JEP_20170131T084118.SAFE/GRANULE/L1C_T34JEP_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.dir42 = paste0(data.dir,"/S2A_MSIL1C_20170710T082011_N0205_R121_T34JEP_20170710T084244.SAFE/GRANULE/L1C_T34JEP_A010700_20170710T084244/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_T34JEP_20170819T084427.SAFE/GRANULE/L1C_T34JEP_A011272_20170819T084427/IMG_DATA")
img.dir5 = paste0(data.dir,"/S2A_MSIL1C_20170531T082011_N0205_R121_T34JEP_20170531T084246.SAFE/GRANULE/L1C_T34JEP_A010128_20170531T084246/IMG_DATA")
img.dir2 = paste0(data.dir,"/S2A_MSIL1C_20170210T082051_N0204_R121_T34JEP_20170210T083752.SAFE/GRANULE/L1C_T34JEP_A008555_20170210T083752/IMG_DATA")
img.dir66 = paste0(data.dir,"/S2A_MSIL1C_20170620T082011_N0205_R121_T34JEP_20170620T084200.SAFE/GRANULE/L1C_T34JEP_A010414_20170620T084200/IMG_DATA")
img.dir1 = paste0(data.dir,"/S2A_MSIL1C_20170322T081611_N0204_R121_T34JEP_20170322T084728.SAFE/GRANULE/L1C_T34JEP_A009127_20170322T084728/IMG_DATA")
#######
#dir.create(paste0(path.dir,"/tmp"))
save.files.dir = paste0(path.dir,"/tmp")
## sf
# shp = readOGR(paste0(data.dir,"/train2-fixup/new.shp"))
#
# #### sp
# sh = readOGR(paste0(data.dir,"/train/train.shp"),stringsAsFactors = F,dropNULLGeometries = F)
#######
shp2 = readOGR(paste0(data.dir,"/train/train.shp"),stringsAsFactors = F)
sh_trans <- spTransform(shp2, CRS("+proj=longlat +ellps=WGS84 +datum=WGS84"))
df = cbind(shp2 %>% as.data.frame(),coordinates(sh_trans))
colnames(df) = c("Field_Id","Area","Subregion","Crop_Id_Ne","lon","lat")
m = df %>% group_by(Crop_Id_Ne) %>% summarise(mlo =mean(lon),lm = mean(lat))
test = readOGR(paste0(data.dir,"/test/test.shp"),stringsAsFactors = F)
test_trans = spTransform(test,CRS("+proj=longlat +ellps=WGS84 +datum=WGS84"))
test = cbind(test %>% as.data.frame(),coordinates(test_trans))
colnames(test) = c("Field_Id","Area","Subregion","lon","lat")
label = as.numeric(df$Crop_Id_Ne)
# library(CatEncoders)
# label_enc= LabelEncoder.fit(as.factor(label))
# label2 = transform(label_enc,as.factor(label))
train.id = df$Field_Id
test.id = test$Field_Id
df = df %>% within(rm("Crop_Id_Ne","Field_Id"))
test = test %>% within(rm("Field_Id"))
df = rbind(df,test)
########
#@@@@ TODO re-arrange by field id
########
df = df %>%
mutate(
#lat_deg = as.numeric(deg2dms(lat)[,1]),
#lon_deg = as.numeric(deg2dms(lon)[,1]),
latmins = as.numeric(deg2dms(lat)[,2])*60,
lonmins = as.numeric(deg2dms(lon)[,2])*60
) %>%
add_count(latmins) %>%
rename("lat_cnt" = n) %>%
add_count(lonmins) %>%
rename("lon_cnt" = n) %>%
add_count(latmins,lonmins) %>%
rename("lat_lon_cnt" = n) %>%
##############
## add freq encoding
#################
# group_by(lonmins,latmins) %>%
# mutate(m_area = median(Area),
# min = min(Area),
# max = max(Area)) %>%
#ungroup() %>%
# group_by(Subregion) %>%
# mutate(mean_area = median(Area),
# min_area = min(Area),
# max_area = max(Area))
dr.dat = cbind(lon=df$lon,lat=df$lat) %>% as.data.frame()
transform = preProcess(dr.dat, method = c("center","scale","pca"))
pc = predict(transform, dr.dat) %>% as.data.frame()
df$lat_pca = pc$PC1
df$lon_pca = pc$PC2
############
transform = preProcess(dr.dat, method = c("center","scale","ica"),n.comp=2)
ic = predict(transform, dr.dat) %>% as.data.frame()
df$lat_ica = ic$ICA1
df$lon_ica = ic$ICA2
###########
transform = preProcess(dr.dat, method = c("center","scale"))
pc = predict(transform, dr.dat) %>% as.data.frame()
df$lat_scale = pc$lat
df$lon_scale = pc$lon
df$lat1 = cos(df$lat_scale)*cos(df$lon_scale)
df$lat2 = cos(df$lat_scale)*sin(df$lon_scale)
#df$lat3 = sin(df$lat_scale)
#df$rot45_x = 0.707 * df$lat_scale + 0.707*df$lon_scale
#df$rot45_y = 0.707 * df$lon_scale - 0.707*df$lat_scale
df$rot30_x = 0.866 * df$lat_scale + 0.5*df$lon_scale
df$rot30_y = 0.866 * df$lon_scale - 0.5*df$lat_scale
df$rot60_x = 0.5 * df$lat_scale + 0.866*df$lon_scale
df$rot60_y = 0.5 * df$lon_scale - 0.866*df$lat_scale
df$lat = NULL
df$lon = NULL
########
##
#######
df_train = df[1:length(train.id),]
df_test = df[(length(train.id)+1):nrow(df),]
#########
df.dat = cbind(lon = df$lon,lat = df$lat)
sumsq = NULL
for (i in 1:15) {
set.seed(1234)
sumsq[i] = sum(kmeans(j,centers = i, iter.max = 1000,
algorithm = "Lloyd")$withinss)
}
plot(1:15,sumsq,type= "b")
###
set.seed(1234)
kmns = kmeans(dr.dat,8,nstart = 17,iter.max = 1000,
algorithm = "Lloyd",trace = T)
cnts = kmns$centers
kmeans.distance = NULL
for (i in 1:nrow(cnts)) {
kmeans.distance = cbind(kmeans.distance, sqrt(colSums((t(dr.dat)-unlist(cnts[i,]))^2)))
}
table(kmns$cluster)
save(kmeans.distance, file = paste(save.files.dir,"/kmeans_features.RData"))
# ndvi <- VI(, 8, 4) band 8 - band4
###################
###### EXTRACT VALUES FROM IMAGES
#################
train_points = SpatialPoints(sh_trans,CRS("+proj=longlat +ellps=WGS84 +datum=WGS84"))
test_points = SpatialPoints(test_trans,CRS("+proj=longlat +ellps=WGS84 +datum=WGS84"))
####
#filenames = paste0(img.dir,"/T34JEP_20170101T082332_B0",1:9,".jp2")
#### Load Raster
band1 = raster::raster(paste0(img.dir,'/T34JEP_20170101T082332_B01.jp2'))
band2 = raster::raster(paste0(img.dir,'/T34JEP_20170101T082332_B02.jp2'))
band3 = raster::raster(paste0(img.dir,'/T34JEP_20170101T082332_B03.jp2'))
band4 = raster::raster(paste0(img.dir,'/T34JEP_20170101T082332_B04.jp2'))
band5 = raster::raster(paste0(img.dir,'/T34JEP_20170101T082332_B05.jp2'))
band6 = raster::raster(paste0(img.dir,'/T34JEP_20170101T082332_B06.jp2'))
band7 = raster::raster(paste0(img.dir,'/T34JEP_20170101T082332_B07.jp2'))
band8 = raster::raster(paste0(img.dir,'/T34JEP_20170101T082332_B08.jp2'))
band8A = raster::raster(paste0(img.dir,'/T34JEP_20170101T082332_B8A.jp2'))
band9 = raster::raster(paste0(img.dir,'/T34JEP_20170101T082332_B09.jp2'))
band10 = raster::raster(paste0(img.dir,'/T34JEP_20170101T082332_B10.jp2'))
band11 = raster::raster(paste0(img.dir,'/T34JEP_20170101T082332_B11.jp2'))
band12 = raster::raster(paste0(img.dir,'/T34JEP_20170101T082332_B12.jp2'))
bandTC = raster::stack(paste0(img.dir,'/T34JEP_20170101T082332_TCI.jp2'))
# ##### RASTER PCA
# library(doSNOW)
# # 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))
# # Processing is done, stop cluster.
# stopCluster(cl)
# Sys.time() - start.time
####
#ndvi_11_5 = NDVI(band11,band5)
#ndvi_8_2 = NDVI(band8,band2)
## band 1,band 6, band 3
df_train$band1 = raster::extract(band1,train_points)
#df_train$a = extract(band1,train_points,cellnumbers=T)[,1]/10000
df_train$band2 = raster::extract(band2,train_points)
#df_train$a2 = extract(band2,train_points,cellnumbers=T)[,1]/10000
df_train$band3 = raster::extract(band3,train_points)
#df_train$a3 = extract(band3,train_points,cellnumbers=T)[,1]/10000
df_train$band4 = raster::extract(band4,train_points)
#df_train$a4 = extract(band4,train_points,cellnumbers=T)[,1]/10000
df_train$band5 = raster::extract(band5,train_points)
#df_train$a5 = extract(band5,train_points,cellnumbers=T)[,1]/10000
df_train$band6 = raster::extract(band6,train_points)
#df_train$a6 = extract(band6,train_points,cellnumbers=T)[,1]/10000
df_train$band7 = raster::extract(band7,train_points)
#df_train$a7 = extract(band7,train_points,cellnumbers=T)[,1]/10000
df_train$band8 = raster::extract(band8,train_points)
df_train$band8A = raster::extract(band8A,train_points)
#df_train$a8 = extract(band8,train_points,cellnumbers=T)[,1]/10000
df_train$band9 = raster::extract(band9,train_points)
df_train$band10 = raster::extract(band10,train_points)
df_train$band11 = raster::extract(band11,train_points)
df_train$band12 = raster::extract(band12,train_points)
df_train$bandtc = extract(bandTC,train_points)
#df_train$ndvi_1 = extract(ndvi_11_5,train_points)
# 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))
#
# df_train = cbind(df_train,pc2,pc3)
df_train$stack3_PC2=NULL
df_train$stack2_PC1=NULL
df_train$stack2_PC3=NULL
###
df_test$band1 = raster::extract(band1, test_points)
df_test$band2 = raster::extract(band2, test_points)
df_test$band3 = raster::extract(band3, test_points)
df_test$band4 = raster::extract(band4, test_points)
df_test$band5 = raster::extract(band5, test_points)
df_test$band6 = raster::extract(band6, test_points)
df_test$band7 = raster::extract(band7, test_points)
df_test$band8A = raster::extract(band8A, test_points)
df_test$band9 = raster::extract(band9, test_points)
df_test$band10 = raster::extract(band10, test_points)
df_test$band11 = raster::extract(band11, test_points)
df_test$bandtc = extract(bandTC,test_points)
#df_test$ndvi_1 = extract(ndvi_11_5,test_points)
# 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))
# #df_train$ndvi = raster::extract(ndvi, train_points)
# #df_train$band1 = raster::extract(band1,train_points)
# df_test = cbind(df_test,pc2,pc3)
# df_test$stack3_PC2=NULL
# df_test$stack2_PC1=NULL
# df_test$stack2_PC3=NULL
df_train$band1 = ifelse(df_train$band1==1682,1591,df_train$band1)
df_test$band1 = ifelse(df_test$band1==1040,1037,df_test$band1)
df_train$band6 = ifelse(df_train$band6 >4136,4136,df_train$band6) #band6
df_train$band3 = ifelse(df_train$band3 >2262,2262,df_train$band3)
#df_train$ndvi = ifelse(df_train$ndvi >4706, 4706,df_train$ndvi)
df_test$band4 = ifelse(df_test$band4 <483,483,df_test$band4)
#df_test$band2 = ifelse(df_test$band2 )
## read raster files
###### 2017-05-31
band1 = raster::raster(paste0(img.dir5,'/T34JEP_20170531T082011_B01.jp2'))
band2 = raster::raster(paste0(img.dir5,'/T34JEP_20170531T082011_B02.jp2'))
band3 = raster::raster(paste0(img.dir5,'/T34JEP_20170531T082011_B03.jp2'))
band4 = raster::raster(paste0(img.dir5,'/T34JEP_20170531T082011_B04.jp2'))
band5 = raster::raster(paste0(img.dir5,'/T34JEP_20170531T082011_B05.jp2'))
band6 = raster::raster(paste0(img.dir5,'/T34JEP_20170531T082011_B06.jp2'))
band7 = raster::raster(paste0(img.dir5,'/T34JEP_20170531T082011_B07.jp2'))
band8 = raster::raster(paste0(img.dir5,'/T34JEP_20170531T082011_B08.jp2'))
band8A = raster::raster(paste0(img.dir5,'/T34JEP_20170531T082011_B8A.jp2'))
band9 = raster::raster(paste0(img.dir5,'/T34JEP_20170531T082011_B09.jp2'))
band10 = raster::raster(paste0(img.dir5,'/T34JEP_20170531T082011_B10.jp2'))
band11 = raster::raster(paste0(img.dir5,'/T34JEP_20170531T082011_B11.jp2'))
band12 = raster::raster(paste0(img.dir5,'/T34JEP_20170531T082011_B12.jp2'))
bandTC5 = raster::stack(paste0(img.dir5,'/T34JEP_20170531T082011_TCI.jp2'))
##### RASTER PCA
library(doSNOW)
# Time the code execution
start.time <- Sys.time()
# Create a cluster to work on 10 logical cores.
cl <- makeCluster(8, type = "PSOCK")
registerDoParallel(cl)
stack12 = rasterPCA(raster::stack(band1,band9,band10))
stack2 = rasterPCA(raster::stack(band2,band3,band4))
stack3 = rasterPCA(raster::stack(band5,band6,band7))
# Processing is done, stop cluster.
stopCluster(cl)
Sys.time() - start.time
pc = extract(stack1$map,train_points)
colnames(pc) = c(paste0("05_stack1_PC",1:3))
pc2 = extract(stack2$map,train_points)
colnames(pc2) = c(paste0("05_stack2_PC",1:3))
pc3 = extract(stack3$map,train_points)
colnames(pc3) = c(paste0("05_stack3_PC",1:3))
# ndvi_11_5 = NDVI(band11,band5)
df_train$band1_5 = raster::extract(band1,train_points)
df_train$band2_5 = raster::extract(band2,train_points)
df_train$band3_5 = raster::extract(band3,train_points)
df_train$band4_5 = raster::extract(band4,train_points)
df_train$band5_5 = raster::extract(band5,train_points)
df_train$band6_5 = raster::extract(band6,train_points)
df_train$band7_5 = raster::extract(band7,train_points)
df_train$band8_5 = raster::extract(band8,train_points)
df_train$band8A_5 = raster::extract(band8A,train_points)
df_train$band9_5 = raster::extract(band9,train_points)
df_train$band10_5 = raster::extract(band10,train_points)
df_train$band11_5 = raster::extract(band11,train_points)
df_train$band12_5 = raster::extract(band12,train_points)
df_train$bandtc_5 = extract(bandTC,train_points)
df_test$band1_5 = raster::extract(band1,test_points)
df_test$band2_5 = raster::extract(band2,test_points)
df_test$band3_5 = raster::extract(band3,test_points)
df_test$band4_5 = raster::extract(band4,test_points)
df_test$band5_5 = raster::extract(band5,test_points)
df_test$band6_5 = raster::extract(band6,test_points)
df_test$band7_5 = raster::extract(band7,test_points)
df_test$band8_5 = raster::extract(band8,test_points)
df_test$band9_5 = raster::extract(band9,test_points)
df_test$band10_5 = raster::extract(band10,test_points)
df_test$band11_5 = raster::extract(band11,test_points)
df_test$band12_5 = raster::extract(band12,test_points)
df_test$bandtc_5 = extract(bandTC,test_points)
###
raster::pairs(raster::stack(band5,band6,band7,band11))
raster::pairs(raster::stack(band1,band9,band10))
raster::pairs(raster::stack(band2,band3,band4,band8))
###### 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'))
#band8A = raster::raster(paste0(img.dir8,'/T34JEP_20170819T082011_B8A.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'))
band12 = raster::raster(paste0(img.dir8,'/T34JEP_20170819T082011_B12.jp2'))
bandTC = raster::raster(paste0(img.dir8,'/T34JEP_20170819T082011_TCI.jp2'))
#ndvi_11_5 = NDVI(band11,band5)
df_train$band1_8 = raster::extract(band1,train_points)
df_train$band2_8 = raster::extract(band2,train_points)
df_train$band3_8 = raster::extract(band3,train_points)
df_train$band4_8 = raster::extract(band4,train_points)
df_train$band5_8 = raster::extract(band5,train_points)
df_train$band6_8 = raster::extract(band6,train_points)
df_train$band7_8 = raster::extract(band7,train_points)
df_train$band8_8 = raster::extract(band8,train_points)
#df_train$band8A_8 = raster::extract(band8A,train_points)
df_train$band9_8 = raster::extract(band9,train_points)
#df_train$band10_8 = raster::extract(band10,train_points)
df_train$band11_8 = raster::extract(band11,train_points)
df_train$bandtc_8 = extract(bandTC,train_points)
df_test$band1_8 = raster::extract(band1,test_points)
df_test$band2_8 = raster::extract(band2,test_points)
df_test$band3_8 = raster::extract(band3,test_points)
df_test$band4_8 = raster::extract(band4,test_points)
df_test$band5_8 = raster::extract(band5,test_points)
df_test$band6_8 = raster::extract(band6,test_points)
df_test$band7_8 = raster::extract(band7,test_points)
df_test$band8_8 = raster::extract(band8,test_points)
df_test$band9_8 = raster::extract(band9,test_points)
df_test$band10_8 = raster::extract(band10,test_points)
df_test$band11_8 = raster::extract(band11,test_points)
df_test$bandtc_8 = extract(bandTC,test_points)
###### 2017-03-12
band1 = raster::raster(paste0(img.dir3,'/T34JEP_20170312T082001_B01.jp2'))
band2 = raster::raster(paste0(img.dir3,'/T34JEP_20170312T082001_B02.jp2'))
band3 = raster::raster(paste0(img.dir3,'/T34JEP_20170312T082001_B03.jp2'))
band4 = raster::raster(paste0(img.dir3,'/T34JEP_20170312T082001_B04.jp2'))
band5 = raster::raster(paste0(img.dir3,'/T34JEP_20170312T082001_B05.jp2'))
band6 = raster::raster(paste0(img.dir3,'/T34JEP_20170312T082001_B06.jp2'))
band7 = raster::raster(paste0(img.dir3,'/T34JEP_20170312T082001_B07.jp2'))
band8 = raster::raster(paste0(img.dir3,'/T34JEP_20170312T082001_B08.jp2'))
band8A = raster::raster(paste0(img.dir3,'/T34JEP_20170312T082001_B8A.jp2'))
band9 = raster::raster(paste0(img.dir3,'/T34JEP_20170312T082001_B09.jp2'))
band10 = raster::raster(paste0(img.dir3,'/T34JEP_20170312T082001_B10.jp2'))
band11 = raster::raster(paste0(img.dir3,'/T34JEP_20170312T082001_B11.jp2'))
band12 = raster::raster(paste0(img.dir3,'/T34JEP_20170312T082001_B12.jp2'))
bandTC3 = raster::stack(paste0(img.dir3,'/T34JEP_20170312T082001_TCI.jp2'))
#ndvi_11_5 = NDVI(band11,band5)
df_train$band1_12 = raster::extract(band1,train_points)
df_train$band2_12= raster::extract(band2,train_points)
df_train$band3_12 = raster::extract(band3,train_points)
df_train$band4_12 = raster::extract(band4,train_points)
df_train$band5_12 = raster::extract(band5,train_points)
df_train$band6_12 = raster::extract(band6,train_points)
df_train$band7_12 = raster::extract(band7,train_points)
df_train$band8_12 = raster::extract(band8,train_points)
df_train$band8A_12 = raster::extract(band8A,train_points)
df_train$band9_12 = raster::extract(band9,train_points)
df_train$band10_12 = raster::extract(band10,train_points)
df_train$band11_12 = raster::extract(band11,train_points)
df_train$band12_12 = raster::extract(band12,train_points)
df_train$bandtc_12 = extract(bandTC,train_points)
#df_train$ndvi_1_12 = extract(ndvi_11_5,train_points)
#
df_test$band1_12 = raster::extract(band1,test_points)
df_test$band2_12 = raster::extract(band2,test_points)
df_test$band3_12 = raster::extract(band3,test_points)
df_test$band4_12 = raster::extract(band4,test_points)
df_test$band5_12 = raster::extract(band5,test_points)
df_test$band6_12 = raster::extract(band6,test_points)
df_test$band7_12 = raster::extract(band7,test_points)
df_test$band8_12 = raster::extract(band8,test_points)
df_test$band9_12 = raster::extract(band9,test_points)
df_test$band10_12 = raster::extract(band10,test_points)
df_test$band11_12 = raster::extract(band11,test_points)
df_test$band12_12 = raster::extract(band12,test_points)
df_test$bandtc_12 = extract(bandTC,test_points)
#df_test$ndvi_1_12 = extract(ndvi_11_5,test_points)
###### 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'))
band8A = raster::raster(paste0(img.dir7,'/T34JEP_20170804T081559_B8A.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'))
band12 = raster::raster(paste0(img.dir7,'/T34JEP_20170804T081559_B12.jp2'))
bandTC = raster::raster(paste0(img.dir7,'/T34JEP_20170804T081559_TCI.jp2'))
#ndvi_11_5 = NDVI(band11,band5)
df_train$band1_7 = raster::extract(band1,train_points)
df_train$band2_7 = raster::extract(band2,train_points)
df_train$band3_7 = raster::extract(band3,train_points)
df_train$band4_7 = raster::extract(band4,train_points)
df_train$band5_7 = raster::extract(band5,train_points)
df_train$band6_7 = raster::extract(band6,train_points)
df_train$band7_7 = raster::extract(band7,train_points)
df_train$band8_7 = raster::extract(band8,train_points)
df_train$band8A_7 = raster::extract(band8A,train_points)
df_train$band9_7 = raster::extract(band9,train_points)
df_train$band10_7 = raster::extract(band10,train_points)
df_train$band11_7 = raster::extract(band11,train_points)
df_train$band12_7 = raster::extract(band12,train_points)
df_train$bandtc_7 = extract(bandTC,train_points)
#df_train$ndvi_1_7 = extract(ndvi_11_5,train_points)
df_test$band1_7 = raster::extract(band1,test_points)
df_test$band2_7 = raster::extract(band2,test_points)
df_test$band3_7 = raster::extract(band3,test_points)
df_test$band4_7 = raster::extract(band4,test_points)
df_test$band5_7 = raster::extract(band5,test_points)
df_test$band6_7 = raster::extract(band6,test_points)
df_test$band7_7 = raster::extract(band7,test_points)
df_test$band8_7 = raster::extract(band8,test_points)
df_test$band9_7 = raster::extract(band9,test_points)
df_test$band10_7 = raster::extract(band10,test_points)
df_test$band11_7 = raster::extract(band11,test_points)
df_test$bandtc_7 = extract(bandTC,test_points)
#df_test$ndvi_1_7 = extract(ndvi_11_5,test_points)
###### 2017-01-31
band1 = raster::raster(paste0(img.dir6,'/T34JEP_20170131T082151_B01.jp2'))
band2 = raster::raster(paste0(img.dir6,'/T34JEP_20170131T082151_B02.jp2'))
band3 = raster::raster(paste0(img.dir6,'/T34JEP_20170131T082151_B03.jp2'))
band4 = raster::raster(paste0(img.dir6,'/T34JEP_20170131T082151_B04.jp2'))
band5 = raster::raster(paste0(img.dir6,'/T34JEP_20170131T082151_B05.jp2'))
band6 = raster::raster(paste0(img.dir6,'/T34JEP_20170131T082151_B06.jp2'))
band7 = raster::raster(paste0(img.dir6,'/T34JEP_20170131T082151_B07.jp2'))
band8 = raster::raster(paste0(img.dir6,'/T34JEP_20170131T082151_B08.jp2'))
band8A = raster::raster(paste0(img.dir6,'/T34JEP_20170131T082151_B8A.jp2'))
band9 = raster::raster(paste0(img.dir6,'/T34JEP_20170131T082151_B09.jp2'))
band10 = raster::raster(paste0(img.dir6,'/T34JEP_20170131T082151_B10.jp2'))
band11 = raster::raster(paste0(img.dir6,'/T34JEP_20170131T082151_B11.jp2'))
bandTC = raster::raster(paste0(img.dir6,'/T34JEP_20170131T082151_TCI.jp2'))
df_train$band1_1 = raster::extract(band1,train_points)
df_train$band2_1 = raster::extract(band2,train_points)
#df_train$band3_1 = raster::extract(band3,train_points)#
df_train$band4_1 = raster::extract(band4,train_points)
df_train$band5_1 = raster::extract(band5,train_points)
df_train$band6_1 = raster::extract(band6,train_points)
#df_train$band7_1 = raster::extract(band7,train_points)#
df_train$band8_1 = raster::extract(band8,train_points)
df_train$band8A_1 = raster::extract(band8A,train_points)
df_train$band9_1 = raster::extract(band9,train_points)
df_train$band10_1 = raster::extract(band10,train_points)
df_train$band11_1 = raster::extract(band11,train_points)
df_train$bandtc_1 = raster::extract(bandTC,train_points)
df_test$band1_1 = raster::extract(band1,test_points)
df_test$band2_1 = raster::extract(band2,test_points)
df_test$band3_1 = raster::extract(band3,test_points)
df_test$band4_1 = raster::extract(band4,test_points)
df_test$band5_1 = raster::extract(band5,test_points)
df_test$band6_1 = raster::extract(band6,test_points)
df_test$band7_1 = raster::extract(band7,test_points)
df_test$band8_1 = raster::extract(band8,test_points)
df_test$band9_1 = raster::extract(band9,test_points)
df_test$band10_1 = raster::extract(band10,test_points)
df_test$band11_1 = raster::extract(band11,test_points)
df_test$bandtc_1 = raster::extract(bandTC,test_points)
#### 10-02-2017
band1 = raster::raster(paste0(img.dir2,'/T34JEP_20170210T082051_B01.jp2'))
band2 = raster::raster(paste0(img.dir2,'/T34JEP_20170210T082051_B02.jp2'))
band3 = raster::raster(paste0(img.dir2,'/T34JEP_20170210T082051_B03.jp2'))
band4 = raster::raster(paste0(img.dir2,'/T34JEP_20170210T082051_B04.jp2'))
band5 = raster::raster(paste0(img.dir2,'/T34JEP_20170210T082051_B05.jp2'))
band6 = raster::raster(paste0(img.dir2,'/T34JEP_20170210T082051_B06.jp2'))
band7 = raster::raster(paste0(img.dir2,'/T34JEP_20170210T082051_B07.jp2'))
band8 = raster::raster(paste0(img.dir2,'/T34JEP_20170210T082051_B08.jp2'))
band8A = raster::raster(paste0(img.dir2,'/T34JEP_20170210T082051_B8A.jp2'))
band9 = raster::raster(paste0(img.dir2,'/T34JEP_20170210T082051_B09.jp2'))
band10 = raster::raster(paste0(img.dir2,'/T34JEP_20170210T082051_B10.jp2'))
band11 = raster::raster(paste0(img.dir2,'/T34JEP_20170210T082051_B11.jp2'))
bandTC = raster::stack(paste0(img.dir2,'/T34JEP_20170210T082051_TCI.jp2'))
df_train$band1_6 = raster::extract(band1,train_points)
df_train$band2_6 = raster::extract(band2,train_points)
df_train$band3_6 = raster::extract(band3,train_points)
df_train$band4_6 = raster::extract(band4,train_points)
df_train$band5_6 = raster::extract(band5,train_points)
df_train$band6_6 = raster::extract(band6,train_points)
df_train$band7_6 = raster::extract(band7,train_points)
df_train$band8_6 = raster::extract(band8,train_points)
df_train$band8A_6 = raster::extract(band8A,train_points)
df_train$band9_6 = raster::extract(band9,train_points)
df_train$band10_6 = raster::extract(band10,train_points)
df_train$band11_6 = raster::extract(band11,train_points)
df_train$bandtc_6 = raster::extract(bandTC,train_points)
df_test$band1_6 = raster::extract(band1,test_points)
df_test$band2_6 = raster::extract(band2,test_points)
df_test$band3_6 = raster::extract(band3,test_points)
df_test$band4_6 = raster::extract(band4,test_points)
df_test$band5_6 = raster::extract(band5,test_points)
df_test$band6_6 = raster::extract(band6,test_points)
df_test$band7_6 = raster::extract(band7,test_points)
df_test$band8_6 = raster::extract(band8,test_points)
df_test$band9_6 = raster::extract(band8,test_points)
df_test$band10_6 = raster::extract(band8,test_points)
df_test$band11_6 = raster::extract(band8,test_points)
df_test$bandtc_6 = raster::extract(bandTC,test_points)
#### 15-07-2017
band1 = raster::raster(paste0(img.dir4,'/T34JEP_20170715T081609_B01.jp2'))
band2 = raster::raster(paste0(img.dir4,'/T34JEP_20170715T081609_B02.jp2'))
band3 = raster::raster(paste0(img.dir4,'/T34JEP_20170715T081609_B03.jp2'))
band4 = raster::raster(paste0(img.dir4,'/T34JEP_20170715T081609_B04.jp2'))
band5 = raster::raster(paste0(img.dir4,'/T34JEP_20170715T081609_B05.jp2'))
band6 = raster::raster(paste0(img.dir4,'/T34JEP_20170715T081609_B06.jp2'))
band7 = raster::raster(paste0(img.dir4,'/T34JEP_20170715T081609_B07.jp2'))
band8 = raster::raster(paste0(img.dir4,'/T34JEP_20170715T081609_B08.jp2'))
band8A = raster::raster(paste0(img.dir4,'/T34JEP_20170715T081609_B8A.jp2'))
band9 = raster::raster(paste0(img.dir4,'/T34JEP_20170715T081609_B09.jp2'))
band10 = raster::raster(paste0(img.dir4,'/T34JEP_20170715T081609_B10.jp2'))
band11 = raster::raster(paste0(img.dir4,'/T34JEP_20170715T081609_B11.jp2'))
bandTC = raster::raster(paste0(img.dir4,'/T34JEP_20170715T081609_TCI.jp2'))
df_train$band1_4 = raster::extract(band1,train_points)
df_train$band2_4 = raster::extract(band2,train_points)
df_train$band3_4 = raster::extract(band3,train_points)
df_train$band4_4 = raster::extract(band4,train_points)
df_train$band5_4 = raster::extract(band5,train_points)
df_train$band6_4 = raster::extract(band6,train_points)
df_train$band7_4 = raster::extract(band7,train_points)
df_train$band8_4 = raster::extract(band8,train_points)
df_train$band9_4 = raster::extract(band9,train_points)
df_train$band10_4 = raster::extract(band10,train_points)
df_train$band11_4 = raster::extract(band11,train_points)
df_train$bandtc_4 = raster::extract(bandTC,train_points)
df_test$band1_4 = raster::extract(band1,test_points)
df_test$band2_4 = raster::extract(band2,test_points)
df_test$band3_4 = raster::extract(band3,test_points)
df_test$band4_4 = raster::extract(band4,test_points)
df_test$band5_4 = raster::extract(band5,test_points)
df_test$band6_4 = raster::extract(band6,test_points)
df_test$band7_4 = raster::extract(band7,test_points)
df_test$band8_4 = raster::extract(band8,test_points)
df_test$band9_4 = raster::extract(band9,test_points)
df_test$band10_4 = raster::extract(band10,test_points)
df_test$band11_4 = raster::extract(band11,test_points)
df_test$bandtc_4 = raster::extract(bandTC,test_points)
#### 10-07-2017
band1 = raster::raster(paste0(img.dir42,'/T34JEP_20170710T082011_B01.jp2'))
band2 = raster::raster(paste0(img.dir42,'/T34JEP_20170710T082011_B02.jp2'))
band3 = raster::raster(paste0(img.dir42,'/T34JEP_20170710T082011_B03.jp2'))
band4 = raster::raster(paste0(img.dir42,'/T34JEP_20170710T082011_B04.jp2'))
band5 = raster::raster(paste0(img.dir42,'/T34JEP_20170710T082011_B05.jp2'))
band6 = raster::raster(paste0(img.dir42,'/T34JEP_20170710T082011_B06.jp2'))
band7 = raster::raster(paste0(img.dir42,'/T34JEP_20170710T082011_B07.jp2'))
band8 = raster::raster(paste0(img.dir42,'/T34JEP_20170710T082011_B08.jp2'))
band8A = raster::raster(paste0(img.dir42,'/T34JEP_20170710T082011_B8A.jp2'))
band9 = raster::raster(paste0(img.dir42,'/T34JEP_20170710T082011_B09.jp2'))
band10 = raster::raster(paste0(img.dir42,'/T34JEP_20170710T082011_B10.jp2'))
band11 = raster::raster(paste0(img.dir42,'/T34JEP_20170710T082011_B11.jp2'))
bandTC = raster::raster(paste0(img.dir42,'/T34JEP_20170710T082011_TCI.jp2'))
df_train$band1_42 = raster::extract(band1,train_points)
df_train$band2_42 = raster::extract(band2,train_points)
df_train$band3_42 = raster::extract(band3,train_points)
df_train$band4_42 = raster::extract(band4,train_points)
df_train$band5_42 = raster::extract(band5,train_points)
df_train$band6_42 = raster::extract(band6,train_points)
df_train$band7_42 = raster::extract(band7,train_points)
df_train$band8_42 = raster::extract(band8,train_points)
df_train$band8A_42 = raster::extract(band8A,train_points)
df_train$band9_42 = raster::extract(band9,train_points)
df_train$band10_42 = raster::extract(band10,train_points)
df_train$band11_42 = raster::extract(band11,train_points)
df_train$bandtc_42 = raster::extract(bandTC,train_points)
df_test$band1_42 = raster::extract(band1,test_points)
df_test$band2_42 = raster::extract(band2,test_points)
df_test$band3_42 = raster::extract(band3,test_points)
df_test$band4_42 = raster::extract(band4,test_points)
df_test$band5_42 = raster::extract(band5,test_points)
df_test$band6_42 = raster::extract(band6,test_points)
df_test$band7_42 = raster::extract(band7,test_points)
df_test$band8_42 = raster::extract(band8,test_points)
#df_test$band8A_42 = raster::extract(band8A,test_points)
df_test$band9_42 = raster::extract(band9,test_points)
df_test$band10_42 = raster::extract(band10,test_points)
df_test$band11_42 = raster::extract(band11,test_points)
df_test$bandtc_42 = raster::extract(bandTC,test_points)
#### 20-06-2017
band1 = raster::raster(paste0(img.dir66,'/T34JEP_20170620T082011_B01.jp2'))
band2 = raster::raster(paste0(img.dir66,'/T34JEP_20170620T082011_B02.jp2'))
band3 = raster::raster(paste0(img.dir66,'/T34JEP_20170620T082011_B03.jp2'))
band4 = raster::raster(paste0(img.dir66,'/T34JEP_20170620T082011_B04.jp2'))
band5 = raster::raster(paste0(img.dir66,'/T34JEP_20170620T082011_B05.jp2'))
band6 = raster::raster(paste0(img.dir66,'/T34JEP_20170620T082011_B06.jp2'))
band7 = raster::raster(paste0(img.dir66,'/T34JEP_20170620T082011_B07.jp2'))
band8 = raster::raster(paste0(img.dir66,'/T34JEP_20170620T082011_B08.jp2'))
band8A = raster::raster(paste0(img.dir66,'/T34JEP_20170620T082011_B8A.jp2'))
band9 = raster::raster(paste0(img.dir66,'/T34JEP_20170620T082011_B09.jp2'))
band10 = raster::raster(paste0(img.dir66,'/T34JEP_20170620T082011_B10.jp2'))
band11 = raster::raster(paste0(img.dir66,'/T34JEP_20170620T082011_B11.jp2'))
bandTC = raster::raster(paste0(img.dir66,'/T34JEP_20170620T082011_TCI.jp2'))
df_train$band1_66 = raster::extract(band1,train_points)
df_train$band2_66 = raster::extract(band2,train_points)
df_train$band3_66 = raster::extract(band3,train_points)
df_train$band4_66 = raster::extract(band4,train_points)
df_train$band5_66 = raster::extract(band5,train_points)
df_train$band6_66 = raster::extract(band6,train_points)
df_train$band7_66 = raster::extract(band7,train_points)
df_train$band8_66 = raster::extract(band8,train_points)
df_train$band8A_66 = raster::extract(band8A,train_points)
df_train$band9_66 = raster::extract(band9,train_points)
df_train$band10_66 = raster::extract(band10,train_points)
df_train$band11_66 = raster::extract(band11,train_points)
df_train$bandtc_66 = raster::extract(bandTC,train_points)
df_test$band1_66 = raster::extract(band1,test_points)
df_test$band2_66 = raster::extract(band2,test_points)
df_test$band3_66 = raster::extract(band3,test_points)
df_test$band4_66 = raster::extract(band4,test_points)
df_test$band5_66 = raster::extract(band5,test_points)
df_test$band6_66 = raster::extract(band6,test_points)
df_test$band7_66 = raster::extract(band7,test_points)
df_test$band8_66 = raster::extract(band8,test_points)
df_test$band9_66 = raster::extract(band9,test_points)
df_test$band10_66 = raster::extract(band10,test_points)
df_test$band11_66 = raster::extract(band11,test_points)
df_test$bandtc_66 = raster::extract(bandTC,test_points)
#### 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'))
band12 = raster::raster(paste0(img.dir1,'/T34JEP_20170322T081611_B12.jp2'))
bandTC = raster::raster(paste0(img.dir1,'/T34JEP_20170322T081611_TCI.jp2'))
df_train$band1_22 = raster::extract(band1,train_points)
df_train$band2_22 = raster::extract(band2,train_points)
df_train$band3_22 = raster::extract(band3,train_points)
df_train$band4_22 = raster::extract(band4,train_points)
df_train$band5_22 = raster::extract(band5,train_points)
df_train$band6_22 = raster::extract(band6,train_points)
df_train$band7_22 = raster::extract(band7,train_points)
df_train$band8_22 = raster::extract(band8,train_points)
df_train$band8A_22 = raster::extract(band8A,train_points)
df_train$band9_22 = raster::extract(band9,train_points)
df_train$band10_22 = raster::extract(band10,train_points)
df_train$band11_22 = raster::extract(band11,train_points)
df_train$band12_22 = raster::extract(band12,train_points)
df_train$bandtc_22 = raster::extract(bandTC,train_points)
df_test$band1_22 = raster::extract(band1,test_points)
df_test$band2_22 = raster::extract(band2,test_points)
df_test$band3_22 = raster::extract(band3,test_points)
df_test$band4_22 = raster::extract(band4,test_points)
df_test$band5_22 = raster::extract(band5,test_points)
df_test$band6_22 = raster::extract(band6,test_points)
df_test$band7_22 = raster::extract(band7,test_points)
df_test$band8_22 = raster::extract(band8,test_points)
df_test$band9_22 = raster::extract(band8,test_points)
df_test$band10_22 = raster::extract(band8,test_points)
df_test$band11_22 = raster::extract(band8,test_points)
df_test$bandtc_22 = raster::extract(band8,test_points)
## extract values from values
## extract values from values
features = raster::extract(a, trans)
###################
### XGB MODEL
###################
#####
duplicated.x = findCorrelation(cor(randomForest::na.roughfix(df_train)),
cutoff = 0.999, names = T, verbose = F)
length(duplicated.x)
df_train = df_train[, !names(df_train) %in% duplicated.x]
#### PCA Reduction
#run loop to create RasterStack with weighted mean of
#neighbour values for each pixel and each band
#img_2017_2dates_4band is my raster file with several dates and bands for the same location
# library(raster)
# m=matrix(c(1, 1, 1, 1, 1, 1, 1, 1, 1), ncol=3, nrow=3)
# p = raster::stack(band5,band6,band7,band11)
# names(p) = c("B5","B6","B7","B8")
# rs_focal=focal(p[[1]], fun=mean, w=a)
# for (i in 2:nlayers(p))
# {
# rl=focal(p[[i]], fun=mean, w=m)
# rs_focal=stack(rs_focal, rl)
# }
# names(rs_focal)=paste(names(p), "f")
# #end of loop
# k nearest neighbors
d = cbind(dr.dat$lon, dr.dat$lat)
k <- 8
nn <- knearneigh(d, k)
orstationc.neighbors.knn <- knn2nb(nn)
plot(orstationc.neighbors.knn, d, add=TRUE, col="green")
knn.feature = matrix(unlist(orstationc.neighbors.knn), ncol = 8)
####
mat = matrix(0,ncol = 4, nrow = 3568)
nbl = tri2nb(d)
a = c()
b = c()
for (i in 1:3568) {
a[i] = length(nbl[[i]])
b[i] = sum(nbl[[i]])
}
for (i in 1:3568) {
mat[i,] = nbl[[i]][1:4]
}
### Kmeans on Top features.
library(Rtsne)
#dat = df %>% dplyr::select(starts_with("stack"))
set.seed(1234)
tsne.res = Rtsne(umap_data2, check_duplicates = T, max_iter=1000,
perplexity = 50, #pca = TRUE,
theta = 0.5, dims = 2, verbose =T)
save(tsne.res, file = paste0(save.files.dir,"/tsne.RData"))
### Kmeans on top features
sumsq = NULL
for (i in 1:15) {
set.seed(1234)
sumsq[i] = sum(kmeans(umap_data2,centers = i, iter.max = 1000,
algorithm = "Lloyd")$withinss)
}
plot(1:15,sumsq,type= "b")
###
set.seed(1234)
kmns = kmeans(dr.dat,8,nstart = 17,iter.max = 1000,
algorithm = "Lloyd",trace = T)
#######
df_train$stack2_PC1 = ifelse(is.na(df_train$stack2_PC1),trainpca2$stack2_PC1,df_train$stack2_PC1)
df_train$stack2_PC2 = ifelse(is.na(df_train$stack2_PC2),trainpca2$stack2_PC2,df_train$stack2_PC2)
df_train$stack2_PC3 = ifelse(is.na(df_train$stack2_PC3),trainpca2$stack2_PC3,df_train$stack2_PC3)
df_train$stack2_PC4 = ifelse(is.na(df_train$stack2_PC4),trainpca2$stack2_PC4,df_train$stack2_PC4)
df_train$stack3_PC1 = ifelse(is.na(df_train$stack3_PC1),trainpca2$stack3_PC1,df_train$stack3_PC1)
df_train$stack3_PC2 = ifelse(is.na(df_train$stack3_PC2),trainpca2$stack3_PC2,df_train$stack3_PC2)
df_train$stack3_PC3 = ifelse(is.na(df_train$stack3_PC3),trainpca2$stack3_PC3,df_train$stack3_PC3)
df_train$stack3_PC4 = ifelse(is.na(df_train$stack3_PC4),trainpca2$stack3_PC4,df_train$stack3_PC4)
####
df_train$stack2_5_PC1 = ifelse(is.na(df_train$stack2_5_PC1),trainpca2$stack2_5_PC1,df_train$stack2_5_PC1)
df_train$stack2_5_PC2 = ifelse(is.na(df_train$stack2_5_PC2),trainpca2$stack2_5_PC2,df_train$stack2_5_PC2)
df_train$stack2_5_PC3 = ifelse(is.na(df_train$stack2_5_PC3),trainpca2$stack2_5_PC3,df_train$stack2_5_PC3)
df_train$stack3_5_PC1 = ifelse(is.na(df_train$stack3_5_PC1),trainpca2$stack3_5_PC1,df_train$stack3_5_PC1)
df_train$stack3_5_PC2 = ifelse(is.na(df_train$stack3_5_PC2),trainpca2$stack3_5_PC2,df_train$stack3_5_PC2)
df_train$stack3_5_PC3 = ifelse(is.na(df_train$stack3_5_PC3),trainpca2$stack3_5_PC3,df_train$stack3_5_PC3)
#####
df_train$stack2_12_PC1 = ifelse(is.na(df_train$stack2_12_PC1),trainpca2$stack2_12_PC1,df_train$stack2_12_PC1)
df_train$stack2_12_PC2 = ifelse(is.na(df_train$stack2_12_PC2),trainpca2$stack2_12_PC2,df_train$stack2_12_PC2)
df_train$stack2_12_PC3 = ifelse(is.na(df_train$stack2_12_PC3),trainpca2$stack2_12_PC3,df_train$stack2_12_PC3)
df_train$stack3_12_PC1 = ifelse(is.na(df_train$stack3_12_PC1),trainpca2$stack3_12_PC1,df_train$stack3_12_PC1)
df_train$stack3_12_PC2 = ifelse(is.na(df_train$stack3_12_PC2),trainpca2$stack3_12_PC2,df_train$stack3_12_PC2)
df_train$stack3_12_PC3 = ifelse(is.na(df_train$stack3_12_PC3),trainpca2$stack3_12_PC3,df_train$stack3_12_PC3)
####
df_train$stack2_1_PC1 = ifelse(is.na(df_train$stack2_1_PC1),trainpca6$stack2_1_PC1,df_train$stack2_1_PC1)
df_train$stack2_1_PC2 = ifelse(is.na(df_train$stack2_1_PC2),trainpca6$stack2_1_PC2,df_train$stack2_1_PC2)
df_train$stack3_1_PC1 = ifelse(is.na(df_train$stack3_1_PC1),trainpca6$stack3_1_PC1,df_train$stack3_1_PC1)
df_train$stack3_1_PC2 = ifelse(is.na(df_train$stack3_1_PC2),trainpca6$stack3_1_PC2,df_train$stack3_1_PC2)
########
df_train$stack2_2_PC1 = ifelse(is.na(df_train$stack2_2_PC1),trainpca6$stack2_2_PC1,df_train$stack2_2_PC1)
df_train$stack2_2_PC2 = ifelse(is.na(df_train$stack2_2_PC2),trainpca6$stack2_2_PC2,df_train$stack2_2_PC2)
df_train$stack3_2_PC1 = ifelse(is.na(df_train$stack3_2_PC1),trainpca6$stack3_2_PC1,df_train$stack3_2_PC1)
df_train$stack3_2_PC2 = ifelse(is.na(df_train$stack3_2_PC2),trainpca6$stack3_2_PC2,df_train$stack3_2_PC2)
#########
df_train$stack2_6_PC1 = ifelse(is.na(df_train$stack2_6_PC1),trainpca6$stack2_6_PC1,df_train$stack2_6_PC1)
df_train$stack2_6_PC2 = ifelse(is.na(df_train$stack2_6_PC2),trainpca6$stack2_6_PC2,df_train$stack2_6_PC2)
df_train$stack3_6_PC1 = ifelse(is.na(df_train$stack3_6_PC1),trainpca6$stack3_6_PC1,df_train$stack3_6_PC1)
df_train$stack3_6_PC2 = ifelse(is.na(df_train$stack3_6_PC2),trainpca6$stack3_6_PC2,df_train$stack3_6_PC2)
#######
df_test$stack2_PC1 = ifelse(is.na(df_test$stack2_PC1),testpca2$stack2_PC1,df_test$stack2_PC1)
df_test$stack2_PC2 = ifelse(is.na(df_test$stack2_PC2),testpca2$stack2_PC2,df_test$stack2_PC2)
df_test$stack2_PC3 = ifelse(is.na(df_test$stack2_PC3),testpca2$stack2_PC3,df_test$stack2_PC3)
df_test$stack2_PC4 = ifelse(is.na(df_test$stack2_PC4),testpca2$stack2_PC4,df_test$stack2_PC4)
df_test$stack3_PC1 = ifelse(is.na(df_test$stack3_PC1),testpca2$stack3_PC1,df_test$stack3_PC1)
df_test$stack3_PC2 = ifelse(is.na(df_test$stack3_PC2),testpca2$stack3_PC2,df_test$stack3_PC2)
df_test$stack3_PC3 = ifelse(is.na(df_test$stack3_PC3),testpca2$stack3_PC3,df_test$stack3_PC3)
df_test$stack3_PC4 = ifelse(is.na(df_test$stack3_PC4),testpca2$stack3_PC4,df_test$stack3_PC4)
####
df_test$stack2_5_PC1 = ifelse(is.na(df_test$stack2_5_PC1),testpca2$stack2_5_PC1,df_test$stack2_5_PC1)
df_test$stack2_5_PC2 = ifelse(is.na(df_test$stack2_5_PC2),testpca2$stack2_5_PC2,df_test$stack2_5_PC2)
df_test$stack2_5_PC3 = ifelse(is.na(df_test$stack2_5_PC3),testpca2$stack2_5_PC3,df_test$stack2_5_PC3)
df_test$stack3_5_PC1 = ifelse(is.na(df_test$stack3_5_PC1),testpca2$stack3_5_PC1,df_test$stack3_5_PC1)
df_test$stack3_5_PC2 = ifelse(is.na(df_test$stack3_5_PC2),testpca2$stack3_5_PC2,df_test$stack3_5_PC2)
df_test$stack3_5_PC3 = ifelse(is.na(df_test$stack3_5_PC3),testpca2$stack3_5_PC3,df_test$stack3_5_PC3)
#####
df_test$stack2_12_PC1 = ifelse(is.na(df_test$stack2_12_PC1),testpca2$stack2_12_PC1,df_test$stack2_12_PC1)
df_test$stack2_12_PC2 = ifelse(is.na(df_test$stack2_12_PC2),testpca2$stack2_12_PC2,df_test$stack2_12_PC2)
df_test$stack2_12_PC3 = ifelse(is.na(df_test$stack2_12_PC3),testpca2$stack2_12_PC3,df_test$stack2_12_PC3)
df_test$stack3_12_PC1 = ifelse(is.na(df_test$stack3_12_PC1),testpca2$stack3_12_PC1,df_test$stack3_12_PC1)
df_test$stack3_12_PC2 = ifelse(is.na(df_test$stack3_12_PC2),testpca2$stack3_12_PC2,df_test$stack3_12_PC2)
df_test$stack3_12_PC3 = ifelse(is.na(df_test$stack3_12_PC3),testpca2$stack3_12_PC3,df_test$stack3_12_PC3)
# Ht = function(DF1, ntimes = 1){
# w = function(k){
# s1 = dwt(k, filter = "haar")
# return(s1@V[[1]])
# }
# smt = DF1
# for (i in 1:ntimes) {
# smt = t(apply(smt,1,w))
# }
# return(data.frame(smt))
# }
#
# a = Ht(as.matrix(dat,9))
dat = df %>% dplyr::select(ends_with("_5"))
der = function(x,d=1){
df = t(diff(t(x), differences = d))
return(df)
}
a5 = der(a,1)
colnames(a5) = paste0("der5_",1:ncol(a5))
s = ssvd(as.matrix(j), k = 4,n =4,maxit = 1000)
colnames(s$u) = paste0("svd",1:ncol(s$u))
svd.features = data.frame(s$u)