Skip to content
GitLab
Menu
Projects
Groups
Snippets
Loading...
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in / Register
Toggle navigation
Menu
Open sidebar
Feret Jean-Baptiste
biodivMapR
Commits
c016a605
Commit
c016a605
authored
Jul 02, 2019
by
Florian de Boissieu
Browse files
clean code and update imports/exports for package
parent
bde6155f
Changes
86
Expand all
Hide whitespace changes
Inline
Side-by-side
DESCRIPTION
View file @
c016a605
...
...
@@ -26,5 +26,6 @@ Imports:
dissUtils,
labdsv,
raster,
rgdal
rgdal,
tools
RoxygenNote: 6.1.1
NAMESPACE
View file @
c016a605
...
...
@@ -3,7 +3,24 @@
S3method(split,line)
export(Check.Data.Format)
export(Convert.Raster2BIL)
export(Get.Diversity.From.Plots)
export(Get.List.Shp)
export(Get.Projection)
export(Map.Alpha.Diversity)
export(Map.Beta.Diversity)
export(Map.Spectral.Species)
export(Perform.PCA.Image)
export(Perform.Radiometric.Filtering)
export(Select.Components)
import(future)
import(raster)
import(tools)
importFrom(dissUtils,diss)
importFrom(fields,rdist)
importFrom(future,plan)
importFrom(future.apply,future_lapply)
importFrom(labdsv,pco)
importFrom(matlab,padarray)
importFrom(matrixStats,rowSds)
importFrom(rgdal,readOGR)
importFrom(snow,splitRows)
R/Lib_CheckConvertData.R
View file @
c016a605
...
...
@@ -12,138 +12,142 @@
#' converts a raster into BIL format as expected by DivMapping codes
#'
#' @param Raster.Path
f
ull path for the raster to be converted
#' @param Sensor
n
ame of the sensor
#' @param Convert.Integer
s
hould data be converted into integer ?
#' @param Multiplying.Factor
m
ultiplying factor (eg convert real reflectance values between 0 and 1 into integer between 0 and 10000)
#' @param Output.Directory
#' @param Multiplying.Factor.Last
m
ultiplying factor for last band
#' @param Raster.Path
character. F
ull path for the raster to be converted
#' @param Sensor
character. N
ame of the sensor
#' @param Convert.Integer
boolean. S
hould data be converted into integer ?
#' @param Multiplying.Factor
numeric. M
ultiplying factor (eg convert real reflectance values between 0 and 1 into integer between 0 and 10000)
.
#' @param Output.Directory
character. Path to output directory.
#' @param Multiplying.Factor.Last
numeric. M
ultiplying factor for last band
.
#'
#' @return Output.Path path for the image converted into ENVI BIL format
#' @import raster
#' @export
Convert.Raster2BIL
=
function
(
Raster.Path
,
Sensor
=
'
unknown
'
,
Output.Directory
=
FALSE
,
Convert.Integer
=
TRUE
,
Multiplying.Factor
=
1
,
Multiplying.Factor.Last
=
1
){
Convert.Raster2BIL
<-
function
(
Raster.Path
,
Sensor
=
"
unknown
"
,
Output.Directory
=
FALSE
,
Convert.Integer
=
TRUE
,
Multiplying.Factor
=
1
,
Multiplying.Factor.Last
=
1
)
{
# get directory and file name of original image
Input.File
=
basename
(
Raster.Path
)
Input.Dir
=
dirname
(
Raster.Path
)
Input.File
<-
basename
(
Raster.Path
)
Input.Dir
<-
dirname
(
Raster.Path
)
# define path where data will be stored
if
(
Output.Directory
==
FALSE
){
Output.Path
=
paste
(
Input.Dir
,
'
/Converted_DivMAP/
'
,
file_path_sans_ext
(
Input.File
),
sep
=
''
)
if
(
Output.Directory
==
FALSE
)
{
Output.Path
<-
paste
(
Input.Dir
,
"
/Converted_DivMAP/
"
,
file_path_sans_ext
(
Input.File
),
sep
=
""
)
}
else
{
dir.create
(
Output.Directory
,
showWarnings
=
FALSE
,
recursive
=
TRUE
)
Output.Path
=
paste
(
Output.Directory
,
'/'
,
file_path_sans_ext
(
Input.File
),
sep
=
''
)
dir.create
(
Output.Directory
,
showWarnings
=
FALSE
,
recursive
=
TRUE
)
Output.Path
<-
paste
(
Output.Directory
,
"/"
,
file_path_sans_ext
(
Input.File
),
sep
=
""
)
}
message
(
'
The converted file will be written in the following location:
'
)
message
(
"
The converted file will be written in the following location:
"
)
print
(
Output.Path
)
Output.Dir
=
dirname
(
Output.Path
)
dir.create
(
Output.Dir
,
showWarnings
=
FALSE
,
recursive
=
TRUE
)
Output.Dir
<-
dirname
(
Output.Path
)
dir.create
(
Output.Dir
,
showWarnings
=
FALSE
,
recursive
=
TRUE
)
# apply multiplying factors
message
(
'
reading initial file
'
)
Output.Img
=
Multiplying.Factor
*
brick
(
Raster.Path
)
Last.Band.Name
=
Output.Img
@
data
@
names
[
length
(
Output.Img
@
data
@
names
)]
Output.Img
[[
Last.Band.Name
]]
=
Multiplying.Factor.Last
*
Output.Img
[[
Last.Band.Name
]]
message
(
"
reading initial file
"
)
Output.Img
<-
Multiplying.Factor
*
brick
(
Raster.Path
)
Last.Band.Name
<-
Output.Img
@
data
@
names
[
length
(
Output.Img
@
data
@
names
)]
Output.Img
[[
Last.Band.Name
]]
<-
Multiplying.Factor.Last
*
Output.Img
[[
Last.Band.Name
]]
# convert into integer
if
(
Convert.Integer
==
TRUE
){
Output.Img
=
round
(
Output.Img
)
if
(
Convert.Integer
==
TRUE
)
{
Output.Img
<-
round
(
Output.Img
)
}
# write raster
message
(
'
writing converted file
'
)
if
(
Convert.Integer
==
TRUE
){
r
=
writeRaster
(
Output.Img
,
filename
=
Output.Path
,
format
=
"EHdr"
,
overwrite
=
TRUE
,
datatype
=
'
INT2S
'
)
message
(
"
writing converted file
"
)
if
(
Convert.Integer
==
TRUE
)
{
r
<-
writeRaster
(
Output.Img
,
filename
=
Output.Path
,
format
=
"EHdr"
,
overwrite
=
TRUE
,
datatype
=
"
INT2S
"
)
}
else
{
r
=
writeRaster
(
Output.Img
,
filename
=
Output.Path
,
format
=
"EHdr"
,
overwrite
=
TRUE
)
r
<-
writeRaster
(
Output.Img
,
filename
=
Output.Path
,
format
=
"EHdr"
,
overwrite
=
TRUE
)
}
hdr
(
r
,
format
=
"ENVI"
)
hdr
(
r
,
format
=
"ENVI"
)
# remove unnecessary files
File2Remove
=
paste
(
Output.Path
,
'
.aux.xml
'
,
sep
=
''
)
File2Remove2
=
paste
(
file_path_sans_ext
(
Output.Path
),
'
.aux.xml
'
,
sep
=
''
)
File2Remove
<-
paste
(
Output.Path
,
"
.aux.xml
"
,
sep
=
""
)
File2Remove2
<-
paste
(
file_path_sans_ext
(
Output.Path
),
"
.aux.xml
"
,
sep
=
""
)
file.remove
(
File2Remove
)
File2Remove
=
paste
(
Output.Path
,
'
.prj
'
,
sep
=
''
)
File2Remove2
=
paste
(
file_path_sans_ext
(
Output.Path
),
'
.prj
'
,
sep
=
''
)
if
(
file.exists
(
File2Remove
)){
File2Remove
<-
paste
(
Output.Path
,
"
.prj
"
,
sep
=
""
)
File2Remove2
<-
paste
(
file_path_sans_ext
(
Output.Path
),
"
.prj
"
,
sep
=
""
)
if
(
file.exists
(
File2Remove
))
{
file.remove
(
File2Remove
)
}
else
if
(
file.exists
(
File2Remove2
)){
}
else
if
(
file.exists
(
File2Remove2
))
{
file.remove
(
File2Remove2
)
}
File2Remove
=
paste
(
Output.Path
,
'
.sta
'
,
sep
=
''
)
File2Remove
=
paste
(
file_path_sans_ext
(
Output.Path
),
'
.sta
'
,
sep
=
''
)
if
(
file.exists
(
File2Remove
)){
File2Remove
<-
paste
(
Output.Path
,
"
.sta
"
,
sep
=
""
)
File2Remove
<-
paste
(
file_path_sans_ext
(
Output.Path
),
"
.sta
"
,
sep
=
""
)
if
(
file.exists
(
File2Remove
))
{
file.remove
(
File2Remove
)
}
else
if
(
file.exists
(
File2Remove2
)){
}
else
if
(
file.exists
(
File2Remove2
))
{
file.remove
(
File2Remove2
)
}
File2Remove
=
paste
(
Output.Path
,
'
.stx
'
,
sep
=
''
)
File2Remove2
=
paste
(
file_path_sans_ext
(
Output.Path
),
'
.stx
'
,
sep
=
''
)
if
(
file.exists
(
File2Remove
)){
File2Remove
<-
paste
(
Output.Path
,
"
.stx
"
,
sep
=
""
)
File2Remove2
<-
paste
(
file_path_sans_ext
(
Output.Path
),
"
.stx
"
,
sep
=
""
)
if
(
file.exists
(
File2Remove
))
{
file.remove
(
File2Remove
)
}
else
if
(
file.exists
(
File2Remove2
)){
}
else
if
(
file.exists
(
File2Remove2
))
{
file.remove
(
File2Remove2
)
}
File2Rename
=
paste
(
file_path_sans_ext
(
Output.Path
),
'
.hdr
'
,
sep
=
''
)
File2Rename2
=
paste
(
Output.Path
,
'
.hdr
'
,
sep
=
''
)
if
(
file.exists
(
File2Rename
)){
file.rename
(
from
=
File2Rename
,
to
=
File2Rename2
)
File2Rename
<-
paste
(
file_path_sans_ext
(
Output.Path
),
"
.hdr
"
,
sep
=
""
)
File2Rename2
<-
paste
(
Output.Path
,
"
.hdr
"
,
sep
=
""
)
if
(
file.exists
(
File2Rename
))
{
file.rename
(
from
=
File2Rename
,
to
=
File2Rename2
)
}
# change dot into underscore
Output.Path.US
=
gsub
(
Output.Path
,
pattern
=
'[.]'
,
replacement
=
'_'
)
if
(
!
Output.Path.US
==
Output.Path
){
file.rename
(
from
=
Output.Path
,
to
=
Output.Path.US
)
Output.Path.US
<-
file.path
(
dirname
(
Output.Path
),
gsub
(
basename
(
Output.Path
),
pattern
=
"[.]"
,
replacement
=
"_"
)
)
if
(
!
Output.Path.US
==
Output.Path
)
{
file.rename
(
from
=
Output.Path
,
to
=
Output.Path.US
)
}
Output.Path.US.HDR
=
paste
(
gsub
(
Output.Path
,
pattern
=
'[.]'
,
replacement
=
'_'
),
'.hdr'
,
sep
=
''
)
if
(
!
Output.Path.US.HDR
==
paste
(
Output.Path
,
'.hdr'
,
sep
=
''
)){
file.rename
(
from
=
paste
(
Output.Path
,
'.hdr'
,
sep
=
''
),
to
=
Output.Path.US.HDR
)
file.rename
(
from
=
Output.Path.US.HDR
,
to
=
Output.Path.US.HDR
)
Output.Path.US.HDR
<-
paste0
(
Output.Path.US
,
".hdr"
)
if
(
!
Output.Path.US.HDR
==
paste0
(
Output.Path
,
".hdr"
))
{
file.rename
(
from
=
paste0
(
Output.Path
,
".hdr"
),
to
=
Output.Path.US.HDR
)
### UTILITY?? ###
file.rename
(
from
=
Output.Path.US.HDR
,
to
=
Output.Path.US.HDR
)
}
if
(
!
Sensor
==
'
unknown
'
)
{
HDR.Temp.Path
=
system.file
(
'
extdata
'
,
'HDR'
,
paste0
(
Sensor
,
'
.hdr
'
),
package
=
'
DiversityMappR
'
)
if
(
file.exists
(
HDR.Temp.Path
)){
message
(
'
reading header template corresponding to the sensor located here:
'
)
if
(
!
Sensor
==
"
unknown
"
)
{
HDR.Temp.Path
<-
system.file
(
"
extdata
"
,
"HDR"
,
paste0
(
Sensor
,
"
.hdr
"
),
package
=
"
DiversityMappR
"
)
if
(
file.exists
(
HDR.Temp.Path
))
{
message
(
"
reading header template corresponding to the sensor located here:
"
)
print
(
HDR.Temp.Path
)
# get raster template corresponding to the sensor
HDR.Template
=
read.ENVI.header
(
HDR.Temp.Path
)
HDR.Template
<-
read.ENVI.header
(
HDR.Temp.Path
)
# get info to update hdr file
# read hdr
HDR.input
=
read.ENVI.header
(
Get.HDR.Name
(
Output.Path
))
if
(
!
is.null
(
HDR.Template
$
wavelength
)){
HDR.input
$
wavelength
=
HDR.Template
$
wavelength
HDR.input
<-
read.ENVI.header
(
Get.HDR.Name
(
Output.Path
))
if
(
!
is.null
(
HDR.Template
$
wavelength
))
{
HDR.input
$
wavelength
<-
HDR.Template
$
wavelength
}
if
(
!
is.null
(
HDR.Template
$
`sensor type`
)){
HDR.input
$
`sensor type`
=
HDR.Template
$
`sensor type`
if
(
!
is.null
(
HDR.Template
$
`sensor type`
))
{
HDR.input
$
`sensor type`
<-
HDR.Template
$
`sensor type`
}
if
(
!
is.null
(
HDR.Template
$
`band names`
)){
HDR.input
$
`band names`
=
HDR.Template
$
`band names`
if
(
!
is.null
(
HDR.Template
$
`band names`
))
{
HDR.input
$
`band names`
<-
HDR.Template
$
`band names`
}
if
(
!
is.null
(
HDR.Template
$
`wavelength units`
)){
HDR.input
$
`wavelength units`
=
HDR.Template
$
`wavelength units`
if
(
!
is.null
(
HDR.Template
$
`wavelength units`
))
{
HDR.input
$
`wavelength units`
<-
HDR.Template
$
`wavelength units`
}
# define visual stretch in the VIS domain
HDR.input
$
`default stretch`
=
'
0 1000 linear
'
HDR.input
$
`default stretch`
<-
"
0 1000 linear
"
# write corresponding hdr file
write.ENVI.header
(
HDR.input
,
Get.HDR.Name
(
Output.Path
))
}
else
if
(
!
file.exists
(
HDR.Temp.Path
)){
message
(
'Header template corresponding to the sensor expected to be found here'
)
}
else
if
(
!
file.exists
(
HDR.Temp.Path
))
{
message
(
"Header template corresponding to the sensor expected to be found here"
)
print
(
HDR.Temp.Path
)
message
(
'
please provide this header template in order to write info in HDR file
'
)
message
(
"
please provide this header template in order to write info in HDR file
"
)
print
(
Get.HDR.Name
(
Output.Path
))
message
(
'
or manually add wavelength location in HDR file, if relevant
'
)
message
(
"
or manually add wavelength location in HDR file, if relevant
"
)
}
}
else
if
(
Sensor
==
'
unknown
'
)
{
message
(
'
please make sure that the follozing header file contains information required
'
)
}
else
if
(
Sensor
==
"
unknown
"
)
{
message
(
"
please make sure that the follozing header file contains information required
"
)
print
(
Get.HDR.Name
(
Output.Path
))
message
(
'
or manually add wavelength location in HDR file, if relevant
'
)
message
(
"
or manually add wavelength location in HDR file, if relevant
"
)
}
return
(
Output.Path
)
}
...
...
@@ -153,78 +157,76 @@ Convert.Raster2BIL = function(Raster.Path,Sensor = 'unknown',Output.Directory =
#' @param Raster.Path full path for the raster to be converted
#' @param Mask is the raster a mask?
#'
#' @return Output.Path path for the image converted into ENVI BIL format
#' @export
Check.Data.Format
=
function
(
Raster.Path
,
Mask
=
FALSE
){
HDR.Path
=
Get.HDR.Name
(
Raster.Path
)
Check.Data.Format
<-
function
(
Raster.Path
,
Mask
=
FALSE
)
{
HDR.Path
<-
Get.HDR.Name
(
Raster.Path
)
# check if the hdr file exists
if
(
file.exists
(
HDR.Path
)){
HDR
=
read.ENVI.header
(
HDR.Path
)
if
(
Mask
==
FALSE
&
(
!
HDR
$
interleave
==
'
bil
'
)
&
(
!
HDR
$
interleave
==
'
BIL
'
)){
message
(
''
)
message
(
'
*********************************************************
'
)
message
(
'
The image format may not compatible with the processing chain
'
)
message
(
'
Image format expected:
'
)
message
(
'
ENVI hdr file with band interleaved by line (BIL) file format
'
)
message
(
''
)
message
(
'
Current Image format
'
)
if
(
file.exists
(
HDR.Path
))
{
HDR
<-
read.ENVI.header
(
HDR.Path
)
if
(
Mask
==
FALSE
&
(
!
HDR
$
interleave
==
"
bil
"
)
&
(
!
HDR
$
interleave
==
"
BIL
"
))
{
message
(
""
)
message
(
"
*********************************************************
"
)
message
(
"
The image format may not compatible with the processing chain
"
)
message
(
"
Image format expected:
"
)
message
(
"
ENVI hdr file with band interleaved by line (BIL) file format
"
)
message
(
""
)
message
(
"
Current Image format
"
)
print
(
HDR
$
interleave
)
message
(
'
Please run the function named
'
)
print
(
'
Convert.Raster2BIL
'
)
message
(
'
in order to convert your raster data
'
)
message
(
'
or use appropriate software
'
)
message
(
'
*********************************************************
'
)
message
(
''
)
message
(
"
Please run the function named
"
)
print
(
"
Convert.Raster2BIL
"
)
message
(
"
in order to convert your raster data
"
)
message
(
"
or use appropriate software
"
)
message
(
"
*********************************************************
"
)
message
(
""
)
stop
()
}
else
if
(
Mask
==
FALSE
&
((
HDR
$
interleave
==
'
bil
'
)
|
(
HDR
$
interleave
==
'
BIL
'
))){
if
(
HDR
$
`wavelength units`
==
'
Unknown
'
)
{
message
(
''
)
message
(
'
*********************************************************
'
)
message
(
'
Please make sure the wavelengths are in nanometers
'
)
message
(
'
if not, stop processing and convert wavelengths in nanometers in HDR file
'
)
message
(
'
*********************************************************
'
)
message
(
''
)
}
else
if
(
Mask
==
FALSE
&
((
HDR
$
interleave
==
"
bil
"
)
|
(
HDR
$
interleave
==
"
BIL
"
)))
{
if
(
HDR
$
`wavelength units`
==
"
Unknown
"
)
{
message
(
""
)
message
(
"
*********************************************************
"
)
message
(
"
Please make sure the wavelengths are in nanometers
"
)
message
(
"
if not, stop processing and convert wavelengths in nanometers in HDR file
"
)
message
(
"
*********************************************************
"
)
message
(
""
)
}
if
((
!
HDR
$
`wavelength units`
==
'
Nanometers
'
)
&
(
!
HDR
$
`wavelength units`
==
'
nanometers
'
)){
message
(
''
)
message
(
'
*********************************************************
'
)
message
(
'
Please make sure the wavelengths are in nanometers
'
)
message
(
'
if not, stop processing and convert wavelengths in nanometers in HDR file
'
)
message
(
'
*********************************************************
'
)
message
(
''
)
if
((
!
HDR
$
`wavelength units`
==
"
Nanometers
"
)
&
(
!
HDR
$
`wavelength units`
==
"
nanometers
"
))
{
message
(
""
)
message
(
"
*********************************************************
"
)
message
(
"
Please make sure the wavelengths are in nanometers
"
)
message
(
"
if not, stop processing and convert wavelengths in nanometers in HDR file
"
)
message
(
"
*********************************************************
"
)
message
(
""
)
}
if
(
HDR
$
`wavelength units`
==
'
micrometers
'
)
{
message
(
''
)
message
(
'
*********************************************************
'
)
message
(
'
Please convert wavelengths in nanometers in HDR file
'
)
message
(
'
*********************************************************
'
)
message
(
''
)
if
(
HDR
$
`wavelength units`
==
"
micrometers
"
)
{
message
(
""
)
message
(
"
*********************************************************
"
)
message
(
"
Please convert wavelengths in nanometers in HDR file
"
)
message
(
"
*********************************************************
"
)
message
(
""
)
stop
()
}
if
((
HDR
$
`wavelength units`
==
'
nanometers
'
)
|
(
HDR
$
`wavelength units`
==
'
Nanometers
'
)){
message
(
''
)
message
(
'
*********************************************************
'
)
message
(
'
All information seem OK for image processing
'
)
message
(
'
*********************************************************
'
)
message
(
''
)
if
((
HDR
$
`wavelength units`
==
"
nanometers
"
)
|
(
HDR
$
`wavelength units`
==
"
Nanometers
"
))
{
message
(
""
)
message
(
"
*********************************************************
"
)
message
(
"
All information seem OK for image processing
"
)
message
(
"
*********************************************************
"
)
message
(
""
)
}
}
}
else
{
message
(
''
)
message
(
'
*********************************************************
'
)
message
(
'
The following HDR file was expected, but could not be found:
'
)
message
(
""
)
message
(
"
*********************************************************
"
)
message
(
"
The following HDR file was expected, but could not be found:
"
)
print
(
HDR.Path
)
message
(
'
The image format may not compatible with the processing chain
'
)
message
(
'
Image format expected:
'
)
message
(
'
ENVI hdr file with band interleaved by line (BIL) file format
'
)
message
(
''
)
message
(
'
Please run the function named
'
)
print
(
'
Convert.Raster2BIL
'
)
message
(
'
in order to convert your raster data
'
)
message
(
'
*********************************************************
'
)
message
(
''
)
message
(
"
The image format may not compatible with the processing chain
"
)
message
(
"
Image format expected:
"
)
message
(
"
ENVI hdr file with band interleaved by line (BIL) file format
"
)
message
(
""
)
message
(
"
Please run the function named
"
)
print
(
"
Convert.Raster2BIL
"
)
message
(
"
in order to convert your raster data
"
)
message
(
"
*********************************************************
"
)
message
(
""
)
stop
()
}
return
()
return
(
invisible
()
)
}
R/Lib_ContinuumRemoval.R
View file @
c016a605
...
...
@@ -9,187 +9,189 @@
# This Library is dedicated to the computation of the continuum removal
# ===============================================================================
#' prepares data to run multithreaded continuum removal
#'
#' @param Spectral.Data initial data matrix (nb samples x nb bands)
#' @param nbCPU
#' @param Spectral information about spectral bands
#'
#' @return samples from image and updated number of pixels to sampel if necessary
Apply.Continuum.Removal
=
function
(
Spectral.Data
,
Spectral
,
nbCPU
=
1
){
if
(
!
length
(
Spectral
$
WaterVapor
)
==
0
){
Spectral.Data
=
Spectral.Data
[,
-
Spectral
$
WaterVapor
]
# prepares data to run multithreaded continuum removal
#
# @param Spectral.Data initial data matrix (nb samples x nb bands)
# @param nbCPU
# @param Spectral information about spectral bands
#
# @return samples from image and updated number of pixels to sampel if necessary
#' @importFrom snow splitRows
#' @import future
#' @importFrom future.apply future_lapply
Apply.Continuum.Removal
<-
function
(
Spectral.Data
,
Spectral
,
nbCPU
=
1
)
{
if
(
!
length
(
Spectral
$
WaterVapor
)
==
0
)
{
Spectral.Data
<-
Spectral.Data
[,
-
Spectral
$
WaterVapor
]
}
# split data to perform continuum removal on into reasonable amount of data
nb.Values
=
size
(
Spectral.Data
)[
1
]
*
size
(
Spectral.Data
)[
2
]
if
(
nb.Values
>
0
){
nb.Values
<-
dim
(
Spectral.Data
)[
1
]
*
dim
(
Spectral.Data
)[
2
]
if
(
nb.Values
>
0
)
{
# corresponds to ~ 40 Mb data, but CR tends to requires ~ 10 times memory
# avoids memory crash
Max.nb.Values
=
2e6
nb.CR
=
ceiling
(
nb.Values
/
Max.nb.Values
)
Spectral.Data
=
splitRows
(
Spectral.Data
,
nb.CR
)
Max.nb.Values
<-
2e6
nb.CR
<-
ceiling
(
nb.Values
/
Max.nb.Values
)
Spectral.Data
<-
splitRows
(
Spectral.Data
,
nb.CR
)
# perform multithread continuum removal
plan
(
multiprocess
,
workers
=
nbCPU
)
## Parallelize using four cores
Schedule.Per.Thread
=
ceiling
(
nb.CR
/
nbCPU
)
Spectral.Data.tmp
=
future_lapply
(
Spectral.Data
,
FUN
=
ContinuumRemoval
,
Spectral.Bands
=
Spectral
$
Wavelength
,
future.scheduling
=
Schedule.Per.Thread
)
Schedule.Per.Thread
<-
ceiling
(
nb.CR
/
nbCPU
)
Spectral.Data.tmp
<-
future_lapply
(
Spectral.Data
,
FUN
=
ContinuumRemoval
,
Spectral.Bands
=
Spectral
$
Wavelength
,
future.scheduling
=
Schedule.Per.Thread
)
plan
(
sequential
)
Spectral.Data
=
do.call
(
"rbind"
,
Spectral.Data.tmp
)
Spectral.Data
<-
do.call
(
"rbind"
,
Spectral.Data.tmp
)
rm
(
Spectral.Data.tmp
)
}
else
{
# edit 31-jan-2018
# resize to delete first and last band as in continuum removal
Spectral.Data
=
Spectral.Data
[,
-
c
(
1
,
2
)]
Spectral.Data
<-
Spectral.Data
[,
-
c
(
1
,
2
)]
}
gc
()
return
(
Spectral.Data
)
}
#
'
Computes continuum removal for matrix shaped data: more efficient than
#
'
processing individual spectra
#
'
the convex hull is based on the computation of the derivative between R at a
#
'
given spectral band and R at the following bands
#
'
#
'
@param Minit initial data matrix (nb samples x nb bands)
#
'
@param Spectral.Bands information about spectral bands
#
'
#
'
@return samples from image and updated number of pixels to sampel if necessary
ContinuumRemoval
=
function
(
Minit
,
Spectral.Bands
){
# Computes continuum removal for matrix shaped data: more efficient than
# processing individual spectra
# the convex hull is based on the computation of the derivative between R at a
# given spectral band and R at the following bands
#
# @param Minit initial data matrix (nb samples x nb bands)
# @param Spectral.Bands information about spectral bands
#
# @return samples from image and updated number of pixels to sampel if necessary
ContinuumRemoval
<-
function
(
Minit
,
Spectral.Bands
)
{
# Filter and prepare data prior to continuum removal
CR.data
=
Filter.Prior.CR
(
Minit
,
Spectral.Bands
)
Minit
=
CR.data
$
Minit
nb.Bands
=
size
(
Minit
)[
2
]
CR.data
$
Minit
=
c
()
Spectral.Bands
=
CR.data
$
Spectral.Bands
nb.Samples
=
CR.data
$
nb.Samples
CR.data
<-
Filter.Prior.CR
(
Minit
,
Spectral.Bands
)
Minit
<-
CR.data
$
Minit
nb.Bands
<-
dim
(
Minit
)[
2
]
CR.data
$
Minit
<-
c
()
Spectral.Bands
<-
CR.data
$
Spectral.Bands
nb.Samples
<-
CR.data
$
nb.Samples
# if samples to be considered
if
(
nb.Samples
>
0
){
if
(
nb.Samples
>
0
)
{
# initialization:
# spectral band corresponding to each element of the data matrix
Lambda
=
repmat
(
matrix
(
Spectral.Bands
,
nrow
=
1
),
nb.Samples
,
1
)
Lambda
<-
repmat
(
matrix
(
Spectral.Bands
,
nrow
=
1
),
nb.Samples
,
1
)
# prepare matrices used to check evolution of the CR process
# - elements still not processed through continuum removal: initialization to 1
Still.Need.CR
=
matrix
(
1
,
nrow
=
nb.Samples
,
ncol
=
nb.Bands
)
Still.Need.CR
<-
matrix
(
1
,
nrow
=
nb.Samples
,
ncol
=
nb.Bands
)
# - value of the convex hull: initially set to 0
Convex.Hull
=
matrix
(
0
,
nrow
=
nb.Samples
,
ncol
=
nb.Bands
)
Convex.Hull
<-
matrix
(
0
,
nrow
=
nb.Samples
,
ncol
=
nb.Bands
)
# - reflectance value for latest interception with convex hull:
# initialization to value of the first reflectance measurement
Intercept.Hull
=
repmat
(
matrix
(
Minit
[,
1
],
ncol
=
1
),
1
,
nb.Bands
)
Intercept.Hull
<-
repmat
(
matrix
(
Minit
[,
1
],
ncol
=
1
),
1
,
nb.Bands
)
# - spectral band of latest interception
Latest.Intercept
=
repmat
(
matrix
(
Spectral.Bands
[
1
],
ncol
=
1
),
nb.Samples
,
nb.Bands
)
Latest.Intercept
<-
repmat
(
matrix
(
Spectral.Bands
[
1
],
ncol
=
1
),
nb.Samples
,
nb.Bands
)
# number of spectral bands found as intercept
nb.Intercept
=
0
nb.Intercept
<-
0
# continues until arbitrary stopping criterion:
# stops when reach last spectral band (all values before last = 0)
while
(
max
(
Still.Need.CR
[,
1
:
(
nb.Bands
-
2
)])
==
1
&
(
nb.Intercept
<=
(
nb.Bands
/
2
))){
nb.Intercept
=
nb.Intercept
+
1
while
(
max
(
Still.Need.CR
[,
1
:
(
nb.Bands
-
2
)])
==
1
&
(
nb.Intercept
<=
(
nb.Bands
/
2
)))
{
nb.Intercept
<-
nb.Intercept
+
1
# Mstep give the position of the values to be updated
Update.Data
=
matrix
(
1
,
nrow
=
nb.Samples
,
ncol
=
nb.Bands
)
Update.Data
[,
nb.Bands
]
=
0
Update.Data
<-
matrix
(
1
,
nrow
=
nb.Samples
,
ncol
=
nb.Bands
)
Update.Data
[,
nb.Bands
]
<-
0
# initial step: first column set to 0; following steps: all bands below
# max of the convex hull are set to 0
Update.Data
[
which
((
Lambda
-
Latest.Intercept
)
<
0
)]
=
0
Update.Data
[
which
((
Lambda
-
Latest.Intercept
)
<
0
)]
<-
0
# compute slope for each coordinate
Slope
=
(
Minit
-
Intercept.Hull
)
/
(
Lambda
-
Latest.Intercept
)
*
Still.Need.CR
Slope
<-
(
Minit
-
Intercept.Hull
)
/
(
Lambda
-
Latest.Intercept
)
*
Still.Need.CR
# set current spectral band and previous bands to -9999
if
(
!
length
(
which
(
Still.Need.CR
==
0
))
==
0
){
Slope
[
which
(
Still.Need.CR
==
0
)]
=
-9999
if
(
!
length
(
which
(
Still.Need.CR
==
0
))
==
0
)
{
Slope
[
which
(
Still.Need.CR
==
0
)]
<-
-9999
}
if
(
!
length
(
which
(
is.na
(
Slope
)))
==
0
){
Slope
[
which
(
is.na
(
Slope
))]
=
-9999
if
(
!
length
(
which
(
is.na
(
Slope
)))
==
0
)
{
Slope
[
which
(
is.na
(
Slope
))]
<-
-9999
}
# get max index for each row and convert into linear index
Index.Max.Slope
=
RowToLinear
(
max.col
(
Slope
,
ties.method
=
"last"
),
nb.Samples
,
nb.Bands
)
Index.Max.Slope
<-
RowToLinear
(
max.col
(
Slope
,
ties.method
=
"last"
),
nb.Samples
,
nb.Bands
)
# !!!! OPTIM: replace repmat with column operation
# update coordinates of latest intercept
Latest.Intercept
=
repmat
(
matrix
(
Lambda
[
Index.Max.Slope
],
ncol
=
1
),
1
,
nb.Bands
)
Latest.Intercept
<-
repmat
(
matrix
(
Lambda
[
Index.Max.Slope
],
ncol
=
1
),
1
,
nb.Bands
)
# update latest intercept
Intercept.Hull
=
repmat
(
matrix
(
Minit
[
Index.Max.Slope
],
ncol
=
1
),
1
,
nb.Bands
)
Intercept.Hull
<-
repmat
(
matrix
(
Minit
[
Index.Max.Slope
],
ncol
=
1
),
1
,
nb.Bands
)
# values corresponding to the domain between the two continuum maxima
Update.Data
[
which
((
Lambda
-
Latest.Intercept
)
>=
0
|
Latest.Intercept
==
Spectral.Bands
[
nb.Bands
])]
=
0
Update.Data
[
which
((
Lambda
-
Latest.Intercept
)
>=
0
|
Latest.Intercept
==
Spectral.Bands
[
nb.Bands
])]
<-
0
# values to eliminate for the next analysis: all spectral bands before latest intercept
Still.Need.CR
[
which
((
Lambda
-
Latest.Intercept
)
<
0
)]
=
0
Still.Need.CR
[
which
((
Lambda
-
Latest.Intercept
)
<
0
)]
<-
0
# the max slope is known, as well as the coordinates of the beginning and ending
# a matrix now has to be built