The function to be used to import PRISMA L1 data is pr_convert.
It takes as input the full path of a PRISMA L1 hdf5 image, an output folder name and format, and a series of switches which allow deciding which hyperspectral cubes and ancillary datasets should be crated.
In particular: - the VNIR and SWIR logical arguments allow deciding if importing the VNIR and SWIR hyperspectral cubes; - the FULL logical argument allows deciding if a complete VNIR+SWIR cube has to be created.
In that case (FULL = TRUE): - the join_priority keyword is used to decide if keeping bands from the “VNIR” or the “SWIR” data cube in the wavelengths’ region where they overlap; - the PAN, LATLON, CLOUD, GLINT, LC and ERR_MATRIX logical arguments allow deciding which of the corresponding ancillary datasets should be created (see the PRISMA manual for additional info); - the apply_errmatric logical argument allows deciding if pixels for which the ERR_MATRIX values are above one are to be set to NA; - the base_georef logical argument allows deciding if a “base” georeferencing in Lat/Lon WGS-84 based on the “GLT and Bowtie Correction” technique used in ENVI and described here is applied (if set to FALSE, the original 1000 x 1000 datasets are returned – flipped to orient them “north/south” – without projection.
All logical arguments are set to FALSE by default (with the exception of base_georef), allowing complete customisation of a pr_convert run.
For example, the following code accesses the input L1 file and saves the VNIR and SWIR cubes and the PAN, ANGLES and CLOUD datasets. See pr_convert() documentation for info on available arguments.
IMPORTANT NOTE: to run this, you’d need to download the example data from GitHub. The data is about 1 GB, so it could take some time. The data would be placed in subfolder testdata of the prismaread installation folder.
l1_he5_path <- file.path(
system.file("testdata", package = "prismaread"),
"PRS_L1_STD_OFFL_20200524103704_20200524103708_0001.he5"
)
# The example file is downloaded from Github if not already downloaded.
# WARNING: this may need a long time (1 GB).
if (!file.exists(l1_he5_path)) {
message("Downloading test data - this may need a long time...")
# Install package piggyback if missing - for downloaing GitHub attachments
if (!requireNamespace("piggyback", quietly = TRUE)) {
install.packages("piggyback")
}
l1_zip_path <- file.path(
system.file("testdata", package = "prismaread"),
"PRS_L1_STD_OFFL_20200524103704_20200524103708_0001.zip"
)
piggyback::pb_download(
"PRS_L1_STD_OFFL_20200524103704_20200524103708_0001.zip",
repo = "irea-cnr-mi/prismaread",
dest = dirname(l1_zip_path)
)
unzip(l1_zip_path, exdir = dirname(l1_he5_path))
unlink(l1_zip_path)
}
l1_out_dir <- file.path(tempdir(), "prismaread/L1")
dir.create(dirname(l1_out_dir), showWarnings = FALSE)
pr_convert(
in_file = l1_he5_path,
out_folder = l1_out_dir,
out_format = "GTiff",
VNIR = TRUE,
SWIR = TRUE,
LATLON = TRUE,
PAN = TRUE,
CLOUD = TRUE
)Output files are saved in l1_out_dir according to conventions described here and can then be accessed, visualised and processed using standard R syntax (e.g., with raster or stars).
list.files(l1_out_dir)
#> [1] "PRS_L1_STD_OFFL_20200524103704_20200524103708_0001_HCO_CLD.tif"
#> [2] "PRS_L1_STD_OFFL_20200524103704_20200524103708_0001_HCO_LATLON.tif"
#> [3] "PRS_L1_STD_OFFL_20200524103704_20200524103708_0001_HCO_PAN.tif"
#> [4] "PRS_L1_STD_OFFL_20200524103704_20200524103708_0001_HCO_SWIR.tif"
#> [5] "PRS_L1_STD_OFFL_20200524103704_20200524103708_0001_HCO_SWIR.tif.aux.xml"
#> [6] "PRS_L1_STD_OFFL_20200524103704_20200524103708_0001_HCO_SWIR.wvl"
#> [7] "PRS_L1_STD_OFFL_20200524103704_20200524103708_0001_HCO_VNIR.tif"
#> [8] "PRS_L1_STD_OFFL_20200524103704_20200524103708_0001_HCO_VNIR.wvl"
#> [9] "PRS_L1_STD_OFFL_20200524103704_20200524103708_0001_HCO.ang"
l1_vnir_path <- file.path(
l1_out_dir,
"PRS_L1_STD_OFFL_20200524103704_20200524103708_0001_HCO_VNIR.tif"
)
l1_vnir <- raster::brick(l1_vnir_path)
l1_vnir
#> class : RasterBrick
#> dimensions : 1279, 1265, 1617935, 63 (nrow, ncol, ncell, nlayers)
#> resolution : 0.0003852844, 0.0002593994 (x, y)
#> extent : 8.287386, 8.774771, 45.21347, 45.54524 (xmin, xmax, ymin, ymax)
#> crs : +proj=longlat +datum=WGS84 +no_defs
#> source : /tmp/Rtmp8H1qXw/prismaread/L1/PRS_L1_STD_OFFL_20200524103704_20200524103708_0001_HCO_VNIR.tif
#> names : PRS_L1_ST//HCO_VNIR.1, PRS_L1_ST//HCO_VNIR.2, PRS_L1_ST//HCO_VNIR.3, PRS_L1_ST//HCO_VNIR.4, PRS_L1_ST//HCO_VNIR.5, PRS_L1_ST//HCO_VNIR.6, PRS_L1_ST//HCO_VNIR.7, PRS_L1_ST//HCO_VNIR.8, PRS_L1_ST//HCO_VNIR.9, PRS_L1_ST//CO_VNIR.10, PRS_L1_ST//CO_VNIR.11, PRS_L1_ST//CO_VNIR.12, PRS_L1_ST//CO_VNIR.13, PRS_L1_ST//CO_VNIR.14, PRS_L1_ST//CO_VNIR.15, ...
#> min values : 0.00, 0.00, 0.00, 34.07, 30.83, 33.53, 31.60, 26.88, 26.54, 21.90, 19.04, 17.68, 10.70, 24.78, 31.15, ...
#> max values : 444.13, 490.64, 503.68, 497.34, 556.07, 635.39, 655.35, 655.35, 655.35, 655.35, 655.35, 655.35, 655.35, 655.35, 493.11, ...
mapview::viewRGB(l1_vnir, 40,30,20)The function also exports in text files ancillary data related to: 1. wavelengths and fwhms of the different images; 2. hour and sun geometry at acquisition.
l1_vnir_wvl <- read.table(
file.path(
l1_out_dir,
"PRS_L1_STD_OFFL_20200524103704_20200524103708_0001_HCO_VNIR.wvl"
),
header = TRUE
)
DT::datatable(l1_vnir_wvl)Associating acquisition angles with L1 data
PRISMA L1 data unfortunately does not contain information concerning acquisition angles, that is instead available for all L2 datasets. However, if both the L1 and any L2 dataset area available, prismaread allows associating the ANGLES data retrieved from the L2 dataset to the L1 one.
To do that, the user has to specify the additional argument in_L2_file in the call to pr_convert(). Note that, in this case, also the georeferencing information used for the GLT georeferencing is taken from the L2 dataset.
l2_he5_path <- file.path(
system.file("testdata", package = "prismaread"),
"PRS_L2C_STD_20200524103704_20200524103708_0001.he5"
)
# The example file is downloaded from Github if not already downloaded.
# WARNING: this may need a long time (1 GB).
if (!file.exists(l2_he5_path)) {
message("Downloading test data - this may need a long time...")
l2_zip_path <- file.path(
system.file("testdata", package = "prismaread"),
"PRS_L2C_STD_20200524103704_20200524103708_0001.zip"
)
piggyback::pb_download(
"PRS_L2C_STD_20200524103704_20200524103708_0001.zip",
repo = "irea-cnr-mi/prismaread",
dest = dirname(l2_zip_path)
)
unzip(l2_zip_path, exdir = dirname(l2_zip_path))
unlink(l2_zip_path)
}
l12_out_dir <- file.path(tempdir(), "prismaread/L12")
# Save only VNIR, also including the ANGLES dataset, in ENVI format,
# taking ANGLES and georeferencing info from a L2 file.
pr_convert(
in_file = l1_he5_path,
in_L2_file = l2_he5_path,
out_folder = l12_out_dir,
out_format = "ENVI",
VNIR = TRUE,
ANGLES = TRUE
)Importing only selected bands
Arguments selbands_vnir and selbands_swir allow selecting a specified subset of PRISMA bands, by specifying an array of required wavelengths. For example, the code below creates a 3-band VNIR cube, a 2-band SWIR and a 5-band FULL dataset, by selecting the original PRISMA bands with wavelengths closer to the requested ones.
l1sel_out_dir <- file.path(tempdir(), "prismaread/L1sel")
pr_convert(
in_file = l1_he5_path,
out_folder = l1sel_out_dir,
out_format = "GTiff",
VNIR = TRUE,
SWIR = TRUE,
FULL = TRUE,
selbands_vnir = c(450,550,650,750),
selbands_swir = c(1000,1330)
)
list.files(l1sel_out_dir)
#> [1] "PRS_L1_STD_OFFL_20200524103704_20200524103708_0001_HCO_FULL.tif"
#> [2] "PRS_L1_STD_OFFL_20200524103704_20200524103708_0001_HCO_FULL.wvl"
#> [3] "PRS_L1_STD_OFFL_20200524103704_20200524103708_0001_HCO_SWIR.tif"
#> [4] "PRS_L1_STD_OFFL_20200524103704_20200524103708_0001_HCO_SWIR.wvl"
#> [5] "PRS_L1_STD_OFFL_20200524103704_20200524103708_0001_HCO_VNIR.tif"
#> [6] "PRS_L1_STD_OFFL_20200524103704_20200524103708_0001_HCO_VNIR.wvl"
#> [7] "PRS_L1_STD_OFFL_20200524103704_20200524103708_0001_HCO.ang"
l1sel_vnir <- raster::brick(
file.path(
l1sel_out_dir,
"PRS_L1_STD_OFFL_20200524103704_20200524103708_0001_HCO_VNIR.tif"
)
)
l1sel_vnir
#> class : RasterBrick
#> dimensions : 1279, 1265, 1617935, 4 (nrow, ncol, ncell, nlayers)
#> resolution : 0.0003852844, 0.0002593994 (x, y)
#> extent : 8.287386, 8.774771, 45.21347, 45.54524 (xmin, xmax, ymin, ymax)
#> crs : +proj=longlat +datum=WGS84 +no_defs
#> source : /tmp/Rtmp8H1qXw/prismaread/L1sel/PRS_L1_STD_OFFL_20200524103704_20200524103708_0001_HCO_VNIR.tif
#> names : PRS_L1_STD_OFFL_20200524103704_20200524103708_0001_HCO_VNIR.1, PRS_L1_STD_OFFL_20200524103704_20200524103708_0001_HCO_VNIR.2, PRS_L1_STD_OFFL_20200524103704_20200524103708_0001_HCO_VNIR.3, PRS_L1_STD_OFFL_20200524103704_20200524103708_0001_HCO_VNIR.4
#> min values : 31.60, 28.61, 8.16, 5.96
#> max values : 655.35, 551.45, 655.35, 551.30
raster::plotRGB(l1sel_vnir, 4,3,2, stretch = "lin")
l1sel_vnir_wvl <- read.table(
file.path(
l1sel_out_dir,
"PRS_L1_STD_OFFL_20200524103704_20200524103708_0001_HCO_VNIR.wvl"
),
header = TRUE
)
DT::datatable(l1sel_vnir_wvl)Creation of ATCOR files
When working on L1 data, function pr_convert() also allows automatic creation of text files required to run an atmospheric correction using ATCOR. Those files are saved in the “ATCOR” subfolder of the main output folder.
In default behaviour, only the three “standard” ATCOR files (.wvl, .dat and .cal) are created in the “ATCOR” subfolder of the main output folder, with the .wvl file containing nominal wavelengths and FWHMs derived from the cw and fwhm attributes of the .he5 file.
For example, the code below creates input files for ATCOR useful for correction of the full hyperspectral cube, and places them in the “ATCOR” subfolder of the main output folder.
# Save a full image, prioritising the VNIR spectrometer and saving in ENVI format
pr_convert(
in_file = l1_he5_path,
out_folder = l1_out_dir,
out_format = "GTiff",
join_priority = "VNIR",
ATCOR = TRUE,
FULL = TRUE,
PAN = TRUE,
CLOUD = TRUE
)The user can also choose to generate additional ATCOR files, containing data about wavelengths and FWHMs related to different columns of the data cube, as derived from the KDP_AUX/Cw_Vnir_Matrix, KDP_AUX/Cw_Swir_Matrix, KDP_AUX/Cw_Fwhm_Matrix, KDP_AUX/Cw_Fwhm_Matrix HDF layers. This could allow running different atmospheric corrections for different columns of the data, potentially allowing compensating “smile” effects on the retrieved surface reflectances.
For example, the code below creates additional ATCOR input files with wavelengths corresponding to those of the columns 200 and 800.
l1sel2_out_dir <- file.path(tempdir(), "prismaread/L1sel2")
# Save a full image, prioritising the VNIR spectrometer and saving in ENVI format
pr_convert(
in_file = l1_he5_path,
out_folder = l1sel2_out_dir,
out_format = "ENVI",
join_priority = "SWIR",
ATCOR = TRUE,
ATCOR_wls = c(200,800),
FULL = TRUE,
PAN = TRUE,
CLOUD = TRUE
)IMPORTANT NOTE: the latter functionality is only appliable to “HRC” L1 data cubes.