From a0692a0f7791e1212a9d9a6ca576afedcc088154 Mon Sep 17 00:00:00 2001 From: lbusett Date: Thu, 17 Sep 2020 11:25:44 +0200 Subject: [PATCH] fix bug if index uses only one band + small improvements --- DESCRIPTION | 13 ++++++++----- R/pr_compute_indexes.R | 6 +++++- R/pr_convert.R | 15 +++------------ R/pr_create_vnir.R | 2 -- inst/extdata/indexes_list.txt | 2 ++ tests/testthat/test_pr_convert.R | 8 ++++---- vignettes/Importing-Level-1-Data.Rmd | 3 ++- vignettes/Importing-Level-2-Data.Rmd | 3 ++- .../{prismaread.Rmd => prismaread_vignette.Rmd} | 0 9 files changed, 26 insertions(+), 26 deletions(-) rename vignettes/{prismaread.Rmd => prismaread_vignette.Rmd} (100%) diff --git a/DESCRIPTION b/DESCRIPTION index a5b9aa5..1687124 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -7,11 +7,14 @@ Authors@R: comment = c(ORCID = '0000-0001-9634-6038')), person("Luigi", "Ranghetti", email = "ranghetti.l@irea.cnr.it", role = c("aut"), comment = c(ORCID = '0000-0001-6207-5188'))) -Description: `prismaread` allows easily importing PRISMA hyperspectral data (http://www.prisma-i.it/index.php/it/) - from the original data provided by ASI in HDF format, and convert them to a easier to use format (ENVI or GeoTiff). - It also provides functionality for automatically computing Spectral Indexes from either the original HDF data - or from hyperspectral data already converted using function `pr_convert`, and for easily and quickly extracting - data and computing statistics for the different bands over areas of interest. +Description: `prismaread` allows easily importing PRISMA (http://www.prisma-i.it/index.php/it/) + satellite hyperspectral data cubes and ancillary information from the original + data provided by ASI in HDF format, and convert them to a easier to use format + (ENVI or GeoTiff). + It also provides functionality for automatically computing Spectral Indexes from + either the original HDF data or from hyperspectral data already converted using + function `pr_convert`, and for easily and quickly extracting data and computing + statistics for the different bands over areas of interest. License: GPL-3 Encoding: UTF-8 LazyData: true diff --git a/R/pr_compute_indexes.R b/R/pr_compute_indexes.R index 939ba3b..806349f 100644 --- a/R/pr_compute_indexes.R +++ b/R/pr_compute_indexes.R @@ -142,7 +142,11 @@ pr_compute_indexes <- function(in_file, message("Writing Block: ", i, " of: ", bs$n) x <- raster::getValues(in_rast_comp, row = bs$row[i], nrows = bs$nrows[i]) - x <- eval(parse(text = indform)) + if (inherits(x, "numeric")) { + x <- x + } else { + x <- eval(parse(text = indform)) + } out <- raster::writeValues(out, x, bs$row[i]) } out <- raster::writeStop(out) diff --git a/R/pr_convert.R b/R/pr_convert.R index 84f3bc8..facc5cf 100644 --- a/R/pr_convert.R +++ b/R/pr_convert.R @@ -465,10 +465,6 @@ pr_convert <- function(in_file, wl_swir <- wl_swir[seqbands_swir] fwhm_swir <- fwhm_swir[seqbands_swir] } - # be sure to remove zeroes also if swir already present to avoid ---- - # errors on creation of FULL if recycling an existing cube from file - # fwhm_swir <- fwhm_swir[wl_swir != 0] - # wl_swir <- wl_swir[wl_swir != 0] # create FULL data cube and convert to raster ---- if (FULL & is.null(indexes)) { @@ -494,6 +490,7 @@ pr_convert <- function(in_file, } } + # Create and write the FULL hyperspectral cube if needed ---- if (FULL) { if (file.exists(out_file_full) & !overwrite) { message("FULL file already exists - use overwrite = TRUE or change @@ -556,7 +553,7 @@ pr_convert <- function(in_file, } } - + # Write ENVI header if needed ---- if (out_format == "ENVI") { cat("band names = {", paste(names(rast_tot),collapse=","), "}", "\n", file=raster::extension(out_file_full, "hdr"), append = TRUE) @@ -606,12 +603,6 @@ pr_convert <- function(in_file, overwrite = overwrite) } - # if (CONTREM) { - # - # # TODO: Add functionality for continuum removal computatation ? - # - # } - } # Save PAN if requested ---- @@ -789,7 +780,7 @@ pr_convert <- function(in_file, # second run: create indexes ----- # in this way we can use the same function, - # and when creating indexes create a temporary full raster, containgn only + # and when creating indexes create a temporary full raster, containing only # bands required for the selected indexes if (!is.null(indexes)) { diff --git a/R/pr_create_vnir.R b/R/pr_create_vnir.R index 508eb08..d4e5176 100644 --- a/R/pr_create_vnir.R +++ b/R/pr_create_vnir.R @@ -177,8 +177,6 @@ pr_create_vnir <- function(f, fwhm_sub <- fwhm_vnir[seqbands] } - # wl_vnir <- wl_vnir[wl_vnir != 0] - # fwhm_vnir <- fwhm_vnir[fwhm_vnir != 0] rm(vnir_cube) rm(band) gc() diff --git a/inst/extdata/indexes_list.txt b/inst/extdata/indexes_list.txt index 0613be1..90754f8 100644 --- a/inst/extdata/indexes_list.txt +++ b/inst/extdata/indexes_list.txt @@ -129,3 +129,5 @@ "128" "ZTM" "Zarco-Tejada & Miller" "R750 / R710" "Chlorophyll" "Zarco-Tejada et al. (2001)" "129" "myindex" "My custom Index" "R600 / R700" NA "Me (2020)" "130" "009066L" "My custom Index" "R600 / R700" NA "Me (2020)" +"131" "009449D" "My custom Index" "R600 / R700" NA "Me (2020)" +"132" "008530N" "My custom Index" "R600 / R700" NA "Me (2020)" diff --git a/tests/testthat/test_pr_convert.R b/tests/testthat/test_pr_convert.R index 34b035b..6f51f7d 100644 --- a/tests/testthat/test_pr_convert.R +++ b/tests/testthat/test_pr_convert.R @@ -356,19 +356,19 @@ test_that( vnir <- raster::brick( file.path(out_folder_L1, "testL1_1_HCO_VNIR.tif")) means_vnir <- as.numeric(raster::cellStats(vnir, "mean", na.rm = TRUE)) - testthat::expect_equal(means_vnir, c(71.72868 , 45.73275 , 58.20885), + testthat::expect_equal(means_vnir, c(71.74590 , 45.75219 , 58.21456), tolerance = 0.0001) swir <- raster::brick( file.path(out_folder_L1, "testL1_1_HCO_SWIR.tif")) means_swir <- as.numeric(raster::cellStats(swir, "mean", na.rm = TRUE)) - testthat::expect_equal(means_swir, c(8.921804), tolerance = 0.0001) + testthat::expect_equal(means_swir, c(8.926156), tolerance = 0.0001) full <- raster::brick( file.path(out_folder_L1, "testL1_1_HCO_FULL.tif")) means_full <- as.numeric(raster::cellStats(full, "mean", na.rm = TRUE)) - testthat::expect_equal(means_full, c(71.72868, 45.73275, 58.20885, - 8.921804), tolerance = 0.0001) + testthat::expect_equal(means_full, c(71.74590, 45.75219, 58.21456, + 8.926156), tolerance = 0.0001) # Use L2file to georeference ---- in_L2_file <- file.path( diff --git a/vignettes/Importing-Level-1-Data.Rmd b/vignettes/Importing-Level-1-Data.Rmd index 0d7d79d..de454c6 100644 --- a/vignettes/Importing-Level-1-Data.Rmd +++ b/vignettes/Importing-Level-1-Data.Rmd @@ -154,7 +154,8 @@ raster::plotRGB(r, 4,3,2, stretch = "lin") DT::datatable(read.table(file.path( out_folder, - "PRS_L1_STD_OFFL_20200524103704_20200524103708_0001_HCO_VNIR.wvl"))) + "PRS_L1_STD_OFFL_20200524103704_20200524103708_0001_HCO_VNIR.wvl")), + header = TRUE) ``` diff --git a/vignettes/Importing-Level-2-Data.Rmd b/vignettes/Importing-Level-2-Data.Rmd index 048b8ae..f04c1e8 100644 --- a/vignettes/Importing-Level-2-Data.Rmd +++ b/vignettes/Importing-Level-2-Data.Rmd @@ -79,7 +79,8 @@ The function also saves ancillary data related to wavelengths and fwhms of the d ```{r echo=TRUE, message=FALSE, warning=FALSE} wvls <- read.table(file.path( out_folder_L2D, - "PRS_L2D_STD_20200524103704_20200524103708_0001_HCO_VNIR.wvl")) + "PRS_L2D_STD_20200524103704_20200524103708_0001_HCO_VNIR.wvl"), + header = TRUE) DT::datatable(wvls) ``` diff --git a/vignettes/prismaread.Rmd b/vignettes/prismaread_vignette.Rmd similarity index 100% rename from vignettes/prismaread.Rmd rename to vignettes/prismaread_vignette.Rmd