Add SpectralIndex transform (#300)

* Add SpectralIndex transform

* Use _assert utility

* Add tests for auxiliary functions

* Update src/transforms/spectralindex.jl

Co-authored-by: Elias Carvalho <[email protected]>

---------

Co-authored-by: Elias Carvalho <[email protected]>

Author: juliohm

Project.toml

@@ -16,6 +16,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 NelderMead = "2f6b4ddb-b4ff-44c0-b59b-2ab99302f970"
 PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+SpectralIndices = "df0093a1-273d-40bc-819a-796ec3476907"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 TableDistances = "e5d66e97-8c70-46bb-8b66-04a2d73ad782"
@@ -36,6 +37,7 @@ LinearAlgebra = "1.9"
 NelderMead = "0.4"
 PrettyTables = "2"
 Random = "1.9"
+SpectralIndices = "0.2"
 Statistics = "1.9"
 StatsBase = "0.33, 0.34"
 TableDistances = "1.0"

docs/src/transforms.md

@@ -242,6 +242,12 @@ SDS
 ProjectionPursuit
 ```
 
+## SpectralIndex
+
+```@docs
+SpectralIndex
+```
+
 ## KMedoids
 
 ```@docs

src/TableTransforms.jl

@@ -25,6 +25,7 @@ using Distributions: ContinuousUnivariateDistribution, Normal
 using InverseFunctions: NoInverse, inverse as invfun
 using StatsBase: AbstractWeights, Weights, sample
 using Distributed: CachingPool, pmap, workers
+using SpectralIndices: indices, bands, compute
 using NelderMead: optimise
 
 import Distributions: quantile, cdf
@@ -91,6 +92,7 @@ export
   DRS,
   SDS,
   ProjectionPursuit,
+  SpectralIndex,
   KMedoids,
   Closure,
   Remainder,
@@ -101,6 +103,10 @@ export
   RowTable,
   ColTable,
   →,
-  ⊔
+  ⊔,
+
+  # utilities
+  spectralindices,
+  spectralbands
 
 end

src/transforms.jl

@@ -286,6 +286,7 @@ include("transforms/quantile.jl")
 include("transforms/functional.jl")
 include("transforms/eigenanalysis.jl")
 include("transforms/projectionpursuit.jl")
+include("transforms/spectralindex.jl")
 include("transforms/kmedoids.jl")
 include("transforms/closure.jl")
 include("transforms/remainder.jl")

src/transforms/spectralindex.jl

@@ -0,0 +1,110 @@
+# ------------------------------------------------------------------
+# Licensed under the MIT License. See LICENSE in the project root.
+# ------------------------------------------------------------------
+
+"""
+    SpectralIndex(name, [bands...])
+
+Compute the spectral index of given `name` from SpectralIndices.jl.
+Optionally, specify the column names corresponding to spectral `bands`
+as keyword arguments.
+
+# Examples
+
+```julia
+# vegetation index
+SpectralIndex("NDVI")
+
+# water index
+SpectralIndex("NDWI")
+
+# specify "R" (red) and "N" (near infra red) columns
+SpectralIndex("NDVI", R="col1", N="col4")
+```
+
+### Notes
+
+The list of supported indices can be obtained with `spectralindices()`
+and the list of spectral bands for a given index can be obtained with
+`spectralbands(name)`.
+"""
+struct SpectralIndex{B} <: StatelessFeatureTransform
+  name::String
+  bands::B
+
+  function SpectralIndex(name, bands)
+    sname = string(name)
+    _assert(sname ∈ keys(indices), "$sname not found in SpectralIndices.jl")
+    sbands = if isempty(bands)
+      nothing
+    else
+      skeys = string.(keys(bands))
+      vkeys = Tuple(indices[sname].bands)
+      _assert(skeys ⊆ vkeys, "bands $skeys are not valid for spectral index $sname, please choose from $vkeys")
+      svals = string.(values(values(bands)))
+      (; zip(Symbol.(skeys), svals)...)
+    end
+    new{typeof(sbands)}(sname, sbands)
+  end
+end
+
+SpectralIndex(name; bands...) = SpectralIndex(name, bands)
+
+function applyfeat(transform::SpectralIndex, feat, prep)
+  # retrieve spectral index
+  iname = transform.name
+  index = indices[iname]
+
+  # extract band names from feature table
+  cols = Tables.columns(feat)
+  names = Tables.columnnames(cols)
+  bnames = Symbol.(index.bands)
+  tbands = transform.bands
+  snames = map(bnames) do b
+    if !isnothing(tbands) && b ∈ keys(tbands)
+      Symbol(tbands[b])
+    else
+      b
+    end
+  end
+
+  # throw helpful error message in case of invalid names
+  if !(snames ⊆ names)
+    notfound = setdiff(snames, names)
+    required = ["$(b.short_name): $(b.long_name)" for b in spectralbands(iname)]
+    pprint(names) = "\"" * join(string.(names), "\", \"", "\" and \"") * "\""
+    throw(ArgumentError("""columns $(pprint(notfound)) not found in table.
+
+      Please specify valid columns names for the spectral bands as keyword arguments.
+
+      Required bands for $iname:
+
+        $(join(required, "\n    "))
+
+      Available column names: $(pprint(names))
+    """))
+  end
+
+  # compute index for all locations
+  icols = [b => Tables.getcolumn(cols, n) for (b, n) in zip(bnames, snames)]
+  ivals = compute(index; icols...)
+
+  # new table with index feature
+  newfeat = (; Symbol(iname) => ivals) |> Tables.materializer(feat)
+
+  newfeat, nothing
+end
+
+"""
+    spectralindices()
+
+List of spectral indices supported by SpectralIndices.jl.
+"""
+spectralindices() = indices
+
+"""
+    spectralbands(name)
+
+List of spectral bands for spectral index of given `name`.
+"""
+spectralbands(name) = [bands[b] for b in indices[name].bands]

test/transforms.jl

@@ -31,6 +31,7 @@ transformfiles = [
   "functional.jl",
   "eigenanalysis.jl",
   "projectionpursuit.jl",
+  "spectralindex.jl",
   "kmedoids.jl",
   "closure.jl",
   "remainder.jl",

test/transforms/spectralindex.jl

@@ -0,0 +1,41 @@
+@testset "SpectralIndex" begin
+  @test !isrevertible(SpectralIndex("NDVI"))
+
+  # standard names
+  t = Table(R=[1,2,3], N=[4,5,6])
+  T = SpectralIndex("NDVI")
+  @test T(t).NDVI ≈ [0.6, 0.42857142857142855, 0.3333333333333333]
+  T = SpectralIndex("NDVI", R="R")
+  @test T(t).NDVI ≈ [0.6, 0.42857142857142855, 0.3333333333333333]
+  T = SpectralIndex("NDVI", N="N")
+  @test T(t).NDVI ≈ [0.6, 0.42857142857142855, 0.3333333333333333]
+  T = SpectralIndex("NDVI", R="R", N="N")
+  @test T(t).NDVI ≈ [0.6, 0.42857142857142855, 0.3333333333333333]
+
+  # some non-standard names
+  t = Table(R=[1,2,3], NIR=[4,5,6])
+  T = SpectralIndex("NDVI", N="NIR")
+  @test T(t).NDVI ≈ [0.6, 0.42857142857142855, 0.3333333333333333]
+  T = SpectralIndex("NDVI")
+  @test_throws ArgumentError T(t)
+  T = SpectralIndex("NDVI", R="RED")
+  @test_throws ArgumentError T(t)
+
+  # all non-standard names
+  t = Table(RED=[1,2,3], NIR=[4,5,6])
+  T = SpectralIndex("NDVI", R="RED", N="NIR")
+  @test T(t).NDVI ≈ [0.6, 0.42857142857142855, 0.3333333333333333]
+  T = SpectralIndex("NDVI", R="RED")
+  @test_throws ArgumentError T(t)
+  T = SpectralIndex("NDVI", N="NIR")
+  @test_throws ArgumentError T(t)
+
+  # bands as symbols
+  t = Table(RED=[1,2,3], NIR=[4,5,6])
+  T = SpectralIndex("NDVI", R=:RED, N=:NIR)
+  @test T(t).NDVI ≈ [0.6, 0.42857142857142855, 0.3333333333333333]
+
+  # auxiliary functions
+  @test spectralindices() isa Dict
+  @test spectralbands("NDVI") isa Vector
+end