## ----load_packages------------------------------------------------------- library(OpenImageR) library(tidyverse) ## ----load_image_example-------------------------------------------------- path <- "../donnees/dataset-parkinson/spiral/training/healthy/V01HE02.png" image <- imager::load.image(path) class(image) image ## ----structure_image----------------------------------------------------- str(image) ## ----image_example_plot_healthy------------------------------------------ plot(image, main = "Healthy subject") ## ----image_example_plot_parkinson---------------------------------------- path <- "../donnees/dataset-parkinson/spiral/training/parkinson/V03PE05.png" image <- imager::load.image(path) plot(image, main = "Subject with Parkinson's disease") ## ----image_example_grayscale--------------------------------------------- image_grey <- imager::grayscale(image) ## ----image_example_rm_transparency--------------------------------------- if(dim(image)[4] > 3) image <- imager::channel(image, 1:3) dim(image_grey) ## ----create_image_test--------------------------------------------------- image_test <- imager::load.example("parrots") image_test plot(image_test) ## ----plot_larger_image_test---------------------------------------------- larger_image_test <- imager::resize(im = image_test, size_x = 1200, size_y = 512, interpolation_type = 1) larger_image_test_no <- imager::resize(im = image_test, size_x = 1200, size_y = 512, interpolation_type = 0) op <- par(mfrow = c(1,2)) plot(larger_image_test, main = "Nearest-Neighbor interpolation") plot(larger_image_test_no, main = "Without interpolation") ## ----image_example_resize------------------------------------------------ image_resized <- imager::resize(image_grey, size_x = 256, size_y = 256) ## ----create_load_image_gs------------------------------------------------ #' load_image_gs #' For a given image, reads it, then return the grey-scale values #' for each pixel #' @param path path to the image load_image_gs <- function(path, size_x, size_y){ image <- imager::load.image(path) if(dim(image)[4] > 3) image <- imager::channel(image, 1:3) image_grey <- imager::grayscale(image) # Resize the image image_resized <- imager::resize(image_grey, size_x = size_x, size_y = size_y) image_resized } ## ----test_load_image_gs-------------------------------------------------- path <- "../donnees/dataset-parkinson/spiral/training/parkinson/V03PE05.png" test <- load_image_gs(path, size_x = 256, size_y = 256) plot(test) ## ----path_to_spiral_training_healthy------------------------------------- folder_path <- "../donnees/dataset-parkinson/spiral/training/healthy/" files <- list.files(folder_path, full.names = TRUE, pattern = "\\.png") head(files) ## ----create_images------------------------------------------------------- images <- lapply(files, load_image_gs, size_x = 256, size_y = 256) ## ------------------------------------------------------------------------ as.vector(images[[1]]) %>% head() ## ----create_values------------------------------------------------------- values <- lapply(images, function(x) as.vector(x)) %>% do.call("cbind", .) %>% t() dim(values) ## ----create_load_images_sample------------------------------------------- #' load_images_sample #' Computes the histogram of oriented gradient for the images #' of the spiral or waves, for either healthy or patients with Parkinson's disease #' in either the training or the testing set #' The result is a matrix whose lines correspond to pictures and columns correspond to extracted values #' @param drawing_type type of the drawing (`spiral` or `wave`) #' @param sample_type sample type (`training` or `testing`) #' @param patient_condition patient condition (`healthy` or `parkinson`) #' drawing_type <- "wave" ; sample_type <- "training" ; patient_condition <- "healthy" load_images_sample <- function(drawing_type, sample_type, patient_condition){ folder_path <- str_c("../donnees/dataset-parkinson/", drawing_type, "/", sample_type, "/", patient_condition) files <- list.files(folder_path, full.names = TRUE, pattern = "\\.png") size_x <- ifelse(drawing_type == "wave", yes = 512, no = 256) images <- lapply(files, load_image_gs, size_x = size_x, size_y = 256) values <- lapply(images, function(x) as.vector(x)) %>% do.call("cbind", .) %>% t() values } ## ----create_spiral_training_sample--------------------------------------- spiral_training_sample <- load_images_sample(drawing_type = "spiral", sample_type = "training", patient_condition = "healthy") %>% rbind( load_images_sample(drawing_type = "spiral", sample_type = "training", patient_condition = "parkinson") ) dim(spiral_training_sample) ## ----create spiral_training_sample_y------------------------------------- spiral_training_sample_y <- c(rep("healthy", 36), rep("parkinson", 36)) ## ----create_spiral_testing_sample---------------------------------------- # Test sample for the spiral images spiral_testing_sample <- load_images_sample(drawing_type = "spiral", sample_type = "testing", patient_condition = "healthy") %>% rbind( load_images_sample(drawing_type = "spiral", sample_type = "testing", patient_condition = "parkinson") ) dim(spiral_testing_sample) spiral_testing_sample_y <- c(rep("healthy", 15), rep("parkinson", 15)) ## ----create_pca_spiral--------------------------------------------------- pca_spiral <- prcomp(spiral_training_sample) ## ----pca_spiral_eigenvalues---------------------------------------------- pca_spiral$sdev^2 ## ----pca_spiral_eigenvectors--------------------------------------------- # pca_spiral$rotation ## ----create_coordinates_pca_spiral--------------------------------------- library(factoextra) coordinates_pca_spiral <- predict(pca_spiral) ## ----create_coordinates_pca_spiral_test---------------------------------- coordinates_pca_spiral_test <- predict(pca_spiral, newdata = spiral_testing_sample) ## ----create_pca_spiral_results_df---------------------------------------- pca_spiral_results_df <- tibble( factor_1 = coordinates_pca_spiral[,1], factor_2 = coordinates_pca_spiral[,2], y = as.factor(spiral_training_sample_y), type = "spiral") pca_spiral_results_df ## ------------------------------------------------------------------------ p_1 <- ggplot(data = pca_spiral_results_df, aes(x = factor_1, y = factor_2, colour = y)) + geom_point() + scale_colour_discrete("True value") + labs(x = "Factor 1", y = "Factor 2") + stat_ellipse(level = .95, type = "t", segments = 250) p_1 ## ----create_df_spiral---------------------------------------------------- df_spiral <- data.frame( y = as.factor(spiral_training_sample_y), x = coordinates_pca_spiral[,1:3]) %>% tbl_df() df_spiral ## ----create_df_test_spiral----------------------------------------------- df_test_spiral <- data.frame( y = as.factor(spiral_testing_sample_y), x = coordinates_pca_spiral_test[,1:3]) %>% tbl_df() df_test_spiral ## ----model_estimation---------------------------------------------------- library(nnet) reg <- multinom(y ~ ., data = df_spiral, trace = FALSE) ## ----df_spiral_predict--------------------------------------------------- df_spiral$pred <- predict(reg, type = "probs") df_spiral$pred_class <- predict(reg, type = "class") df_spiral ## ----create_confusion_matrix--------------------------------------------- confusion_matrix <- table(df_spiral$y, df_spiral$pred_class) confusion_matrix ## ----df_test_spiral_predict---------------------------------------------- df_test_spiral$pred <- predict(reg, newdata = df_test_spiral, type="probs") df_test_spiral$pred_class <- predict(reg, newdata = df_test_spiral, type="class") ## ----create_confusion_matrix_test---------------------------------------- confusion_matrix_test <- table(df_test_spiral$y, df_test_spiral$pred_class) confusion_matrix_test ## ----define_conf_mat_metrics--------------------------------------------- #' conf_mat_metrics #' Extract metrics from a confusion matrix #' @param conf_mat confusion matrix (2x2 matrix) #' @param ind_positive index of the column and row to consider as "positive" conf_mat_metrics <- function(conf_mat, ind_positive = 2){ ind_negative <- ifelse(ind_positive == 1, yes = 2, no = 1) TP <- conf_mat[ind_positive, ind_positive] TN <- conf_mat[ind_negative, ind_negative] FP <- conf_mat[ind_negative, ind_positive] FN <- conf_mat[ind_positive, ind_negative] Pos <- sum(conf_mat[ind_positive, ]) Neg <- sum(conf_mat[ind_negative, ]) # Accuracy acc <- (TP+TN)/(Pos+Neg) # Missclassification rate missclass_rate <- 1-acc # Sensitivity (True positive rate) # proportion of actual positives that are correctly identified as such TPR <- TP/Pos # Specificity (True negative rate) # proportion of actual negatives that are correctly identified as such TNR <- TN/Neg tibble(accuracy = acc, missclass_rate = missclass_rate, sensitivity = TPR, specificity = TNR) } ## ----display_confusion_matrices_metrics---------------------------------- conf_mat_metrics(conf_mat = confusion_matrix, ind_positive = 2) conf_mat_metrics(conf_mat = confusion_matrix_test, ind_positive = 2) ## ----define_classif_performance_pcs-------------------------------------- #' classif_performance_pcs #' Train a log-linear model on the patient condition #' using a subet of the first principal components only #' @param nb_pcs number of principal components to keep classif_performance_pcs <- function(nb_pcs){ # Data frame with the true value and the first k components of the PCA df <- data.frame( y = as.factor(spiral_training_sample_y), x = coordinates_pca_spiral[,1:nb_pcs]) %>% tbl_df() df_test <- data.frame( y = as.factor(spiral_testing_sample_y), x = coordinates_pca_spiral_test[,1:nb_pcs]) %>% tbl_df() # Regression of the values of the images on the first k components of the PCA reg <- multinom(y ~ ., data = df, trace = FALSE) # Training sample # Adding the predicted value and the predicted class to the data frame df$pred <- predict(reg, type="probs") df$pred_class <- predict(reg, type="class") # And then compute the contingency matrix confusion_matrix <- table(df$y, df$pred_class) # Testing sample df_test$pred <- predict(reg, newdata = df_test, type="probs") df_test$pred_class <- predict(reg, newdata = df_test, type="class") # Contingency table confusion_matrix_test <- table(df_test$y, df_test$pred_class) ind_positive <- which(colnames(confusion_matrix) == "parkinson") conf_mat_metrics(conf_mat = confusion_matrix, ind_positive = ind_positive) %>% mutate(sample_type = "training") %>% bind_rows( conf_mat_metrics(conf_mat = confusion_matrix_test, ind_positive = ind_positive) %>% mutate(sample_type = "testing") ) %>% mutate(nb_pcs = !!nb_pcs) }# End of classif_performance_pcs() ## ----example_classif_performance_pcs_3----------------------------------- classif_performance_pcs(3) ## ----create_missclass_rate----------------------------------------------- performances <- pbapply::pblapply(1:50, classif_performance_pcs) %>% bind_rows() performances ## ----create_p_2_--------------------------------------------------------- p_2 <- performances %>% ggplot(data = ., aes(x = nb_pcs, y = missclass_rate, colour = sample_type)) + geom_line() + geom_point() + labs(x = "Number of Factors", y = "Misclassification Rate") + scale_color_manual("Sample", values = c("training" = "orange", "testing" = "blue")) cowplot::plot_grid(p_1, p_2) ## ----facet_plot_metrics-------------------------------------------------- performances %>% gather(metrics, value, -sample_type, -nb_pcs) %>% ggplot(data = ., aes(x = nb_pcs, y = value, colour = sample_type)) + geom_line() + geom_point() + facet_wrap(~metrics) + labs(x = "Number of Factors", y = NULL) + scale_color_manual("Sample", values = c("training" = "orange", "testing" = "blue")) ## ----load_packages_satellite--------------------------------------------- library(tidyverse) library(raster) library(rgeos) library(rgdal) library(RStoolbox) ## ----list_files_2018072601T1--------------------------------------------- list.files("../donnees/landsat/LC080450332018072601T1-SC20190704113634/") ## ----create_files_bands_20180726----------------------------------------- files_bands_20180726 <- list.files("../donnees/landsat/LC080450332018072601T1-SC20190704113634/", pattern = "band[[:digit:]]{1}\\.tif$", full.names = TRUE) files_bands_20180726 ## ----create_bands_20180726_stack----------------------------------------- bands_20180726_stack <- raster::stack(files_bands_20180726) bands_20180726_stack ## ----coordinates_Mendocino----------------------------------------------- x_min <- 485000 x_max <- 550000 y_min <- 4310000 y_max <- 4370000 ## ----create_bounding_box_Mendocino--------------------------------------- p <- cbind(c(x_min, x_min, x_max, x_max), c(y_min, y_max, y_max, y_min)) %>% Polygon() %>% list() %>% Polygons(ID = "1") %>% list() %>% SpatialPolygons() p_df <- data.frame( ID=1:length(p) ) ## ----create_bounding_box_Mendocino_sp------------------------------------ p_sp <- SpatialPolygonsDataFrame(p, p_df) crs(p_sp) <- crs(bands_20180726_stack) plot(p_sp) ## ----plot_area_to_crop--------------------------------------------------- #' theme_map #' Theme for maps created using ggplot2 theme_map <- function(...) { theme_minimal() + theme( text = element_text(color = "#22211d"), axis.line = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), axis.title.x = element_blank(), axis.title.y = element_blank(), # panel.grid.minor = element_line(color = "#ebebe5", size = 0.2), panel.grid.major = element_line(color = "#ebebe5", size = 0.2), panel.grid.minor = element_blank(), plot.background = element_rect(fill = "#f5f5f2", color = NA), panel.background = element_rect(fill = "#f5f5f2", color = NA), legend.key = element_rect(colour = "black"), legend.background = element_rect(fill = "#f5f5f2", color = NA), panel.border = element_blank(), ... ) } # The landsat scene and the area of interest RStoolbox::ggRGB(bands_20180726_stack, r = 4, g = 3, b = 2, stretch = "lin") + # Highlight the area of interest geom_polygon(data = fortify(p_sp), aes(x = long, y = lat), fill = "red", alpha = .3, col = "red") + theme_map() + labs(title = "Landsat scene (2018-07-26)") ## ----create_bands_20180726_brick_crop------------------------------------ bands_20180726_brick_crop <- raster::crop(bands_20180726_stack, p_sp) bands_20180726_brick_crop ## ----plot_bands_20180726_brick_crop-------------------------------------- ggRGB(bands_20180726_brick_crop, r = 4, g = 3, b = 2, stretch = "lin") + theme_map() + labs(title = "Landsat scene (2018-07-26)", subtitle = "Perimeter of the Mendocino Complex fire") ## ----create_cloud_mask_20180726------------------------------------------ # Import the layer cloud_mask_20180726 <- raster("../donnees/landsat/LC080450332018072601T1-SC20190704113634/LC08_L1TP_045033_20180726_20180731_01_T1_pixel_qa.tif") # Extract the area of interest cloud_mask_20180726_crop <- crop(cloud_mask_20180726, p_sp) cloud_mask_20180726_crop ## ----create_cloud_mask_20180726_df--------------------------------------- cloud_mask_20180726_df <- as.data.frame(cloud_mask_20180726_crop, xy=T) %>% as_tibble() %>% # Rename the columns of the table magrittr::set_colnames(c("long", "lat", "pixel_qa")) %>% # Set the column `pixel_qa` as a factor mutate(pixel_qa = factor(pixel_qa)) cloud_mask_20180726_df ## ----table_cloud_mask_20180726_df_pixel_qa------------------------------- table(cloud_mask_20180726_df$pixel_qa) ## ----plot_cloud_mask_20180726_df----------------------------------------- ggplot(data = cloud_mask_20180726_df, aes(x = long, y = lat, fill = pixel_qa)) + geom_raster() + theme_map() + labs(title = "Landsat 2018-07-26 - Cloud mask/Water layer") ## ----cloud_mask_20180726_crop_remove_water------------------------------- cloud_mask_20180726_crop[ cloud_mask_20180726_crop == 324 ] <- NA as.data.frame(cloud_mask_20180726_crop, xy=TRUE) %>% as_tibble() %>% magrittr::set_colnames(c("long", "lat", "pixel_qa")) %>% mutate(pixel_qa = factor(pixel_qa)) %>% ggplot(data = ., aes(x = long, y = lat, fill = pixel_qa)) + geom_raster() + theme_map() + labs(title = "Landsat 2018-07-26 - Cloud/Water mask layer") ## ----cereate_bands_20180726_brick_crop_noclouds-------------------------- bands_20180726_brick_crop_noclouds <- raster::mask(bands_20180726_brick_crop, mask = cloud_mask_20180726_crop) ## ----plot_bands_20180726_brick_crop_noclouds----------------------------- ggRGB(bands_20180726_brick_crop_noclouds, r = 4, g = 3, b = 2, stretch = "lin") + labs(title = "Landsat 2018-07-26 - Without water") + theme_map() ## ----create_files_bands_20180827----------------------------------------- files_bands_20180827 <- list.files("../donnees/landsat/LC080450332018082701T1-SC20190704113725/", pattern = "band[[:digit:]]{+}\\.tif$", full.names = TRUE) ## ----create_bands_20180827_stack----------------------------------------- bands_20180827_stack <- raster::stack(files_bands_20180827) ## ------------------------------------------------------------------------ bands_20180827_brick_crop <- crop(bands_20180827_stack, p_sp) ## ------------------------------------------------------------------------ RStoolbox::ggRGB(bands_20180827_brick_crop, r = 4, g = 3, b = 2, stretch = "lin") + theme_map() + labs(title = "Landsat scene (2018-08-27)", subtitle = "Perimeter of the Mendocino Complex fire") ## ------------------------------------------------------------------------ cloud_mask_20180827 <- raster("../donnees/landsat/LC080450332018082701T1-SC20190704113725/LC08_L1TP_045033_20180827_20180911_01_T1_pixel_qa.tif") cloud_mask_20180827_crop <- crop(cloud_mask_20180827, p_sp) cloud_mask_20180827_df <- as.data.frame(cloud_mask_20180827_crop, xy=T) %>% as_tibble() %>% # Rename the columns of the table magrittr::set_colnames(c("long", "lat", "pixel_qa")) cloud_mask_20180827_df cloud_mask_20180827_df$pixel_qa %>% table() ## ----create_pixel_descriptions------------------------------------------- pixel_descriptions <- tibble(pixel_qa = c( 322, 386, 834, 324, 388, 328, 392, 352, 416, 480), description = c( rep("clear", 3), rep("water", 2), rep("cloud shadow", 2), rep("cloud", 3) )) cloud_mask_20180827_df <- cloud_mask_20180827_df %>% left_join( pixel_descriptions ) cloud_mask_20180827_df ## ----plot_cloud_mask_20180827_df----------------------------------------- ggplot(data = cloud_mask_20180827_df, aes(x = long, y = lat, fill = description)) + geom_raster() + theme_map() + labs(title = "Landsat 2018-08-27 - Cloud mask/Water layer") ## ----create_cloud_mask_20180827_crop_clouds------------------------------ mask_20180827_crop_clouds <- cloud_mask_20180827_crop mask_20180827_crop_clouds[ mask_20180827_crop_clouds %in% c(328, 392, 352, 416, 480) ] <- NA ## ----create_bands_20180827_brick_crop_noclouds--------------------------- bands_20180827_brick_crop_noclouds <- raster::mask(bands_20180827_brick_crop, mask = mask_20180827_crop_clouds) ## ----plot_create_bands_20180827_brick_crop_noclouds---------------------- ggRGB(bands_20180827_brick_crop_noclouds, r = 4, g = 3, b = 2, stretch = "lin") + labs(title = "Landsat 2018-08-27 - Without water") + theme_map() ## ----create_files_bands_20180912----------------------------------------- files_bands_20180912 <- list.files("../donnees/landsat/LC080450332018091201T1-SC20190704113621/", pattern = "band[[:digit:]]{+}\\.tif$", full.names = TRUE) bands_20180912_stack <- raster::stack(files_bands_20180912) # Cropping the raster bands_20180912_brick_crop <- crop(bands_20180912_stack, p_sp) # Plotting the image RStoolbox::ggRGB(bands_20180912_brick_crop, r = 4, g = 3, b = 2, stretch = "lin") + theme_map() + labs(title = "Landsat scene (2018-09-12)", subtitle = "Perimeter of the Mendocino Complex fire") ## ----create_cloud_mask_20180912------------------------------------------ cloud_mask_20180912 <- raster("../donnees/landsat/LC080450332018091201T1-SC20190704113621/LC08_L1TP_045033_20180912_20180927_01_T1_pixel_qa.tif") cloud_mask_20180912_crop <- crop(cloud_mask_20180912, p_sp) ## ----create_cloud_mask_20180912_df--------------------------------------- cloud_mask_20180912_df <- as.data.frame(cloud_mask_20180912_crop, xy=T) %>% as_tibble() %>% # Rename the columns of the table magrittr::set_colnames(c("long", "lat", "pixel_qa")) cloud_mask_20180912_df cloud_mask_20180912_df$pixel_qa %>% table() ## ----create_pixel_descriptions_2----------------------------------------- # Description table pixel_descriptions <- tibble(pixel_qa = c( 322, 386, 834, 898, 324, 388, 836, 328, 392, 836, 904, 352, 416, 480, 864, 928, 992), description = c( rep("clear", 4), rep("water", 3), rep("cloud shadow", 4), rep("cloud", 6) )) # Then joining the description to `cloud_mask_20180827_df` cloud_mask_20180912_df <- cloud_mask_20180912_df %>% left_join( pixel_descriptions ) cloud_mask_20180912_df ## ----plot_cloud_mask_20180912_df----------------------------------------- ggplot(data = cloud_mask_20180912_df, aes(x = long, y = lat, fill = description)) + geom_raster() + theme_map() + labs(title = "Landsat 2018-09-12 - Cloud mask/Water layer") ## ----create_bands_20180912_brick_crop_clear------------------------------ bands_20180912_brick_crop_clear <- mask(bands_20180912_brick_crop, mask = cloud_mask_20180912_crop) ggRGB(bands_20180912_brick_crop_clear, r = 4, g = 3, b = 2, stretch = "lin") + labs(title = "Landsat 2018-09-27 - Without water and clouds") + theme_map() ## ----create_post_fire_brick---------------------------------------------- post_fire_brick <- raster::cover(bands_20180827_brick_crop_noclouds, bands_20180912_brick_crop_clear) ## ----plot_post_fire_brick------------------------------------------------ ggRGB(post_fire_brick, r = 4, g = 3, b = 2, stretch = "lin") + labs(title = "Landsat 2018-09-27 - Without clouds") + theme_map() ## ----create_mask_20180827_crop_water------------------------------------- mask_20180827_crop_water <- cloud_mask_20180827_crop mask_20180827_crop_water[ mask_20180827_crop_water %in% c(324, 388, 836) ] <- NA ## ----post_fire_brick_rm_water-------------------------------------------- post_fire_brick <- raster::mask(post_fire_brick, mask = mask_20180827_crop_water) ## ----plot_post_fire_brick_without_water---------------------------------- ggRGB(post_fire_brick, r = 4, g = 3, b = 2, stretch = "lin") + labs(title = "Landsat 2018-08-27 - Without water") + theme_map() ## ----define_compute_nbr-------------------------------------------------- #' compute_nbr #' @param band_nir NIR band #' @param band_swir SWIR band compute_nbr <- function(band_nir, band_swir){ # this function calculates the difference between two rasters of the same CRS and extent # input: 2 raster layers of the same extent, crs that can be subtracted # output: a single different raster of the same extent, crs of the input rasters (band_nir - band_swir)/(band_nir + band_swir) } ## ----create_landsat_prefire_nbr------------------------------------------ landsat_prefire_nbr <- overlay(bands_20180726_brick_crop_noclouds[[5]], bands_20180726_brick_crop_noclouds[[7]], fun = compute_nbr) landsat_prefire_nbr landsat_postfire_nbr <- overlay(post_fire_brick[[5]], post_fire_brick[[7]], fun = compute_nbr) landsat_postfire_nbr ## ----plot_pre_post_fire_nbr---------------------------------------------- as.data.frame(landsat_prefire_nbr, xy=TRUE) %>% as_tibble() %>% rename("nbr" = "layer") %>% mutate(period = "pre") %>% bind_rows( as.data.frame(landsat_postfire_nbr, xy=TRUE) %>% as_tibble() %>% rename("nbr" = "layer") %>% mutate(period = "post") ) %>% mutate(period = factor(period, levels = c("pre", "post"), labels = c("Pre-fire NBR", "Post-fire NBR"))) %>% ggplot(data = ., aes(x = x, y = y, fill = nbr)) + geom_raster() + scale_fill_gradient2("NBR", low = "red", mid = "white", high = "darkgreen", limits = c(-1, 1)) + facet_wrap(~period) + theme_map() + coord_quickmap() ## ----create_landsat_diff_nbr--------------------------------------------- landsat_diff_nbr <- landsat_prefire_nbr - landsat_postfire_nbr ## ----create_landsat_postfire_dnbr_df------------------------------------- landsat_postfire_dnbr_df <- as.data.frame(landsat_diff_nbr, xy=TRUE) %>% as_tibble() %>% rename("dNBR" = "layer") %>% mutate( dNBR_class = cut(dNBR, c(-Inf, -.1, .1, .27, .66, Inf), right=TRUE, labels = c("Enhanced Regrowth", "Unburned", "Low Severity", "Moderate Severity", "High Severity")) ) landsat_postfire_dnbr_df ## ----create_colors_corresp----------------------------------------------- colours <- tibble(name = c("dark green", "green", "yellow", "orange", "red"), colour = c("#49934b", "#b1d778", "#fffdc6", "#f1b16e", "#c4352a"), dNBR_class = c("Enhanced Regrowth", "Unburned", "Low Severity", "Moderate Severity", "High Severity")) ## ----plot_landsat_postfire_dnbr_df--------------------------------------- ggplot(data = landsat_postfire_dnbr_df, aes(x = x, y = y, fill = dNBR_class)) + geom_raster() + coord_quickmap() + scale_fill_manual(NULL, values = colours$colour) + labs(x = NULL, y = NULL, title = "Landsat dNBR - Cold Spring fire site (July 2018)") + theme_map() ## ----landsat_postfire_dnbr_df_classes------------------------------------ landsat_postfire_dnbr_df$dNBR_class %>% table() ## ----landsat_postfire_dnbr_df_class_hist--------------------------------- landsat_postfire_dnbr_df %>% ggplot(data = ., aes(x = dNBR_class, fill = dNBR_class)) + geom_histogram(stat = "count") + scale_fill_manual(NULL, values = colours$colour) + labs(x = NULL) + theme(axis.text.x = element_text(angle = 45, hjust = 1)) ## ----create_nb_cell_high_severity---------------------------------------- nb_cell_high_severity <- landsat_postfire_dnbr_df %>% filter(dNBR_class %in% c("High Severity")) %>% nrow() # Squares kilometres nb_cell_high_severity * 30 * 30 / 1000000