-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathMosaicFreeCloud.R
71 lines (57 loc) · 2.14 KB
/
MosaicFreeCloud.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
#' MosaicFreeCloud
#'
#' @author Yonatan Tarazona, Vasco M. Mantas, A.J.S.C. Pereira
#'
#' @description MasaicFreeCloud is an algorithm for the generation of cloud-free images through an indeterminate number of optical images of a single band.
#'
#' @param list.r is a list where all the images are.
#'
#' @param MosFCL It's the raster after taking good pixels from the remaining images
#'
#' @details Check Readme.md
lista.r<-list.files("~/", ".tif", all.files=T, recursive=T, full.names=F)
## Extracting the dimension and resolution of a raster
raster <- raster(lista.r[1])
ext<- extent(raster@extent[1],raster@extent[2],raster@extent[3],raster@extent[4]) #extension
res<- raster(ext, nrows=raster@nrows, ncols=raster@ncols) # spatial resolution
## Equal dimensions for all images
vfs <- c() # empty vector
for(k in 1:length(lista.r)){
band <- raster(lista.r[k])
band.r<-resample(band, res, method="ngb")
band.r[is.na(band.r) | band.r==-1] <- 0
vfs <- c(vfs,band.r) # all images with same resolution and extension
}
## Getting the order of the images according to the atmospheric noise
lista <- vfs # copying
suma <- c() # empty vector
for(j in 1:length(lista.r)){
val <- getValues(lista[[j]])
sum <- sum(na.omit(val))
suma <- c(suma, sum)
orden <- order(suma, decreasing = TRUE)
bmax <- lista[[orden[1]]] # The most free image of atmospheric noise
}
## Calculating the distances (time) of the selected image with respect to the others
time <- as.numeric(substring(lista.r,10,17)) # from 10 to 17 is the time of acquisition of the image
t.max <- time[orden[1]]
t.dif <- order(abs(t.max - time),decreasing = FALSE)
## Generating the cloud-free mosaic
MosFCL <- bmax # copying
for(l in t.dif){
MosFCL <- MosFCL + lista[[l]]
MosFCL[bmax > 0] <- 0
MosFCL[bmax < 0] <- 0
MosFCL<- MosFCL+bmax
bmax <- MosFCL
}
## Assigning NA to pixels equal to 0
MosFCL[MosFCL==0]<-NA
## Assigning projection system to the final mosaic
projection(MosFCL)<-CRS(as.character(NA))
proj <- CRS(projection(raster))
projection(MosFCL)<-proj
## Saving the final mosaic in .tif
writeRaster(MosFCL,'Name.tif',drivername="GTiff", datatype = "FLT4S")
#' @ That's all!
#'