diff --git a/Project.toml b/Project.toml index 6b59126f..31706f98 100644 --- a/Project.toml +++ b/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" diff --git a/docs/src/transforms.md b/docs/src/transforms.md index b0c9bfbc..f5d3dae7 100644 --- a/docs/src/transforms.md +++ b/docs/src/transforms.md @@ -242,6 +242,12 @@ SDS ProjectionPursuit ``` +## SpectralIndex + +```@docs +SpectralIndex +``` + ## KMedoids ```@docs diff --git a/src/TableTransforms.jl b/src/TableTransforms.jl index e624c1d3..101e5b31 100644 --- a/src/TableTransforms.jl +++ b/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 diff --git a/src/transforms.jl b/src/transforms.jl index e0c6082a..adbebdb4 100644 --- a/src/transforms.jl +++ b/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") diff --git a/src/transforms/spectralindex.jl b/src/transforms/spectralindex.jl new file mode 100644 index 00000000..46863fce --- /dev/null +++ b/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] diff --git a/test/transforms.jl b/test/transforms.jl index 77461353..da052c5c 100644 --- a/test/transforms.jl +++ b/test/transforms.jl @@ -31,6 +31,7 @@ transformfiles = [ "functional.jl", "eigenanalysis.jl", "projectionpursuit.jl", + "spectralindex.jl", "kmedoids.jl", "closure.jl", "remainder.jl", diff --git a/test/transforms/spectralindex.jl b/test/transforms/spectralindex.jl new file mode 100644 index 00000000..63e2e9ca --- /dev/null +++ b/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