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