Merge pull request #2 from gaelforget/gfdev01

add classifiers & recipes

.travis.yml

@@ -2,7 +2,6 @@
 language: julia
 os:
   - linux
-  - osx
 julia:
   - 1.2
   - 1.3
@@ -23,3 +22,5 @@ jobs:
                                     Pkg.develop(PackageSpec(path=pwd()))'
         - julia --project=docs/ docs/make.jl
       after_success: skip
+  allow_failures:
+    - julia: nightly

Project.toml

@@ -1,7 +1,13 @@
 name = "OceanColorData"
 uuid = "39357346-8f20-4302-be7b-c20dcd116b7a"
 authors = ["gaelforget <gforget@mit.edu>"]
-version = "0.1.0"
+version = "0.1.1"
+
+[deps]
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
+NetCDF = "30363a11-5582-574a-97bb-aa9a979735b9"
 
 [compat]
+Distributions = "0.23"
+NetCDF = "0.10"
 julia = "1"

src/Classifiers.jl

@@ -0,0 +1,121 @@
+
+using Distributions, NetCDF
+
+"""
+    FuzzyClassification(M,Sinv,Rrs)
+
+Apply fuzzy membership classifier (`M` + `Sinv`) to a vector of remotely sensed
+reflectcances (`Rrs`) and returns a membership vector (values between 0 and 1).
+"""
+function FuzzyClassification(M,Sinv,Rrs)
+        f=Array{Any,1}(undef,length(M))
+        for ii=1:length(M)
+                X=vec(Rrs)-M[ii]
+                Z=transpose(X)*Sinv[ii]*X
+                f[ii]=ccdf(Chisq(6),Z)
+        end
+
+        return f
+end
+
+"""
+    Jackson2017()
+
+`Fuzzy logic` classifiers defined in [Moore et al 2009](https://doi.org/10.1016/j.rse.2009.07.016) and
+[Jackson et al 2017](http://dx.doi.org/10.1016/j.rse.2017.03.036) can be used to
+assign optical class memberships from an `Rrs` vector. This function provides
+vector `M` and inverse matrix `Sinv` for the Jackson et al 2017 classifier.
+
+Credits: T. Jackson kindly provided the `J17.nc` classifier file
+"""
+function Jackson2017()
+        fil=joinpath(dirname(pathof(OceanColorData)),"../examples/J17.nc")
+        tmpM = ncread(fil, "cluster_means")
+        tmpSinv = ncread(fil, "inverse_covariance")
+
+        M=Array{Any,1}(undef,14)
+        Sinv=Array{Any,1}(undef,14)
+        for ii=1:length(M)
+                M[ii]=vec(tmpM[ii,:])
+                Sinv[ii]=tmpSinv[1:6,1:6,ii]
+        end
+
+        return M,Sinv
+end
+
+"""
+    Moore2009()
+
+`Fuzzy logic` classifiers defined in [Moore et al 2009](https://doi.org/10.1016/j.rse.2009.07.016) and
+[Jackson et al 2017](http://dx.doi.org/10.1016/j.rse.2017.03.036) can be used to
+assign optical class memberships from an `Rrs` vector. This function provides
+vector `M` and inverse matrix `Sinv` for the Moore et al 2009 classifier.
+"""
+function Moore2009()
+
+        #Moore et al 2009 mean vectors:
+        M=Array{Any,1}(undef,8)
+        M[1]=vec([0.0234 0.0192 0.0129 0.0075 0.0031 0.0002])
+        M[2]=vec([0.0162 0.0141 0.0112 0.0073 0.0034 0.0002])
+        M[3]=vec([0.0107 0.0098 0.0092 0.0070 0.0039 0.0003])
+        M[4]=vec([0.0065 0.0064 0.0070 0.0064 0.0048 0.0006])
+        M[5]=vec([0.0033 0.0034 0.0042 0.0042 0.0043 0.0009])
+        M[6]=vec([0.0064 0.0074 0.0105 0.0116 0.0140 0.0041])
+        M[7]=vec([0.0121 0.0140 0.0192 0.0204 0.0231 0.0084])
+        M[8]=vec([0.0184 0.0230 0.0333 0.0359 0.0409 0.0137])
+
+        #Moore et al 2009 covariance matrices:
+        S=Array{Any,1}(undef,8)
+        S[1]=[0.00000959 0.00000556 0.00000138 -0.00000034 -0.00000024 -0.00000003;
+        0.00000556 0.00000493 0.00000193 0.00000060 0.00000023 0.00000001;
+        0.00000138 0.00000193 0.00000282 0.00000223 0.00000119 0.00000007;
+        -0.00000034 0.00000060 0.00000223 0.00000232 0.00000119 0.00000007;
+        -0.00000024 0.00000023 0.00000119 0.00000119 0.00000071 0.00000005;
+        -0.00000003 0.00000001 0.00000007 0.00000007 0.00000005 0.00000001]
+        S[2]=[0.00000346 0.00000186 -0.00000011 -0.00000060 -0.00000062 -0.00000005;
+        0.00000186 0.00000228 0.00000086 0.00000033 -0.00000007 0.00000003;
+        -0.00000011 0.00000086 0.00000231 0.00000221 0.00000145 0.00000023;
+        -0.00000060 0.00000033 0.00000221 0.00000266 0.00000191 0.00000034;
+        -0.00000062 -0.00000007 0.00000145 0.00000191 0.00000175 0.00000041;
+        -0.00000005 0.00000003 0.00000023 0.00000034 0.00000041 0.00000021]
+        S[3]=[0.00000241 0.00000144 0.00000035 -0.00000031 -0.00000063 -0.00000006;
+        0.00000144 0.00000138 0.00000076 0.00000015 -0.00000021 -0.00000001;
+        0.00000035 0.00000076 0.00000161 0.00000156 0.00000121 0.00000016;
+        -0.00000031 0.00000015 0.00000156 0.00000227 0.00000209 0.00000031;
+        -0.00000063 -0.00000021 0.00000121 0.00000209 0.00000225 0.00000037;
+        -0.00000006 -0.00000001 0.00000016 0.00000031 0.00000037 0.00000013]
+        S[4]=[0.00000166 0.00000091 0.00000034 -0.00000009 -0.00000080 -0.00000025;
+        0.00000091 0.00000097 0.00000071 0.00000025 -0.00000041 -0.00000015;
+        0.00000034 0.00000071 0.00000118 0.00000103 0.00000072 0.00000003;
+        -0.00000009 0.00000025 0.00000103 0.00000137 0.00000162 0.00000025;
+        -0.00000080 -0.00000041 0.00000072 0.00000162 0.00000290 0.00000065;
+        -0.00000025 -0.00000015 0.00000003 0.00000025 0.00000065 0.00000050]
+        S[5]=[0.00000178 0.00000132 0.00000104 0.00000081 0.00000018 -0.00000014;
+        0.00000132 0.00000127 0.00000121 0.00000099 0.00000034 -0.00000010;
+        0.00000104 0.00000121 0.00000150 0.00000142 0.00000110 0.00000013;
+        0.00000081 0.00000099 0.00000142 0.00000158 0.00000177 0.00000042;
+        0.00000018 0.00000034 0.00000110 0.00000177 0.00000351 0.00000131;
+        -0.00000014 -0.00000010 0.00000013 0.00000042 0.00000131 0.00000081]
+        S[6]=[0.00000715 0.00000586 0.00000409 0.00000292 0.00000005 -0.00000075;
+        0.00000586 0.00000589 0.00000520 0.00000398 0.00000027 -0.00000114;
+        0.00000409 0.00000520 0.00000634 0.00000541 0.00000188 -0.00000097;
+        0.00000292 0.00000398 0.00000541 0.00000528 0.00000392 0.00000070;
+        0.00000005 0.00000027 0.00000188 0.00000392 0.00000995 0.00000657;
+        -0.00000075 -0.00000114 -0.00000097 0.00000070 0.00000657 0.00000819]
+        S[7]=[0.00002625 0.00001981 0.00001058 0.00000544 -0.00000654 -0.00001122;
+        0.00001981 0.00001745 0.00001314 0.00000822 -0.00000431 -0.00001228;
+        0.00001058 0.00001314 0.00001629 0.00001226 0.00000035 -0.00001311;
+        0.00000544 0.00000822 0.00001226 0.00001170 0.00000742 -0.00000500;
+        -0.00000654 -0.00000431 0.00000035 0.00000742 0.00002241 0.00001782;
+        -0.00001122 -0.00001228 -0.00001311 -0.00000500 0.00001782 0.00003987]
+        S[8]=[0.00001186 0.00001134 0.00001139 0.00000919 0.00000395 -0.00000186;
+        0.00001134 0.00001484 0.00002034 0.00001907 0.00001531 0.00000087;
+        0.00001139 0.00002034 0.00003467 0.00003546 0.00003555 0.00000708;
+        0.00000919 0.00001907 0.00003546 0.00003907 0.00004604 0.00001733;
+        0.00000395 0.00001531 0.00003555 0.00004604 0.00007306 0.00004953;
+        -0.00000186 0.00000087 0.00000708 0.00001733 0.00004953 0.00006542];
+
+        Sinv=inv.(S)
+
+        return M,Sinv
+end

src/Conversions.jl

@@ -0,0 +1,82 @@
+
+"""
+    RemotelySensedReflectance(rirr_in,wvbd_in,wvbd_out)
+
+Compute remotely sensed reflectance at wvbd_out from
+irradiance reflectance at wvbd_in.
+
+```
+wvbd_out=Float64.([412, 443, 490, 510, 555, 670])
+wvbd_in=Float64.([400,425,450,475,500,525,550,575,600,625,650,675,700])
+rirr_in=1e-3*[23.7641,26.5037,27.9743,30.4914,28.1356,
+    21.9385,18.6545,13.5100,5.6338,3.9272,2.9621,2.1865,1.8015]
+
+rrs_out=RemotelySensedReflectance(rirr_in,wvbd_in,wvbd_out)
+
+using Plots
+plot(vec(rrs_out),linewidth=4,lab="recomputed RRS")
+rrs_ref=1e-3*[4.4099, 4.8533, 5.1247, 4.5137, 3.0864, 0.4064]
+plot!(rrs_ref,linewidth=4,ls=:dash,lab="reference result")
+```
+"""
+function RemotelySensedReflectance(rirr_in,wvbd_in,wvbd_out)
+    rrs_out=similar(wvbd_out)
+
+    nin=length(wvbd_in)
+    nout=length(wvbd_out)
+
+    jj=Array{Int64,1}(undef,nin)
+    ww=Array{Float64,1}(undef,nin)
+    for ii=1:6
+        tmp=wvbd_out[ii].-wvbd_in
+        kk=maximum(findall(tmp.>=0))
+        jj[ii]=kk
+        ww[ii]=tmp[kk]/(wvbd_in[kk+1]-wvbd_in[kk])
+    end
+
+    tmp=rirr_in/3
+    rrs_tmp=(0.52*tmp)./(1.0 .-1.7*tmp);
+
+    rrs_out=Array{Float64,1}(undef,nout)
+    for vv=1:6
+        tmp0=rrs_tmp[jj[vv]]
+        tmp1=rrs_tmp[jj[vv]+1]
+        rrs_out[vv]=tmp0.*(1-ww[vv])+tmp1.*ww[vv]
+    end
+
+    return rrs_out
+end
+
+
+"""
+    RrsToChla(Rrs; Eqn="OC4")
+
+Satellite `Chl_a` estimates typicaly derive from remotely sensed reflectances
+using the blue/green reflectance ratio method and a polynomial formulation.
+
+```
+wvbd_out=Float64.([412, 443, 490, 510, 555, 670])
+wvbd_in=Float64.([400,425,450,475,500,525,550,575,600,625,650,675,700])
+rirr_in=1e-3*[23.7641,26.5037,27.9743,30.4914,28.1356,
+    21.9385,18.6545,13.5100,5.6338,3.9272,2.9621,2.1865,1.8015]
+
+rrs_out=RemotelySensedReflectance(rirr_in,wvbd_in,wvbd_out)
+chla_out=RrsToChla(rrs_out)
+```
+"""
+function RrsToChla(Rrs;Eqn="OC4")
+    RRSB=max.(Rrs[2],Rrs[3]) #blue
+    RRSG=Rrs[5] #green
+    X = log10.(RRSB./RRSG) #ratio of blue to green
+
+    if Eqn=="OC4"
+        C=[0.3272, -2.9940, +2.7218, -1.2259, -0.5683] #OC4 algorithms (SeaWifs, CCI)
+    elseif Eqn=="OC3M-547"
+        C=[0.2424, -2.7423, 1.8017, 0.0015, -1.2280] #OC3M-547 algorithm (MODIS)
+    else
+        error("this case has not been implemented yet...")
+    end
+
+    a0=C[1]; a1=C[2]; a2=C[3]; a3=C[4]; a4=C[5];
+    chld=exp10.(a0.+a1*X+a2*X.^2+a3*X.^3+a4*X.^4); #apply polynomial recipe
+end

src/OceanColorData.jl

@@ -1,54 +1,9 @@
 module OceanColorData
 
-export RemotelySensedReflectance
+include("Conversions.jl")
+include("Classifiers.jl")
 
-"""
-    RemotelySensedReflectance(rirr_in,wvbd_in,wvbd_out)
-
-Compute remotely sensed reflectance at wvbd_out from
-irradiance reflectance at wvbd_in.
-
-```
-wvbd_out=Float64.([412, 443, 490, 510, 555, 670])
-wvbd_in=Float64.([400,425,450,475,500,525,550,575,600,625,650,675,700])
-rirr_in=1e-3*[23.7641,26.5037,27.9743,30.4914,28.1356,
-    21.9385,18.6545,13.5100,5.6338,3.9272,2.9621,2.1865,1.8015]
-
-rrs_out=RemotelySensedReflectance(rirr_in,wvbd_in,wvbd_out)
-
-using Plots
-plot(vec(rrs_out),linewidth=4,lab="recomputed RRS")
-rrs_ref=1e-3*[4.4099, 4.8533, 5.1247, 4.5137, 3.0864, 0.4064]
-plot!(rrs_ref,linewidth=4,ls=:dash,lab="reference result")
-
-```
-"""
-function RemotelySensedReflectance(rirr_in,wvbd_in,wvbd_out)
-    rrs_out=similar(wvbd_out)
-
-    nin=length(wvbd_in)
-    nout=length(wvbd_out)
-
-    jj=Array{Int64,1}(undef,nin)
-    ww=Array{Float64,1}(undef,nin)
-    for ii=1:6
-        tmp=wvbd_out[ii].-wvbd_in
-        kk=maximum(findall(tmp.>=0))
-        jj[ii]=kk
-        ww[ii]=tmp[kk]/(wvbd_in[kk+1]-wvbd_in[kk])
-    end
-
-    tmp=rirr_in/3
-    rrs_tmp=(0.52*tmp)./(1.0 .-1.7*tmp);
-
-    rrs_out=Array{Float64,1}(undef,nout)
-    for vv=1:6
-        tmp0=rrs_tmp[jj[vv]]
-        tmp1=rrs_tmp[jj[vv]+1]
-        rrs_out[vv]=tmp0.*(1-ww[vv])+tmp1.*ww[vv]
-    end
-
-    return rrs_out
-end
+export RemotelySensedReflectance, RrsToChla
+export FuzzyClassification, Moore2009, Jackson2017
 
 end # module

test/runtests.jl

@@ -8,8 +8,18 @@ wvbd_in=[400,425,450,475,500,525,550,575,600,625,650,675,700]
 rirr_in=1e-3*[23.7641,26.5037,27.9743,30.4914,28.1356,
     21.9385,18.6545,13.5100,5.6338,3.9272,2.9621,2.1865,1.8015]
 rrs_ref=1e-3*[4.4099, 4.8533, 5.1247, 4.5137, 3.0864, 0.4064]
+chla_ref=0.6101219875535429
+fcm_ref=[5.683734365402806e-26, 4.958850399758996e-64, 2.0058235898296587e-43,
+  2.399985306394786e-23, 1.0914931954631329e-15, 3.610712128260008e-8,
+  0.06660122938880934, 0.4566726319832639, 0.001153891939362191, 5.901424234631198e-11,
+  1.8433658038301877e-10, 0.4133547888920851, 0.027225883646784382, 0.005595266829573927]
 
 rrs_out=RemotelySensedReflectance(rirr_in,wvbd_in,wvbd_out)
+chla_out=RrsToChla(rrs_out; Eqn="OC4")
+(M,Sinv)=Jackson2017()
+fcm_out=FuzzyClassification(M,Sinv,rrs_out)
 
     @test isapprox(rrs_out,rrs_ref,atol=1e-5)
+    @test isapprox(chla_out,chla_ref,atol=1e-5)
+    @test isapprox(fcm_out,fcm_ref,atol=1e-5)
 end