Mercurial > repos > ecology > srs_pca

comparison functions.r @ 0:79bcb8c6dff0 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:55:34 +0000
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:79bcb8c6dff0
1 #Rscript
2
3 #########################
4 ## Functions ##
5 #########################
6
7 #####Packages : raster
8 # sp
9 # ggplot2
10
11 ####Set paramaters for all tools using BiodivMapR
12
13 # path for the Mask raster corresponding to image to process
14 # expected to be in ENVI HDR format, 1 band, integer 8bits
15 # expected values in the raster: 0 = masked, 1 = selected
16 # set to FALSE if no mask available
17 input_mask_file <- FALSE
18
19 # relative or absolute path for the Directory where results will be stored
20 # For each image processed, a subdirectory will be created after its name
21 output_dir <- "RESULTS"
22
23 # SPATIAL RESOLUTION
24 # resolution of spatial units for alpha and beta diversity maps (in pixels), relative to original image
25 # if Res.Map = 10 for images with 10 m spatial resolution, then spatial units will be 10 pixels x 10m = 100m x 100m surfaces
26 # rule of thumb: spatial units between 0.25 and 4 ha usually match with ground data
27 # too small window_size results in low number of pixels per spatial unit, hence limited range of variation of diversity in the image
28 window_size <- 10
29
30 # PCA FILTERING: Set to TRUE if you want second filtering based on PCA outliers to be processed. Slower
31 filterpca <- TRUE
32
33 ################################################################################
34 ## DEFINE PARAMETERS FOR METHOD ##
35 ################################################################################
36 nbcpu <- 4
37 maxram <- 0.5
38 nbclusters <- 50
39
40 ################################################################################
41 ## PROCESS IMAGE ##
42 ################################################################################
43 # 1- Filter data in order to discard non vegetated / shaded / cloudy pixels
44 ndvi_thresh <- 0.5
45 blue_thresh <- 500
46 nir_thresh <- 1500
47 continuum_removal <- TRUE
48
49
50
51 #### Convert raster to dataframe
52
53 # Convert raster to SpatialPointsDataFrame
54 convert_raster <- function(data_raster) {
55 r_pts <- raster::rasterToPoints(data_raster, spatial = TRUE)
56
57 # reproject sp object
58 geo_prj <- "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"
59 r_pts <- sp::spTransform(r_pts, sp::CRS(geo_prj))
60
61
62 # Assign coordinates to @data slot, display first 6 rows of data.frame
63 r_pts@data <- data.frame(r_pts@data, longitude = sp::coordinates(r_pts)[, 1],
64 latitude = sp::coordinates(r_pts)[, 2])
65
66 return(r_pts@data)
67 }
68
69
70 #### Potting indices
71
72 plot_indices <- function(data, titre) {
73 graph_indices <- ggplot2::ggplot(data) +
74 ggplot2::geom_point(ggplot2::aes_string(x = data[, 2], y = data[, 3], color = data[, titre]), shape = "square", size = 2) + ggplot2::scale_colour_gradient(low = "blue", high = "orange", na.value = "grey50") +
75 ggplot2::xlab("Longitude") + ggplot2::ylab("Latitude") +
76 ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 90, vjust = 0.5, hjust = 1), plot.title = ggplot2::element_text(color = "black", size = 12, face = "bold")) + ggplot2::ggtitle(titre)
77
78 ggplot2::ggsave(paste0(titre, ".png"), graph_indices, width = 12, height = 10, units = "cm")
79 }