Mercurial > repos > ecology > srs_spectral_indices

comparison Lib_preprocess_S2.r @ 0:a8dabbf47e15 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/Marie59/Sentinel_2A/srs_tools commit b32737c1642aa02cc672534e42c5cb4abe0cd3e7
author ecology
date Mon, 09 Jan 2023 13:39:08 +0000
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:a8dabbf47e15
1 # == == == == == == == == == == == == == == == == == == == == == == == == == == ==
2 # preprocS2
3 # Lib_preprocess_S2.R
4 # == == == == == == == == == == == == == == == == == == == == == == == == == == ==
5 # PROGRAMMERS:
6 # Jean-Baptiste FERET <jb.feret@teledetection.fr>
7 # Copyright 2021/08 Jean-Baptiste FERET
8 # == == == == == == == == == == == == == == == == == == == == == == == == == == ==
9 # This Library contains functions to preprocess Sentinel-2 images downloaded from
10 # different data hubs, such as THEIA, PEPS or SCIHUB
11 # == == == == == == == == == == == == == == == == == == == == == == == == == == ==
12
13 #" This function adjusts information from ENVI header
14 #"
15 #" @param dsn character. path where to store the stack
16 #" @param bands list. should include "bandname", and if possible "wavelength"
17 #" @param sensor character. Name of the sensor used to acquire the image
18 #" @param stretch boolean. Set TRUE to get 10% stretching at display for reflectance, mentioned in hdr only
19 #"
20 #" @return None
21 #" @importFrom utils read.table
22 #" @importFrom raster hdr raster
23 #" @export
24 adjust_envi_hdr <- function(dsn, bands, sensor = "Unknown", stretch = FALSE) {
25
26 # Edit hdr file to add metadata
27 hdr <- read_envi_header(get_hdr_name(dsn))
28 hdr$`band names` <- bands$bandname
29 if (length(bands$wavelength) == length(bands$bandname)) {
30 hdr$wavelength <- bands$wavelength
31 }else {
32 hdr$wavelength <- NULL
33 }
34 if (stretch == TRUE) {
35 hdr$`default stretch` <- "0.000000 1000.000000 linear"
36 }
37 hdr$`z plot range` <- NULL
38 hdr$`data ignore value` <- "-Inf"
39 hdr$`sensor type` <- sensor
40 write_envi_header(hdr = hdr, hdrpath = get_hdr_name(dsn))
41
42 # remove unnecessary files
43 file2remove <- paste(dsn, ".aux.xml", sep = "")
44 if (file.exists(file2remove)) file.remove(file2remove)
45 file2remove <- paste(dsn, ".prj", sep = "")
46 if (file.exists(file2remove)) file.remove(file2remove)
47 file2remove <- paste(dsn, ".stx", sep = "")
48 if (file.exists(file2remove)) file.remove(file2remove)
49 return(invisible())
50 }
51
52 #" This function saves reflectance files
53 #"
54 #" @param s2sat character. Sentinel-2 mission ("2A" or "2B")
55 #" @param tile_s2 character. S2 tile name (2 numbers + 3 letters)
56 #" @param dateacq_s2 double. date of acquisition
57 #"
58 #" @return s2mission character. name of the S2 mission (2A or 2B)
59 #" @importFrom sen2r safe_getMetadata check_scihub_connection s2_list
60 #" @export
61 check_s2mission <- function(s2sat, tile_s2, dateacq_s2) {
62
63 # is mission already defined by user?
64 if (!is.null(s2sat)) {
65 if (s2sat == "2A") {
66 s2mission <- "2A"
67 }else if (s2sat == "2B") {
68 s2mission <- "2B"
69 }else {
70 message("Could not identify if image from Sentinel-2A or -2B")
71 message("Defining central wavelength of spectral bands based on S2A")
72 s2mission <- "2A"
73 }
74 }else {
75 message("Could not identify if image from Sentinel-2A or -2B")
76 message("Defining central wavelength of spectral bands based on S2A")
77 s2mission <- "2A"
78 }
79 return(s2mission)
80 }
81
82 #" this function aims at computing directory size
83 #" @param path character. path for directory
84 #" @param recursive boolean . set T if recursive
85 #"
86 #" @return size_files numeric. size in bytes
87 #" - image stack
88 #" - path for individual band files corresponding to the stack
89 #" - path for vector (reprojected if needed)
90 #"
91 #" @importFrom raster raster
92 #" @importFrom tools file_path_sans_ext file_ext
93 #" @export
94 dir_size <- function(path, recursive = TRUE) {
95 stopifnot(is.character(path))
96 files <- list.files(path, full.names = TRUE, recursive = recursive)
97 vect_size <- sapply(files, function(x) file.size(x))
98 size_files <- sum(vect_size)
99 return(size_files)
100 }
101
102 #" This function reads S2 data from L2A directories downloaded from
103 #" various data hubs including THEIA, PEPS & SCIHUB (SAFE format & LaSRC)
104 #" @param path_dir_s2 character. path for S2 directory
105 #" @param path_vector character. path for vector file
106 #" @param s2source character. type of directory format (depends on atmospheric correction: SAFE produced from Sen2Cor)
107 #" @param resolution numeric. buffer applied to vector file (in meters)
108 #" @param interpolation character. method for resampling. default = "bilinear"
109 #" @param fre_sre character. SRE or FRE products from THEIA
110 #"
111 #" @return listout list.
112 #" - image stack
113 #" - path for individual band files corresponding to the stack
114 #" - path for vector (reprojected if needed)
115 #"
116 #" @importFrom raster raster
117 #" @importFrom tools file_path_sans_ext file_ext
118 #" @export
119 extract_from_s2_l2a <- function(path_dir_s2, path_vector = NULL, s2source = "SAFE",
120 resolution = 10, interpolation = "bilinear", fre_sre = "FRE") {
121 # Get list of paths corresponding to S2 bands and depending on S2 directory
122 s2_bands <- get_s2_bands(path_dir_s2 = path_dir_s2,
123 s2source = s2source,
124 resolution = resolution,
125 fre_sre = fre_sre)
126
127 if (length(s2_bands$s2bands_10m) > 0) {
128 rastmp <- raster::raster(s2_bands$s2bands_10m[[1]])
129 } else if (length(s2_bands$s2bands_20m) > 0) {
130 rastmp <- raster::raster(s2_bands$s2bands_20m[[1]])
131 }
132 # check if vector and raster share the same projection. if not, re-project vector
133 if (!is.null(path_vector)) {
134 raster_proj <- raster::projection(rastmp)
135 path_vector_reproj <- paste(tools::file_path_sans_ext(path_vector), "_reprojected.shp", sep = "")
136 path_vector <- reproject_shp(path_vector_init = path_vector,
137 newprojection = raster_proj,
138 path_vector_reproj = path_vector_reproj)
139 }
140 # Extract data corresponding to the vector footprint (if provided) & resample data if needed
141 if (length(s2_bands$s2bands_10m) > 0) {
142 stack_10m <- read_s2bands(s2_bands = s2_bands$s2bands_10m, path_vector = path_vector,
143 resampling = 1, interpolation = interpolation)
144 }
145 if (length(s2_bands$s2bands_20m) > 0) {
146 if (resolution == 10 && s2source != "LaSRC") {
147 resampling <- 2
148 }else {
149 resampling <- 1
150 }
151 stack_20m <- read_s2bands(s2_bands = s2_bands$s2bands_20m, path_vector = path_vector,
152 resampling = resampling, interpolation = interpolation)
153 }
154 # get full stack including 10m and 20m spatial resolution
155 if (length(s2_bands$s2bands_10m) > 0 && length(s2_bands$s2bands_20m) > 0) {
156 diffxstart <- attributes(stack_10m)$dimensions[[1]]$from - attributes(stack_20m)$dimensions[[1]]$from
157 diffxstop <- attributes(stack_10m)$dimensions[[1]]$to - attributes(stack_20m)$dimensions[[1]]$to
158 diffystart <- attributes(stack_10m)$dimensions[[2]]$from - attributes(stack_20m)$dimensions[[2]]$from
159 diffystop <- attributes(stack_10m)$dimensions[[2]]$to - attributes(stack_20m)$dimensions[[2]]$to
160 if (!diffxstop == 0) {
161 # size of 20m > size of 10m --> reduce 20m
162 # size of 10m > size of 20m --> reduce 10m
163 if (diffxstop > 0) {
164 stack_10m <- stack_10m[, 1:(dim(stack_10m)[1] - diffxstop), , ]
165 }else if (diffxstop < 0) {
166 stack_20m <- stack_20m[, 1:(dim(stack_20m)[1] + diffxstop), , ]
167 }
168 }
169 if (!diffystop == 0) {
170 if (diffystop > 0) {
171 stack_10m <- stack_10m[, , 1:(dim(stack_10m)[2] - diffystop), ]
172 }else if (diffystop < 0) {
173 stack_20m <- stack_20m[, , 1:(dim(stack_20m)[2] + diffystop), ]
174 }
175 }
176 if (!diffxstart == 0) {
177 if (diffxstart > 0) {
178 stack_20m <- stack_20m[, (1 + diffxstart):dim(stack_20m)[1], , ]
179 }else if (diffxstart < 0) {
180 stack_10m <- stack_10m[, (1 - diffxstart):dim(stack_10m)[1], , ]
181 }
182 }
183 if (!diffystart == 0) {
184 if (diffystart > 0) {
185 stack_20m <- stack_20m[, , (1 + diffystart):dim(stack_20m)[2], ]
186 }else if (diffystart < 0) {
187 stack_10m <- stack_10m[, , (1 - diffystart):dim(stack_10m)[2], ]
188 }
189 }
190 # reorder bands with increasing wavelength
191 s2bands <- c("B02", "B03", "B04", "B05", "B06", "B07", "B08", "B8A", "B11", "B12", "Cloud")
192 namebands <- c(names(s2_bands$s2bands_10m), names(s2_bands$s2bands_20m))
193 reorder_bands <- match(s2bands, namebands)
194 namebands <- namebands[reorder_bands]
195 listfiles <- c(stack_10m$attr, stack_20m$attr)[reorder_bands]
196
197 # adjust size to initial vector footprint without buffer
198 # --> buffer is needed in order to ensure that extraction following
199 # footprint of vector matches for images of different spatial resolution
200 # get bounding box corresponding to footprint of image or image subset
201 bb_xycoords <- get_bb(path_raster = listfiles[1],
202 path_vector = path_vector, buffer = 0)
203
204 # prepare reading data for extent defined by bounding box
205 nxoff <- bb_xycoords$UL$col
206 nyoff <- bb_xycoords$UL$row
207 nxsize <- bb_xycoords$UR$col - bb_xycoords$UL$col + 1
208 nysize <- bb_xycoords$LR$row - bb_xycoords$UR$row + 1
209 nbufxsize <- nxsize
210 nbufysize <- nysize
211 s2_stack <- stars::read_stars(listfiles, along = "band",
212 RasterIO = list(nXOff = nxoff, nYOff = nyoff,
213 nXSize = nxsize, nYSize = nysize,
214 nBufXSize = nbufxsize, nBufYSize = nbufysize,
215 resample = "nearest_neighbour"), proxy = TRUE)
216
217
218 names(s2_stack$attr) <- namebands
219 }else if (length(s2_bands$s2bands_10m) > 0) {
220 s2_stack <- stack_10m
221 namebands <- names(s2_bands$s2bands_10m)
222 names(s2_stack$attr) <- namebands
223 }else if (length(s2_bands$s2bands_20m) > 0) {
224 s2_stack <- stack_20m
225 namebands <- names(s2_bands$s2bands_20m)
226 names(s2_stack$attr) <- namebands
227 }
228
229 listout <- list("s2_stack" = s2_stack, "s2_bands" = s2_bands, "path_vector" = path_vector,
230 "namebands" = namebands)
231 return(listout)
232 }
233
234 #" This function gets coordinates of a bounding box defined by a vector (optional) and a raster
235 #"
236 #" @param path_raster character. path for raster file
237 #" @param path_vector character. path for vector file
238 #" @param buffer numeric. buffer applied to vector file (in meters)
239 #"
240 #" @return bb_xycoords list. Coordinates (in pixels) of the upper/lower right/left corners of bounding box
241 #" @export
242 get_bb <- function(path_raster, path_vector = NULL, buffer = 0) {
243
244 if (!is.null(path_vector)) {
245 # get bounding box with a 50m buffer in order to allow for interpolation
246 bb_xycoords <- get_bb_from_vector(path_raster = path_raster,
247 path_vector = path_vector,
248 buffer = buffer)
249 }else if (is.null(path_vector)) {
250 bb_xycoords <- get_bb_from_fullimage(path_raster)
251 }
252 return(bb_xycoords)
253 }
254
255 #" This function gets extreme coordinates of a bounding box corresponding to a full image
256 #"
257 #" @param path_raster character. path for raster file
258 #"
259 #" @return bb_xycoords list. Coordinates (in pixels) of the upper/lower right/left corners of bounding box
260 #" @importFrom raster raster
261 #" @export
262 get_bb_from_fullimage <- function(path_raster) {
263 # get raster coordinates corresponding to Full image
264 rasterobj <- raster::raster(path_raster)
265 bb_xycoords <- list()
266 bb_xycoords[["UL"]] <- data.frame("row" = 1, "col" = 1)
267 bb_xycoords[["UR"]] <- data.frame("row" = 1, "col" = dim(rasterobj)[2])
268 bb_xycoords[["LL"]] <- data.frame("row" = dim(rasterobj)[1], "col" = 1)
269 bb_xycoords[["LR"]] <- data.frame("row" = dim(rasterobj)[1], "col" = dim(rasterobj)[2])
270 return(bb_xycoords)
271 }
272
273 #" This gets bounding box corresponding to a vector from a raster (UL, UR, LL, LR corners)
274 #"
275 #" @param path_raster character. path for raster file
276 #" @param path_vector character. path for vector file
277 #" @param buffer numeric. buffer applied to vector file (in meters)
278 #"
279 #" @return bb_xycoords list. Coordinates (in pixels) of the upper/lower right/left corners of bounding box
280 #" @importFrom sf st_read st_bbox st_crop
281 #" @importFrom rgeos gbuffer bbox2SP
282 #" @importFrom sp SpatialPoints bbox
283 #" @importFrom raster projection extract extent raster
284 #" @importFrom methods as
285 #" @export
286 get_bb_from_vector <- function(path_raster, path_vector, buffer = 0) {
287
288 data_raster <- raster::raster(path_raster)
289 # extract BB coordinates from vector
290 bb_vector <- rgeos::gbuffer(spgeom = as(sf::st_read(dsn = path_vector, quiet = TRUE), "Spatial"),
291 width = buffer, byid = TRUE)
292 # extract BB coordinates from raster
293 bb_raster <- rgeos::bbox2SP(bbox = bbox(data_raster))
294 # compute intersection
295 intersect <- rgeos::gIntersection(bb_vector, bb_raster)
296 bbext <- raster::extent(intersect)
297 xmin <- bbext[1]
298 xmax <- bbext[2]
299 ymin <- bbext[3]
300 ymax <- bbext[4]
301 # get coordinates of bounding box corresponding to vector
302 corners <- list()
303 corners[["UR"]] <- sp::SpatialPoints(coords = cbind(xmax, ymax))
304 corners[["LR"]] <- sp::SpatialPoints(coords = cbind(xmax, ymin))
305 corners[["UL"]] <- sp::SpatialPoints(coords = cbind(xmin, ymax))
306 corners[["LL"]] <- sp::SpatialPoints(coords = cbind(xmin, ymin))
307 raster::projection(corners[["UL"]]) <- raster::projection(corners[["UR"]]) <-
308 raster::projection(corners[["LL"]]) <- raster::projection(corners[["LR"]]) <-
309 raster::projection(sf::st_read(dsn = path_vector, quiet = TRUE))
310 # get coordinates for corners of bounding box
311 bb_xycoords <- list()
312 for (corner in names(corners)) {
313 ex_df <- as.data.frame(raster::extract(data_raster, corners[[corner]], cellnumbers = TRUE))
314 colrow <- ind2sub(data_raster, ex_df$cell)
315 bb_xycoords[[corner]] <- data.frame("row" = colrow$row, "col" = colrow$col)
316 }
317 return(bb_xycoords)
318 }
319
320 #" get hdr name from image file name, assuming it is BIL format
321 #"
322 #" @param impath path of the image
323 #"
324 #" @return corresponding hdr
325 #" @import tools
326 #" @export
327 get_hdr_name <- function(impath) {
328 if (tools::file_ext(impath) == "") {
329 impathhdr <- paste(impath, ".hdr", sep = "")
330 }else if (tools::file_ext(impath) == "bil") {
331 impathhdr <- gsub(".bil", ".hdr", impath)
332 }else if (tools::file_ext(impath) == "zip") {
333 impathhdr <- gsub(".zip", ".hdr", impath)
334 }else {
335 impathhdr <- paste(tools::file_path_sans_ext(impath), ".hdr", sep = "")
336 }
337
338 if (!file.exists(impathhdr)) {
339 message("WARNING : COULD NOT FIND hdr FILE")
340 print(impathhdr)
341 message("Process may stop")
342 }
343 return(impathhdr)
344 }
345
346 #" This function returns path for the spectral bands to be used
347 #"
348 #" @param path_dir_s2 character. Path for the directory containing S2 data. either L2A .SAFE S2 file or THEIA directory
349 #" @param s2source character. defines if data comes from SciHub as SAFE directory, from THEIA or from LaSRC
350 #" @param resolution numeric. spatial resolution of the final image: 10m or 20m
351 #" @param fre_sre character. SRE or FRE products from THEIA
352 #"
353 #" @return listbands list. contains path for spectral bands corresponding to 10m and 20m resolution
354 #" @export
355 get_s2_bands <- function(path_dir_s2, s2source = "SAFE", resolution = 10, fre_sre = "FRE") {
356
357 if (s2source == "SAFE" || s2source == "Sen2Cor") {
358 listbands <- get_s2_bands_from_sen2cor(path_dir_s2 = path_dir_s2, resolution = resolution)
359 }else if (s2source == "THEIA") {
360 listbands <- get_s2_bands_from_theia(path_dir_s2 = path_dir_s2, resolution = resolution,
361 fre_sre = fre_sre)
362 }else if (s2source == "LaSRC") {
363 listbands <- get_s2_bands_from_lasrc(path_dir_s2 = path_dir_s2, resolution = resolution)
364 }else {
365 message("The data source (Atmospheric correction) for Sentinel-2 image is unknown")
366 message("Please provide S2 images from one of the following data sources:")
367 message("- LaSRC (atmospheric correction: LaSRC)")
368 message("- THEIA (atmospheric correction: MAJA)")
369 message("- SAFE (atmospheric correction: Sen2Cor)")
370 s2bands_10m <- s2bands_20m <- granule <- mtdfile <- metadata_msi <- metadata_lasrc <- NULL
371 listbands <- list("s2bands_10m" = s2bands_10m, "s2bands_20m" = s2bands_20m, "GRANULE" = granule,
372 "metadata" = mtdfile, "metadata_MSI" = metadata_msi,
373 "metadata_lasrc" = metadata_lasrc)
374 }
375 return(listbands)
376 }
377
378 #" This function returns path for the spectral bands in SAFE / sen2Cor directory
379 #"
380 #" @param path_dir_s2 character. Path for the SAFE directory containing S2 data
381 #" @param resolution numeric. spatial resolution of the final image: 10m or 20m
382 #"
383 #" @return listbands list. contains path for spectral bands corresponding to 10m and 20m resolution, as well name of as granule
384 #" @export
385 get_s2_bands_from_sen2cor <- function(path_dir_s2, resolution = 10) {
386 # build path for all bands
387 if (resolution == 10) {
388 b10m <- c("B02", "B03", "B04", "B08")
389 b20m <- c("B05", "B06", "B07", "B8A", "B11", "B12")
390 }else {
391 b10m <- c()
392 b20m <- c("B02", "B03", "B04", "B05", "B06", "B07", "B08", "B8A", "B11", "B12")
393 }
394 # get granule directory & path for corresponding metadata XML file
395 granule <- list.dirs(list.dirs(path_dir_s2, recursive = FALSE)[grep(pattern = "GRANULE",
396 x = list.dirs(path_dir_s2, recursive = FALSE))], recursive = FALSE)
397 mtdfile <- file.path(granule, "MTD_TL.xml")
398 if (file.exists(file.path(path_dir_s2, "MTD_MSIL2A.xml"))) {
399 mtd_msi_file <- file.path(path_dir_s2, "MTD_MSIL2A.xml")
400 } else {
401 mtd_msi_file <- NULL
402 }
403
404 # Define path for bands
405 s2bands_20m_dir <- file.path(granule, "IMG_DATA", "R20m")
406 s2bands_10m_dir <- file.path(granule, "IMG_DATA", "R10m")
407 s2bands_10m <- s2bands_20m <- list()
408 for (band in b20m) {
409 s2bands_20m[[band]] <- file.path(s2bands_20m_dir, list.files(s2bands_20m_dir, pattern = band))
410 }
411 for (band in b10m) {
412 s2bands_10m[[band]] <- file.path(s2bands_10m_dir, list.files(s2bands_10m_dir, pattern = band))
413 }
414 # get cloud mask
415 cloud <- "MSK_CLDPRB_20m"
416 cloud_20m_dir <- file.path(granule, "QI_DATA")
417 s2bands_20m[["Cloud"]] <- file.path(cloud_20m_dir, list.files(cloud_20m_dir, pattern = cloud))
418 listbands <- list("s2bands_10m" = s2bands_10m,
419 "s2bands_20m" = s2bands_20m,
420 "GRANULE" = granule,
421 "metadata" = mtdfile,
422 "metadata_MSI" = mtd_msi_file,
423 "metadata_lasrc" = NULL)
424 return(listbands)
425 }
426
427 #" This function returns path for the spectral bands in LaSRC directory
428 #"
429 #" @param path_dir_s2 character. Path for the SAFE directory containing S2 data
430 #" @param resolution numeric. spatial resolution of the final image: 10m or 20m
431 #"
432 #" @return listbands list. contains path for spectral bands corresponding to 10m and 20m resolution, as well name of as granule
433 #" @importFrom stringr str_subset
434 #" @export
435 get_s2_bands_from_lasrc <- function(path_dir_s2, resolution = 10) {
436
437 # get granule directory & path for corresponding metadata XML file
438 granule <- path_dir_s2
439 mtdfile <- file.path(granule, "MTD_TL.xml")
440 if (file.exists(file.path(path_dir_s2, "MTD_MSIL1C.xml"))) {
441 mtd_msi_file <- file.path(path_dir_s2, "MTD_MSIL1C.xml")
442 } else {
443 mtd_msi_file <- NULL
444 }
445
446 # build path for all bands
447 b10m <- c("band2", "band3", "band4", "band5", "band6", "band7", "band8", "band8a", "band11", "band12")
448 b10m_standard <- c("B02", "B03", "B04", "B05", "B06", "B07", "B08", "B8A", "B11", "B12")
449 # Define path for bands
450 s2bands_10m <- s2bands_20m <- list()
451 for (i in 1:seq_along(b10m)) {
452 s2bands_10m[[b10m_standard[i]]] <- file.path(path_dir_s2,
453 list.files(path_dir_s2,
454 pattern = paste(b10m[i], ".tif", sep = "")))
455 }
456
457 # get metadata file containing offset
458 mtd_lasrc <- str_subset(list.files(path_dir_s2, pattern = "S2"), ".xml$")
459 if (file.exists(file.path(path_dir_s2, mtd_lasrc))) {
460 metadata_lasrc <- file.path(path_dir_s2, mtd_lasrc)
461 } else {
462 metadata_lasrc <- NULL
463 }
464 # get cloud mask
465 cloud <- "CLM"
466 s2bands_10m[["Cloud"]] <- file.path(path_dir_s2, list.files(path_dir_s2, pattern = cloud))
467 listbands <- list("s2bands_10m" = s2bands_10m,
468 "s2bands_20m" = s2bands_20m,
469 "GRANULE" = granule,
470 "metadata" = mtdfile,
471 "metadata_MSI" = mtd_msi_file,
472 "metadata_lasrc" = metadata_lasrc)
473 return(listbands)
474 }
475
476 #" This function returns path for the spectral bands in THEIA directory
477 #"
478 #" @param path_dir_s2 character. Path for the SAFE directory containing S2 data
479 #" @param resolution numeric. spatial resolution of the final image: 10m or 20m
480 #" @param fre_sre character. SRE or FRE products from THEIA
481 #"
482 #" @return listbands list. contains path for spectral bands corresponding to 10m and 20m resolution, as well name of as granule
483 #" @export
484 get_s2_bands_from_theia <- function(path_dir_s2, resolution = 10, fre_sre = "FRE") {
485
486 # build path for all bands
487 if (resolution == 10) {
488 b10m <- c("B02", "B03", "B04", "B08")
489 b20m <- c("B05", "B06", "B07", "B8A", "B11", "B12")
490 } else {
491 b10m <- c()
492 b20m <- c("B02", "B03", "B04", "B05", "B06", "B07", "B08", "B8A", "B11", "B12")
493 }
494
495 # get path_tile_s2 directory & path for corresponding metadata XML file
496 path_tile_s2 <- list.dirs(path_dir_s2, recursive = FALSE)
497 files_tile_s2 <- list.files(path_tile_s2, recursive = FALSE)
498 mtdfile <- file.path(path_tile_s2, files_tile_s2[grep(pattern = "MTD_ALL.xml", x = files_tile_s2)])
499
500 # Define path for bands
501 s2bands_10m_dir <- s2bands_20m_dir <- path_tile_s2
502 s2bands_10m <- s2bands_20m <- list()
503 for (band in b20m) {
504 band_20m_pattern <- paste0(gsub("0", "", band), ".tif") # for THEAI band 2 is "B2" ("B02" for SAFE)
505 list_files_20m <- list.files(s2bands_20m_dir, pattern = band_20m_pattern)
506 s2bands_20m[[band]] <- file.path(s2bands_20m_dir, list_files_20m)[grep(pattern = fre_sre,
507 x = file.path(s2bands_20m_dir, list_files_20m))]
508 }
509 for (band in b10m) {
510 band_10m_pattern <- paste0(gsub("0", "", band), ".tif") # for THEAI band 2 is "B2" ("B02" for SAFE)
511 list_files_10m <- list.files(s2bands_10m_dir, pattern = band_10m_pattern)
512 s2bands_10m[[band]] <- file.path(s2bands_10m_dir, list_files_10m)[grep(pattern = fre_sre,
513 x = file.path(s2bands_10m_dir, list_files_10m))]
514 }
515
516 # get cloud mask 10m
517 cloud_10m <- "CLM_R1"
518 cloud_10m_dir <- file.path(path_tile_s2, "MASKS")
519 s2bands_10m[["Cloud"]] <- file.path(cloud_10m_dir, list.files(cloud_10m_dir, pattern = cloud_10m))
520
521 # get cloud mask 20m
522 cloud_20m <- "CLM_R2"
523 cloud_20m_dir <- file.path(path_tile_s2, "MASKS")
524 s2bands_20m[["Cloud"]] <- file.path(cloud_20m_dir, list.files(cloud_20m_dir, pattern = cloud_20m))
525
526 # return list bands
527 listbands <- list("s2bands_10m" = s2bands_10m,
528 "s2bands_20m" = s2bands_20m,
529 "path_tile_s2" = path_tile_s2,
530 "metadata" = mtdfile)
531 return(listbands)
532 }
533
534 #" This function check S2 data level:
535 #" - L2A: already atmospherically corrected
536 #" - L1C: requires atmospheric corrections with sen2cor
537 #"
538 #" @param prodname character. original name for the S2 image
539 #"
540 #" @return s2level character. S2 level: L1C or L2A
541 #" @export
542 get_s2_level <- function(prodname) {
543 prodname <- basename(prodname)
544 if (length(grep(pattern = "L1C_", x = prodname)) == 1) {
545 s2level <- "L1C"
546 } else if (length(grep(pattern = "L2A_", x = prodname)) == 1) {
547 s2level <- "L2A"
548 }
549 return(s2level)
550 }
551
552 #" This function gets tile from S2 image
553 #"
554 #" @param prodname character. original name for the S2 image
555 #"
556 #" @return tilename character
557 #" @importFrom tools file_path_sans_ext
558 #" @export
559 get_tile <- function(prodname) {
560 prodname <- basename(prodname)
561 tilename <- tools::file_path_sans_ext(gsub("_.*", "", gsub(".*_T", "", prodname)))
562 return(tilename)
563 }
564
565 #" This function gets acquisition date from S2 image
566 #"
567 #" @param prodname character. original name for the S2 image
568 #"
569 #" @return dateacq character
570 #" @export
571 get_date <- function(prodname) {
572 prodname <- basename(prodname)
573 dateacq <- as.Date(gsub("T.*", "", gsub(".*_20", "20", prodname)), format = "%Y%m%d")
574 return(dateacq)
575 }
576
577 #" download S2 L1C data from Copernicus hub or Google cloud
578 #"
579 #" @param list_safe safe object. produced with sen2r::s2_list
580 #" @param l1c_path character. path for storage of L1C image
581 #" @param path_vector path for a vector file
582 #" @param time_interval dates. time interval for S2 query
583 #" @param googlecloud boolean. set to TRUE if google cloud SDK is installed and
584 #" @param forcegoogle boolean. set to TRUE if only google requested
585 #" sen2r configured as an alternative hub for S2 download
586 #"
587 #" @return prodname character. S2 Product name
588 #" @importFrom sen2r safe_is_online s2_list s2_download s2_order check_gcloud
589 #" @export
590 get_s2_l1c_image <- function(list_safe, l1c_path, path_vector, time_interval,
591 googlecloud = FALSE, forcegoogle = FALSE) {
592 # Check if available from Copernicus hub first
593 copernicus_avail <- sen2r::safe_is_online(list_safe)
594 # if available: download
595 prodname <- attr(list_safe, which = "name")
596 if (file.exists(file.path(l1c_path, prodname))) {
597 message("L1C file already downloaded")
598 message(file.path(l1c_path, prodname))
599 } else {
600 if (copernicus_avail == TRUE && forcegoogle == FALSE) {
601 sen2r::s2_download(list_safe, outdir = l1c_path)
602 } else if (copernicus_avail == FALSE || forcegoogle == TRUE) {
603 # if not available and googlecloud==TRUE
604 if (googlecloud == TRUE) {
605 # check if google cloud SDK available from this computer
606 ggc <- sen2r::check_gcloud()
607 if (ggc == TRUE) {
608 message("downloading from Google cloud")
609 list_safe_ggc <- sen2r::s2_list(spatial_extent = sf::st_read(dsn = path_vector),
610 time_interval = time_interval,
611 server = "gcloud")
612 prodname <- attr(list_safe_ggc, which = "name")
613 if (file.exists(file.path(l1c_path, prodname))) {
614 message("L1C file already downloaded")
615 message(file.path(l1c_path, prodname))
616 } else {
617 sen2r::s2_download(list_safe_ggc, outdir = l1c_path)
618 # check if QI_DATA exists in DATASTRIP, and create it if not the case
619 datastrip_path <- file.path(l1c_path, prodname, "DATASTRIP")
620 dsdir <- list.dirs(datastrip_path, recursive = FALSE)
621 if (length(match(list.dirs(dsdir, recursive = FALSE, full.names = FALSE), "QI_DATA")) == 0) {
622 dir.create(file.path(dsdir, "QI_DATA"))
623 }
624 }
625 } else if (ggc == FALSE) {
626 message("googlecloud set to TRUE but missing")
627 message("Please install Google cloud SDK")
628 message("https://cloud.google.com/sdk/docs/install")
629 message("and/or set configuration of sen2r following instructions")
630 message("https://www.r-bloggers.com/2021/06/downloading-sentinel-2-archives-from-google-cloud-with-sen2r/")
631 }
632 }
633 }
634 if (copernicus_avail == FALSE && googlecloud == FALSE) {
635 message("S2 image in Long Term Archive (LTA)")
636 message("Ordering image from LTA")
637 message("This may take 1 day, please run your script later")
638 orders2 <- sen2r::s2_order(list_safe)
639 message("An alternative is possible with Google cloud SDK")
640 message("https://cloud.google.com/sdk/docs/install")
641 message("and/or set configuration of sen2r following instructions")
642 message("https://www.r-bloggers.com/2021/06/downloading-sentinel-2-archives-from-google-cloud-with-sen2r/")
643 }
644 }
645 return(prodname)
646 }
647
648 #" download S2 L2A data from Copernicus hub or convert L1C to L2A
649 #"
650 #" @param l2a_path character. path for storage of L2A image
651 #" @param spatial_extent path for a vector file
652 #" @param dateacq character. date of acquisition
653 #" @param deletel1c Boolean. set TRUE to delete L1C images
654 #" @param Sen2Cor Boolean. set TRUE to automatically perform atmospheric corrections using sen2Cor
655 #" @param googlecloud boolean. set to TRUE if google cloud SDK is installed and
656 #" sen2r configured as an alternative hub for S2 download
657 #"
658 #" @return pathl2a character. Path for L2A image
659 #" @importFrom sen2r s2_list s2_download
660 #" @importFrom R.utils getAbsolutePath
661
662 #" @export
663 get_s2_l2a_image <- function(l2a_path, spatial_extent, dateacq,
664 deletel1c = FALSE, sen2cor = TRUE,
665 googlecloud = FALSE) {
666
667 # Needs to be updated: define path for L1c data
668 l1c_path <- l2a_path
669 # define time interval
670 time_interval <- as.Date(c(dateacq, dateacq))
671 # get list S2 products corresponding to study area and date of interest using sen2r package
672 if (googlecloud == TRUE) {
673 server <- c("scihub", "gcloud")
674 } else if (googlecloud == FALSE) {
675 server <- "scihub"
676 }
677 list_safe <- sen2r::s2_list(spatial_extent = sf::st_read(dsn = spatial_extent),
678 time_interval = time_interval,
679 server = server, availability = "check")
680 # download products
681 sen2r::s2_download(list_safe, outdir = l2a_path)
682 # name all products
683 prodname <- attr(list_safe, which = "name")
684 prodfullpath <- file.path(l2a_path, prodname)
685 if (sen2cor == TRUE) {
686 for (imgname in prodname) {
687 s2level <- get_s2_level(imgname)
688 if (s2level == "L1C") {
689 datepattern <- gsub(pattern = "-", replacement = "", x = dateacq)
690 pathl2a <- s2_from_l1c_to_l2a(prodname = imgname, l1c_path = l2a_path, l2a_path = l2a_path,
691 datepattern = datepattern, tmp_path = NULL)
692 if (deletel1c == TRUE) {
693 unlink(x = R.utils::getAbsolutePath(file.path(l1c_path, prodname)),
694 recursive = TRUE, force = TRUE)
695 # delete from full path and add atmospherically corrected
696 whichimg <- grep(x = prodfullpath, pattern = imgname)
697 dateacq <- get_date(imgname)
698 tilename <- get_tile(imgname)
699 pathl2a <- list.files(path = l2a_path, pattern = tilename, full.names = TRUE)
700 pathl2a <- pathl2a[grep(x = pathl2a, pattern = dateacq)]
701 pathl2a <- pathl2a[grep(x = basename(pathl2a), pattern = "L2A")]
702 prodfullpath[whichimg] <- pathl2a
703 }
704 }
705 }
706 }
707
708 return(prodfullpath)
709 }
710
711 #" convert image coordinates from index to X-Y
712 #"
713 #" @param Raster image raster object
714 #" @param image_index coordinates corresponding to the raster
715 ind2sub <- function(data_raster, image_index) {
716 c <- ((image_index - 1) %% data_raster@ncols) + 1
717 r <- floor((image_index - 1) / data_raster@ncols) + 1
718 my_list <- list("col" = c, "row" = r)
719 return(my_list)
720 }
721
722 #" mosaicing a set of rasters
723 #"
724 #" @param list_rasters character. list of paths corresponding to rasters to mosaic
725 #" @param dst_mosaic character. path and name of mosaic produced
726 #" @param stretch boolean. Set TRUE to get 10% stretching at display for reflectance, mentioned in hdr only
727 #"
728 #" @return None
729 #" @importFrom gdalUtils mosaic_rasters
730 #" @importFrom raster hdr raster
731 #" @export
732 mosaic_rasters <- function(list_rasters, dst_mosaic, stretch = FALSE) {
733
734 # produce mosaic
735 gdalUtils::mosaic_rasters(gdalfile = list_rasters, dst_dataset = dst_mosaic,
736 separate = FALSE, of = "Ehdr", verbose = TRUE)
737
738 # convert hdr to ENVI format
739 raster::hdr(raster(dst_mosaic), format = "ENVI")
740 # add info to hdr based on initial rasters
741 hdr_init <- read_envi_header(get_hdr_name(list_rasters[1]))
742 hdr <- read_envi_header(get_hdr_name(dst_mosaic))
743 hdr$`band names` <- hdr_init$`band names`
744 hdr$wavelength <- hdr_init$wavelength
745 if (stretch == TRUE) {
746 hdr$`default stretch` <- "0.000000 1000.000000 linear"
747 }
748 hdr$`z plot range` <- NULL
749 hdr$`data ignore value` <- "-Inf"
750 hdr$`sensor type` <- hdr_init$`sensor type`
751 hdr$`coordinate system string` <- read.table(paste(file_path_sans_ext(dst_mosaic), ".prj", sep = ""))
752 write_envi_header(hdr = hdr, hdrpath = get_hdr_name(dst_mosaic))
753 return(invisible())
754 }
755
756 #" Reads ENVI hdr file
757 #"
758 #" @param hdrpath Path of the hdr file
759 #"
760 #" @return list of the content of the hdr file
761 #" @export
762 read_envi_header <- function(hdrpath) {
763 if (!grepl(".hdr$", hdrpath)) {
764 stop("File extension should be .hdr")
765 }
766 hdr <- readLines(hdrpath)
767 ## check ENVI at beginning of file
768 if (!grepl("ENVI", hdr[1])) {
769 stop("Not an ENVI header (ENVI keyword missing)")
770 } else {
771 hdr <- hdr [-1]
772 }
773 ## remove curly braces and put multi-line key-value-pairs into one line
774 hdr <- gsub("\\{([^}]*)\\}", "\\1", hdr)
775 l <- grep("\\{", hdr)
776 r <- grep("\\}", hdr)
777
778 if (length(l) != length(r)) {
779 stop("Error matching curly braces in header (differing numbers).")
780 }
781
782 if (any(r <= l)) {
783 stop("Mismatch of curly braces in header.")
784 }
785
786 hdr[l] <- sub("\\{", "", hdr[l])
787 hdr[r] <- sub("\\}", "", hdr[r])
788
789 for (i in rev(seq_along(l))) {
790 hdr <- c(
791 hdr [seq_len(l [i] - 1)],
792 paste(hdr [l [i]:r [i]], collapse = "\n"),
793 hdr [-seq_len(r [i])]
794 )
795 }
796
797 ## split key = value constructs into list with keys as names
798 hdr <- sapply(hdr, split_line, "=", USE.NAMES = FALSE)
799 names(hdr) <- tolower(names(hdr))
800
801 ## process numeric values
802 tmp <- names(hdr) %in% c(
803 "samples", "lines", "bands", "header offset", "data type",
804 "byte order", "default bands", "data ignore value",
805 "wavelength", "fwhm", "data gain values"
806 )
807 hdr [tmp] <- lapply(hdr [tmp], function(x) {
808 as.numeric(unlist(strsplit(x, ", ")))
809 })
810
811 return(hdr)
812 }
813
814 #" This function reads a list of files corresponding to S2 bands
815 #" S2 bands are expected to have uniform spatial resolution and footprint
816 #" @param s2_bands list. list of S2 bands obtained from get_s2_bands
817 #" @param path_vector path for a vector file
818 #" @param resampling numeric. resampling factor (default = 1, set to resampling = 2 to convert 20m into 10m resolution)
819 #" @param interpolation character. method for resampling. default = "bilinear"
820 #"
821 #" @return stack_s2 list. contains stack of S2 bands
822 #"
823 #" @importFrom stars read_stars
824 #" @importFrom sf st_bbox st_read st_crop
825 #" @export
826
827 read_s2bands <- function(s2_bands, path_vector = NULL,
828 resampling = 1, interpolation = "bilinear") {
829 # get bounding box corresponding to footprint of image or image subset
830 bb_xycoords <- get_bb(path_raster = s2_bands[[1]],
831 path_vector = path_vector, buffer = 50)
832
833 # prepare reading data for extent defined by bounding box
834 nxoff <- bb_xycoords$UL$col
835 nyoff <- bb_xycoords$UL$row
836 nxsize <- bb_xycoords$UR$col - bb_xycoords$UL$col + 1
837 nysize <- bb_xycoords$LR$row - bb_xycoords$UR$row + 1
838 nbufxsize <- resampling * nxsize
839 nbufysize <- resampling * nysize
840 if (resampling == 1) {
841 interpolation <- "nearest_neighbour"
842 }
843 # write interpolated individual bands in temp directory
844 tmpdir <- tempdir()
845 tmpfile <- list()
846 for (band in names(s2_bands)) {
847 stack_s2_tmp <- stars::read_stars(s2_bands[[band]], along = "band",
848 RasterIO = list(nXOff = nxoff, nYOff = nyoff,
849 nXSize = nxsize, nYSize = nysize,
850 nBufXSize = nbufxsize, nBufYSize = nbufysize,
851 resample = interpolation), proxy = FALSE)
852 if (!is.null(path_vector)) {
853 stack_s2_tmp <- sf::st_crop(x = stack_s2_tmp, y = st_bbox(st_read(dsn = path_vector, quiet = TRUE)))
854 }
855 tmpfile[[band]] <- file.path(tmpdir, tools::file_path_sans_ext(basename(s2_bands[[band]])))
856 if (band == "Cloud") {
857 stars::write_stars(obj = stack_s2_tmp, dsn = tmpfile[[band]],
858 driver = "ENVI", type = "Byte", overwrite = TRUE)
859 } else {
860 stars::write_stars(obj = stack_s2_tmp, dsn = tmpfile[[band]],
861 driver = "ENVI", type = "Int16", overwrite = TRUE)
862 }
863 gc()
864 }
865
866 stack_s2 <- stars::read_stars(tmpfile, along = "band", proxy = TRUE)
867 return(stack_s2)
868 }
869
870 #" This function reads a raster stack, and gets footprint as pixel coordinates or vector file as input
871 #" @param path_raster character. path for raster file
872 #" @param path_vector character. path for vector file
873 #" @param bbpix list. coordinates of pixels corresponding to a bounding box
874 #"
875 #" @return starsobj stars object corresponding to raster or raster subset
876 #"
877 #" @importFrom stars read_stars
878 #" @importFrom sf st_bbox st_read st_crop
879 #" @export
880 read_raster <- function(path_raster, path_vector = NULL, bbpix = NULL) {
881 # get bounding box corresponding to footprint of image or image subset
882 if (is.null(bbpix)) {
883 bb_xycoords <- get_bb(path_raster = path_raster,
884 path_vector = path_vector, buffer = 0)
885 } else {
886 bb_xycoords <- bbpix
887 }
888 # prepare reading data for extent defined by bounding box
889 nxoff <- bb_xycoords$UL$col
890 nyoff <- bb_xycoords$UL$row
891 nxsize <- bb_xycoords$UR$col - bb_xycoords$UL$col + 1
892 nysize <- bb_xycoords$LR$row - bb_xycoords$UR$row + 1
893 nbufxsize <- nxsize
894 nbufysize <- nysize
895 starsobj <- stars::read_stars(path_raster, along = "band",
896 RasterIO = list(nXOff = nxoff, nYOff = nyoff,
897 nXSize = nxsize, nYSize = nysize,
898 nBufXSize = nbufxsize, nBufYSize = nbufysize),
899 proxy = FALSE)
900 return(starsobj)
901 }
902
903 #" This function reprojects a shapefile and saves reprojected shapefile
904 #"
905 #" @param path_vector_init character. path for a shapefile to be reprojected
906 #" @param newprojection character. projection to be applied to path_vector_init
907 #" @param path_vector_reproj character. path for the reprojected shapefile
908 #"
909 #" @return path_vector character. path of the shapefile
910 #" - path_vector_init if the vector did not need reprojection
911 #" - path_vector_reproj if the vector needed reprojection
912 #"
913 #" @importFrom rgdal readOGR writeOGR
914 #" @importFrom sp spTransform
915 #" @importFrom raster projection
916 #" @export
917 reproject_shp <- function(path_vector_init, newprojection, path_vector_reproj) {
918
919 dir_vector_init <- dirname(path_vector_init)
920 # shapefile extension
921 fileext <- file_ext(basename(path_vector_init))
922 if (fileext == "shp") {
923 name_vector_init <- file_path_sans_ext(basename(path_vector_init))
924 vector_init_ogr <- rgdal::readOGR(dir_vector_init, name_vector_init, verbose = FALSE)
925 } else if (fileext == "kml") {
926 vector_init_ogr <- rgdal::readOGR(path_vector_init, verbose = FALSE)
927 }
928 vector_init_proj <- raster::projection(vector_init_ogr)
929
930 if (!vector_init_proj == newprojection) {
931 dir_vector_reproj <- dirname(path_vector_reproj)
932 name_vector_reproj <- file_path_sans_ext(basename(path_vector_reproj))
933 vector_reproj <- sp::spTransform(vector_init_ogr, newprojection)
934 rgdal::writeOGR(obj = vector_reproj, dsn = dir_vector_reproj, layer = name_vector_reproj,
935 driver = "ESRI Shapefile", overwrite_layer = TRUE)
936 path_vector <- path_vector_reproj
937 } else {
938 path_vector <- path_vector_init
939 }
940 return(path_vector)
941 }
942
943
944 #" perform atmospheric corrections to convert L1C to L2A data with Sen2cor
945 #"
946 #" @param prodname character. produced with sen2r::s2_list
947 #" @param l1c_path character. path of directory where L1C image is stored
948 #" @param l2a_path character. path of directory where L2A image is stored
949 #" @param datepattern character. pattern corresponding to date of acquisition to identify L2A directory
950 #" @param tmp_path character. path of temporary directory where L2A image is stored
951 #" sen2r configured as an alternative hub for S2 download
952 #"
953 #" @return pathl2a character. S2 Product name
954 #" @importFrom sen2r safe_is_online s2_list s2_download s2_order
955 #" @importFrom R.utils getAbsolutePath
956 #"
957 #" @export
958 s2_from_l1c_to_l2a <- function(prodname, l1c_path, l2a_path, datepattern, tmp_path = NULL) {
959
960 # define path for tmp directory
961 if (is.null(tmp_path)) {
962 tmp_path <- tempdir(check = TRUE)
963 }
964 tmp_prodlist <- prodname
965 # perform Sen2Cor atmospheric corrections
966 binpath <- sen2r::load_binpaths()
967 # 2- open a command prompt and directly run sen2cor with following command line
968 cmd <- paste(binpath$sen2cor,
969 "--output_dir", R.utils::getAbsolutePath(l2a_path),
970 R.utils::getAbsolutePath(file.path(l1c_path, prodname)), sep = " ")
971 system(cmd)
972 pathl2a <- list.files(path = l2a_path, pattern = datepattern, full.names = TRUE)
973
974 return(pathl2a)
975 }
976
977 #" This function saves cloud masks.
978 #" "cloudMask_Binary" is default binary mask with 0 where clouds are detected and 1 for clean pixels
979 #" "cloudMask_RAW" is the original cloud layer produced by atmospheric correction algorithm
980 #" --> may be useful to refine cloud mask
981 #"
982 #" @param s2_stars list. stars object containing raster data. Can be produced with function extract_from_s2_l2a
983 #" @param cloud_path character.
984 #" @param s2source character.
985 #" @param footprint character. path for vector file defining footprint of interest in the image
986 #" @param saveraw boolean. should the original cloud mask layer be saved?
987 #" @param maxchunk numeric. Size of individual chunks to be written (in Mb)
988 #"
989 #" @return list of cloudmasks (binary mask, and raw mask if required)
990 #" @importFrom sf st_read
991 #" @importFrom stars write_stars
992 #" @importFrom raster raster
993 #" @export
994 save_cloud_s2 <- function(s2_stars, cloud_path, s2source = "SAFE",
995 footprint = NULL, saveraw = FALSE, maxchunk = 256) {
996
997 whichcloud <- which(names(s2_stars$attr) == "Cloud")
998 # Save cloud mask
999 if (saveraw == TRUE) {
1000 cloudraw <- file.path(cloud_path, "CloudMask_RAW")
1001 obj <- stars::read_stars(s2_stars$attr[whichcloud], proxy = TRUE)
1002 sizeobj <- dim(obj)[1] * dim(obj)[2] / (1024**2)
1003 nbchunks <- ceiling(sizeobj / maxchunk)
1004 stars::write_stars(obj,
1005 dsn = cloudraw,
1006 driver = "ENVI",
1007 type = "Byte",
1008 chunk_size = c(dim(obj)[1], dim(obj)[2] / nbchunks),
1009 progress = TRUE)
1010 } else {
1011 cloudraw <- NULL
1012 }
1013 # Save cloud mask as in biodivMapR (0 = clouds, 1 = pixel ok)
1014 cloudmask <- stars::read_stars(s2_stars$attr[whichcloud], proxy = FALSE)
1015 if (s2source == "SAFE" || s2source == "THEIA") {
1016 cloudy <- which(cloudmask[[1]] > 0)
1017 sunny <- which(cloudmask[[1]] == 0)
1018 } else if (s2source == "LaSRC") {
1019 cloudy <- which(is.na(cloudmask[[1]]))
1020 sunny <- which(cloudmask[[1]] == 1)
1021 }
1022
1023 cloudmask[[1]][cloudy] <- 0
1024 cloudmask[[1]][sunny] <- 1
1025 cloudbin <- file.path(cloud_path, "CloudMask_Binary")
1026 stars::write_stars(cloudmask, dsn = cloudbin, driver = "ENVI", type = "Byte", overwrite = TRUE)
1027 cloudmasks <- list("BinaryMask" = cloudbin, "RawMask" = cloudraw)
1028 # delete temporary file
1029 file.remove(s2_stars$attr[whichcloud])
1030 if (file.exists(paste(s2_stars$attr[whichcloud], ".hdr", sep = ""))) file.remove(paste(s2_stars$attr[whichcloud], ".hdr", sep = ""))
1031 gc()
1032 return(cloudmasks)
1033 }
1034
1035 #" This function saves reflectance files
1036 #"
1037 #" @param s2_stars list. stars object containing raster data. Can be produced with function extract_from_s2_l2a
1038 #" @param refl_path character. path for reflectance file to be stored
1039 #" @param format character. file format for reflectance data
1040 #" @param datatype character. data type (integer, float, 16bits, 32bits...)
1041 #" @param s2sat character. Sentinel-2 mission ("2A" or "2B")
1042 #" @param tile_s2 character. S2 tile name (2 numbers + 3 letters)
1043 #" @param dateacq_s2 double. date of acquisition
1044 #" @param MTD character. path for metadata file
1045 #" @param MTD_MSI character. path for metadata MSI file
1046 #" @param mtd_lasrc character. path for metadata LaSRC file
1047 #" @param maxchunk numeric. Size of individual chunks to be written (in Mb)
1048 #"
1049 #" @return None
1050 #" @importFrom stars write_stars st_apply
1051 #" @importFrom XML xml
1052 #" @export
1053 save_reflectance_s2 <- function(s2_stars, refl_path, format = "ENVI", datatype = "Int16",
1054 s2sat = NULL, tile_s2 = NULL, dateacq_s2 = NULL,
1055 mtd = NULL, mtd_msi = NULL, mtd_lasrc = NULL,
1056 maxchunk = 256) {
1057 # identify if S2A or S2B, if possible
1058 s2mission <- check_s2mission(s2sat = s2sat, tile_s2 = tile_s2, dateacq_s2 = dateacq_s2)
1059
1060 # define central wavelength corresponding to each spectral band
1061 if (s2mission == "2A") {
1062 wl_s2 <- list("B02" = 496.6, "B03" = 560.0, "B04" = 664.5,
1063 "B05" = 703.9, "B06" = 740.2, "B07" = 782.5, "B08" = 835.1,
1064 "B8A" = 864.8, "B11" = 1613.7, "B12" = 2202.4)
1065 } else if (s2mission == "2B") {
1066 wl_s2 <- list("B02" = 492.1, "B03" = 559.0, "B04" = 665.0,
1067 "B05" = 703.8, "B06" = 739.1, "B07" = 779.7, "B08" = 833.0,
1068 "B8A" = 864.0, "B11" = 1610.4, "B12" = 2185.7)
1069 }
1070 if (s2mission == "2A") {
1071 sensor <- "Sentinel_2A"
1072 } else if (s2mission == "2B") {
1073 sensor <- "Sentinel_2B"
1074 }
1075
1076 # apply offset when necessary
1077 listbands_bis <- c("B2", "B3", "B4", "B5", "B6", "B7", "B8", "B8A", "B11", "B12")
1078 if (!is.null(mtd_msi) && is.null(mtd_lasrc)) {
1079 # read XML file containing info about geometry of acquisition
1080 s2xml <- XML::xmlToList(mtd_msi)
1081 xml_offset <- s2xml$General_Info$Product_Image_Characteristics$BOA_ADD_offset_VALUES_LIST
1082 bands <- lapply(s2xml$General_Info$Product_Image_Characteristics$Spectral_Information_List, "[[", 4)
1083 if (!is.null(xml_offset) && !is.null(bands)) {
1084 bandid <- lapply(bands, "[[", 1)
1085 bandname <- lapply(bands, "[[", 2)
1086 offset <- data.frame("bandname" = unlist(bandname),
1087 "bandid" = unlist(bandid),
1088 "offset" = unlist(lapply(xml_offset, "[[", 1)))
1089 selbands <- match(listbands_bis, offset$bandname)
1090 offset <- offset[selbands, ]
1091 boa_quantval <- as.numeric(s2xml$General_Info$Product_Image_Characteristics$QUANTIFICATION_VALUES_LIST$BOA_QUANTIFICATION_VALUE[1])
1092 } else {
1093 offset <- data.frame("bandname" = listbands_bis,
1094 "bandid" = c(1, 2, 3, 4, 5, 6, 7, 8, 11, 12),
1095 "offset" = 0)
1096 boa_quantval <- 10000
1097 }
1098 } else if (!is.null(mtd_lasrc)) {
1099 # read XML file containing info about geometry of acquisition
1100 s2xml <- XML::xmlToList(mtd_lasrc)
1101 attributes_lasrc <- s2xml$bands[[14]]$.attrs
1102 attributes_lasrc_df <- data.frame(attributes_lasrc)
1103 if (match("add_offset", rownames(attributes_lasrc_df)) > 0 && match("scale_factor", rownames(attributes_lasrc_df)) > 0) {
1104 xml_offset <- as.numeric(attributes_lasrc[["add_offset"]])
1105 boa_quantval <- 1 / as.numeric(attributes_lasrc[["scale_factor"]])
1106 offset <- data.frame("bandname" = listbands_bis,
1107 "bandid" = c(1, 2, 3, 4, 5, 6, 7, 8, 11, 12),
1108 "offset" = xml_offset)
1109 } else {
1110 offset <- data.frame("bandname" = listbands_bis,
1111 "bandid" = c(1, 2, 3, 4, 5, 6, 7, 8, 11, 12),
1112 "offset" = 0)
1113 boa_quantval <- 10000
1114 }
1115 } else {
1116 offset <- data.frame("bandname" = listbands_bis,
1117 "bandid" = c(1, 2, 3, 4, 5, 6, 7, 8, 11, 12),
1118 "offset" = 0)
1119 boa_quantval <- 10000
1120 }
1121
1122 # identify where spectral bands are in the stars object
1123 stars_spectral <- list()
1124 starsnames <- names(s2_stars$attr)
1125 stars_spectral$bandname <- starsnames[which(!starsnames == "Cloud")]
1126 stars_spectral$wavelength <- wl_s2[stars_spectral$bandname]
1127
1128 sortedwl <- names(wl_s2)
1129 reorder <- match(sortedwl, stars_spectral$bandname)
1130 elim <- which(is.na(reorder))
1131 if (length(elim) > 0) {
1132 reorder <- reorder[-elim]
1133 }
1134 pathr <- s2_stars$attr[reorder]
1135
1136 names(pathr) <- NULL
1137 s2_stars2 <- stars::read_stars(pathr, along = "band", proxy = TRUE)
1138 stars_spectral$bandname <- stars_spectral$bandname[reorder]
1139 stars_spectral$wavelength <- stars_spectral$wavelength[reorder]
1140
1141 uniqueoffset <- as.numeric(unique(offset$offset))
1142 if (length(uniqueoffset) > 1) {
1143 message("Warning: BOA offset differs between bands.")
1144 message("offset will not be applied to the final S2 reflectance raster")
1145 message("check metadata file to identify the offset applied on each band")
1146 print(mtd_msi)
1147 } else {
1148 message("applying offset to reflectance data")
1149 if (is.null(mtd_lasrc) || uniqueoffset == 0) {
1150 offsets2 <- function(x) (round(x + uniqueoffset) * (10000 / boa_quantval))
1151 s2_stars2 <- stars::st_apply(X = s2_stars2, MARGIN = "band", FUN = offsets2)
1152 } else {
1153 offsets2 <- function(x) (round(10000 * ((x + uniqueoffset * boa_quantval) / boa_quantval)))
1154 s2_stars2 <- stars::st_apply(X = s2_stars2, MARGIN = "band", FUN = offsets2)
1155 }
1156 }
1157 write_stack_s2(stars_s2 = s2_stars2, stars_spectral = stars_spectral, refl_path = refl_path,
1158 format = format, datatype = datatype, sensor = sensor, maxchunk = maxchunk)
1159 # save metadata file as well if available
1160 if (!is.null(mtd)) {
1161 if (file.exists(mtd)) {
1162 file.copy(from = mtd, to = file.path(dirname(refl_path), basename(mtd)), overwrite = TRUE)
1163 }
1164 }
1165 # save metadata file as well if available
1166 if (!is.null(mtd_msi)) {
1167 if (file.exists(mtd_msi)) {
1168 file.copy(from = mtd_msi, to = file.path(dirname(refl_path), basename(mtd_msi)), overwrite = TRUE)
1169 }
1170 }
1171 # save LaSRC metadata file as well if available
1172 if (!is.null(mtd_lasrc)) {
1173 if (file.exists(mtd_lasrc)) {
1174 file.copy(from = mtd_lasrc, to = file.path(dirname(refl_path), basename(mtd_lasrc)), overwrite = TRUE)
1175 }
1176 }
1177 # delete temporary file
1178 for (pathtemp in pathr) {
1179 file.remove(pathtemp)
1180 if (file.exists(paste(pathtemp, ".hdr", sep = ""))) file.remove(paste(pathtemp, ".hdr", sep = ""))
1181 }
1182 gc()
1183 return(invisible())
1184 }
1185
1186 #" ENVI functions
1187 #"
1188 #" based on https://github.com/cran/hyperSpec/blob/master/R/read.ENVI.R
1189 #" added wavelength, fwhm, ... to header reading
1190 #" Title
1191 #"
1192 #" @param x character.
1193 #" @param separator character
1194 #" @param trim_blank boolean.
1195 #"
1196 #" @return list.
1197 #" @export
1198 split_line <- function(x, separator, trim_blank = TRUE) {
1199 tmp <- regexpr(separator, x)
1200 key <- substr(x, 1, tmp - 1)
1201 value <- substr(x, tmp + 1, nchar(x))
1202 if (trim_blank) {
1203 blank_pattern <- "^[[:blank:]]*([^[:blank:]]+.*[^[:blank:]]+)[[:blank:]]*$"
1204 key <- sub(blank_pattern, "\\1", key)
1205 value <- sub(blank_pattern, "\\1", value)
1206 }
1207 value <- as.list(value)
1208 names(value) <- key
1209 return(value)
1210 }
1211
1212 #" save raster footprint as vector file
1213 #"
1214 #" @param path_raster character. path for a raster file
1215 #" @param path_vector character. path for a vector file
1216 #" @param driver character. driver for vector
1217 #"
1218 #" @return None
1219 #" @importFrom raster raster extent projection
1220 #" @importFrom sf st_as_sf st_write
1221 #" @export
1222 vectorize_raster_extent <- function(path_raster, path_vector, driver = "ESRI Shapefile") {
1223 rast <- raster(path_raster)
1224 e <- extent(rast)
1225 # coerce to a SpatialPolygons object
1226 p <- as(e, "SpatialPolygons")
1227 projection(p) <- projection(rast)
1228 p <- sf::st_as_sf(p)
1229 sf::st_write(obj = p, path_vector, driver = driver) # create to a shapefile
1230 return(invisible())
1231 }
1232
1233 #" writes ENVI hdr file
1234 #"
1235 #" @param hdr content to be written
1236 #" @param hdrpath Path of the hdr file
1237 #"
1238 #" @return None
1239 #" @importFrom stringr str_count
1240 #" @export
1241 write_envi_header <- function(hdr, hdrpath) {
1242 h <- lapply(hdr, function(x) {
1243 if (length(x) > 1 || (is.character(x) && stringr::str_count(x, "\\w+") > 1)) {
1244 x <- paste0("{", paste(x, collapse = ", "), "}")
1245 }
1246 # convert last numerics
1247 x <- as.character(x)
1248 })
1249 writeLines(c("ENVI", paste(names(hdr), h, sep = " = ")), con = hdrpath)
1250 return(invisible())
1251 }
1252
1253 #" This function writes a raster Stack object into a ENVI raster file
1254 #"
1255 #" @param stackobj list. raster stack
1256 #" @param stackpath character. path where to store the stack
1257 #" @param bands list. should include "bandname", and if possible "wavelength"
1258 #" @param datatype character. should be INT2S or FLT4S for example
1259 #" @param sensor character. Name of the sensor used to acquire the image
1260 #" @param stretch boolean. Set TRUE to get 10% stretching at display for reflectance, mentioned in hdr only
1261 #"
1262 #" @return None
1263 #" @importFrom utils read.table
1264 #" @export
1265 write_rasterstack_envi <- function(stackobj, stackpath, bands, datatype = "INT2S",
1266 sensor = "Unknown", stretch = FALSE) {
1267
1268 r <- raster::writeRaster(x = stackobj, filename = stackpath, format = "Ehdr", overwrite = TRUE, datatype = datatype)
1269 raster::hdr(r, format = "ENVI")
1270 # Edit hdr file to add metadata
1271 hdr <- read_envi_header(get_hdr_name(stackpath))
1272 hdr$`band names` <- bands$bandname
1273 if (length(bands$wavelength) == length(bands$bandname)) {
1274 hdr$wavelength <- bands$wavelength
1275 } else {
1276 hdr$wavelength <- NULL
1277 }
1278 if (stretch == TRUE) {
1279 hdr$`default stretch` <- "0.000000 1000.000000 linear"
1280 }
1281 hdr$`z plot range` <- NULL
1282 hdr$`data ignore value` <- "-Inf"
1283 hdr$`coordinate system string` <- read.table(paste(stackpath, ".prj", sep = ""))
1284 proj <- strsplit(x = strsplit(x = projection(stackobj), split = " ")[[1]][1], split = "=")[[1]][2]
1285 zone <- strsplit(x = strsplit(x = projection(stackobj), split = " ")[[1]][2], split = "=")[[1]][2]
1286 datum <- strsplit(x = strsplit(x = projection(stackobj), split = " ")[[1]][3], split = "=")[[1]][2]
1287 oldproj <- hdr$`map info`
1288 newproj <- gsub(pattern = "projection", replacement = proj, x = oldproj)
1289 newproj <- paste(newproj, zone, datum, sep = ", ")
1290 hdr$`map info` <- newproj
1291 hdr$`sensor type` <- sensor
1292 write_envi_header(hdr = hdr, hdrpath = get_hdr_name(stackpath))
1293
1294 # remove unnecessary files
1295 file2remove <- paste(stackpath, ".aux.xml", sep = "")
1296 file.remove(file2remove)
1297 file2remove <- paste(stackpath, ".prj", sep = "")
1298 file.remove(file2remove)
1299 file2remove <- paste(stackpath, ".stx", sep = "")
1300 file.remove(file2remove)
1301 return(invisible())
1302 }
1303
1304
1305 #" This function writes a stars object into a raster file
1306 #"
1307 #" @param stars_s2 list. stars object containing raster data. Can be produced with function Crop_n_resample_S2
1308 #" @param stars_spectral list. band name to be saved in the stack and spectral bands corresponding to the image
1309 #" @param refl_path character. path where to store the image
1310 #" @param format character. default = ENVI BSQ. otherwise use gdal drivers
1311 #" @param datatype character. should be Int16 or Float64 for example
1312 #" @param sensor character. Name of the sensor used to acquire the image
1313 #" @param maxchunk numeric. Size of individual chunks to be written (in Mb)
1314 #"
1315 #" @return None
1316 #" @export
1317 write_stack_s2 <- function(stars_s2, stars_spectral, refl_path, format = "ENVI",
1318 datatype = "Int16", sensor = "Unknown", maxchunk = 256) {
1319
1320 # write raster file from proxy using chunks
1321 sizeobj <- 2 * dim(stars_s2)[1] * dim(stars_s2)[2] * dim(stars_s2)[3] / (1024**2)
1322 nbchunks <- ceiling(sizeobj / maxchunk)
1323 stars::write_stars(obj = stars_s2,
1324 dsn = refl_path,
1325 driver = format,
1326 type = datatype,
1327 chunk_size = c(dim(stars_s2)[1], ceiling(dim(stars_s2)[2] / nbchunks)),
1328 progress = TRUE)
1329
1330 if (format == "ENVI") {
1331 adjust_envi_hdr(dsn = refl_path, bands = stars_spectral,
1332 sensor = sensor, stretch = TRUE)
1333 }
1334 return(invisible())
1335 }