Mercurial > repos > ecology > srs_process_data

comparison pca_raster.r @ 0:82ee94fbefd7 draft default tip

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:53:22 +0000
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:82ee94fbefd7
1 #Rscript
2
3 ###########################################
4 ## Getting PCA raster ##
5 ###########################################
6
7 #####Packages : stars
8 # utils
9 # biodivmapr
10 # raster
11 # sf
12 # mapview
13 # leafpop
14 # RColorBrewer
15 # labdsv
16 # rgdal
17 # ggplot2
18 # gridExtra
19 ## remotes::install_github("jbferet/biodivMapR")
20 #####Load arguments
21
22 args <- commandArgs(trailingOnly = TRUE)
23
24 #####Import the S2 data
25
26 if (length(args) < 1) {
27 stop("This tool needs at least 1 argument")
28 }else {
29 data_raster <- args[1]
30 rasterheader <- args[2]
31 data <- args[3]
32 typepca <- as.character(args[4])
33 source(args[5])
34 }
35
36 ################################################################################
37 ## DEFINE PARAMETERS FOR DATASET TO BE PROCESSED ##
38 ################################################################################
39 # expected to be in ENVI HDR
40
41 if (data_raster == "") {
42 #Create a directory where to unzip your folder of data
43 dir.create("data_dir")
44 unzip(data, exdir = "data_dir")
45 # Path to raster
46 data_raster <- list.files("data_dir/results/Reflectance", pattern = "_Refl")
47 input_image_file <- file.path("data_dir/results/Reflectance", data_raster[1])
48 input_header_file <- file.path("data_dir/results/Reflectance", data_raster[2])
49
50 } else {
51 input_image_file <- file.path(getwd(), data_raster, fsep = "/")
52 input_header_file <- file.path(getwd(), rasterheader, fsep = "/")
53 }
54
55 ################################################################################
56 ## PROCESS IMAGE ##
57 ################################################################################
58 # 1- Filter data in order to discard non vegetated / shaded / cloudy pixels
59 print("PERFORM PCA ON RASTER")
60 pca_output <- biodivMapR::perform_PCA(Input_Image_File = input_image_file, Input_Mask_File = input_mask_file,
61 Output_Dir = output_dir, TypePCA = typepca, FilterPCA = filterpca, nbCPU = nbcpu, MaxRAM = maxram)
62
63
64 pca_path <- file.path(output_dir, basename(data_raster), typepca, "PCA", "OutputPCA_8_PCs")
65 pca_raster <- raster::raster(pca_path)
66 get_pca <- convert_raster(pca_raster)
67
68 colnames(get_pca) <- c("PCA", "longitude", "latitude")
69 plot_indices(get_pca, titre = "PCA")
70
71 write.table(get_pca, file = "PCA.tabular", sep = "\t", dec = ".", na = " ", row.names = FALSE, col.names = TRUE, quote = FALSE)
72 #### Get the raster layer files
73 pca_files <- file.path("RESULTS", basename(data_raster), typepca, "PCA")
74 to_dir_short <- output_dir
75 file.copy(pca_files, to_dir_short) #copy files from long to short paths