diff --git a/vignettes/ars.Rmd b/vignettes/ars.Rmd deleted file mode 100644 index cfc2a2b..0000000 --- a/vignettes/ars.Rmd +++ /dev/null @@ -1,130 +0,0 @@ ---- -title: "ARS - Spectral Unmixing" -author: "Mike Cecil" -date: "2023-02-26" -output: html_document ---- - -```{r setup, include=FALSE} -knitr::opts_chunk$set(echo = TRUE) -``` - -Install `RStoolbox` from source. - -```{r} -library(devtools) -install_github("bleutner/RStoolbox") -``` - -This may throw an error, in which case install any other dependencies then rerun above. - -```{r} -#install.packages("", "", "") ## install dependencies -``` - -Load other libraries -```{r cars} -library(dplyr) -library(raster) -library(sf) -library(RStoolbox) -``` - -Set up file paths -```{r} -taxonomic_path <- "C:/Users/micha/Downloads/Endmembers/TaxonomicEndmembers.shp" -biophysical_path <- "C:/Users/micha/Downloads/Endmembers/BiophysicalEndmembers.shp" -landsat_path <- "C:/Users/micha/Downloads/Landsat8SR_SD_July2019_scale.tif" -``` - -Read in data -```{r} -tax_classes <- st_read (taxonomic_path) -#bio_classes <- st_read(biophysical_path) -lsat <- raster::brick(landsat_path) -``` - -Plot data -```{r} -plot(lsat[[4]]) -plot(tax_classes, add = T) -``` -Extract pixel values at point locations. - -```{r} -tax_em <- extract(lsat, tax_classes %>% st_centroid()) %>% - cbind(tax_classes, .) %>% as_tibble() -tax_em -``` -Because we have two points per land cover class, we will use `group_by` to average -```{r} -tax_group <- tax_em %>% group_by(class) %>% - summarise(mean_b1 = mean(Landsat8SR_SD_July2019_scale_1), - mean_b2 = mean(Landsat8SR_SD_July2019_scale_2), - mean_b3 = mean(Landsat8SR_SD_July2019_scale_3), - mean_b4 = mean(Landsat8SR_SD_July2019_scale_4), - mean_b5 = mean(Landsat8SR_SD_July2019_scale_5)) %>% - st_drop_geometry() -tax_group -``` - -We need our spectral table to only have a column for each band. So we'll put the land cover classes into row names (instead of as a column.) - -```{r} - -class_names <- tax_group$class - -tax_group <- tax_group %>% data.frame() %>% dplyr::select(-class) -rownames(tax_group) <- class_names -tax_group -``` - -Perform spectral unmixing. The issue with this algorithm is that the class percents may not add to 1. (this may take 10-20 min) - -```{r} - -probs_tax <- mesma(lsat, - tax_group, - method = "NNLS", - iterate = 1000) - -plot(probs_tax) -``` - -The code below is used to "rescale" the probabilities so they add to 1. - -```{r} -sum_tax <- sum(probs_tax) - probs_tax[['RMSE']] -probs_tax_adj <- probs_tax -## for loop rescales each band -for(band in class_names){ - # print(band) - probs_tax_adj[[band]] <- probs_tax_adj[[band]]/sum_tax -} - -summary(probs_tax_adj) - - -``` - -Let's check that the adjusted values sum to 1. - -```{r} -sum_tax_adj <- sum(probs_tax_adj) - probs_tax_adj[['RMSE']] ## sum all bands except RMSE -summary(sum_tax_adj) -plot(sum_tax_adj) -``` - -Now let's plot the final class percents - -```{r} -plot(probs_tax_adj, zlim = c(0,1)) -``` - -And we can save the results to a tif file. - -```{r} -output_folder <- "C:/Users/micha/Downloads/" -writeRaster(probs_tax_adj, file = paste0(output_folder, "class_probabilities.tif")) -``` -