0
|
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 }
|