diff --git a/DataSet From DARWIN/AddRrsCorrection/README.md b/DataSet From DARWIN/AddRrsCorrection/README.md new file mode 100644 index 0000000..2fe7407 --- /dev/null +++ b/DataSet From DARWIN/AddRrsCorrection/README.md @@ -0,0 +1,3 @@ +# Add RrsCorrection.jl will open the DARWIN Dataset from the oceanclustering repository (stored here in this folder ) and apply +# the zenith angle correction to it + diff --git a/DataSet From DARWIN/AddRrsCorrection/RrsCorrection.jl b/DataSet From DARWIN/AddRrsCorrection/RrsCorrection.jl new file mode 100644 index 0000000..8209ade --- /dev/null +++ b/DataSet From DARWIN/AddRrsCorrection/RrsCorrection.jl @@ -0,0 +1,10 @@ +using Plots, Distributions, NetCDF, NCDatasets, Interpolations + +filename="HalfDegree.nc" +varnames=["Rirr412", "Rirr443", "Rirr490", "Rirr510", "Rirr555", "Rirr670"] +for n in 1:length(varnames) + Rirr=ncread(filename, varnames[n]) + tmp=Rirr/3 + Rrs0=(0.52*tmp)./(1.0 .-1.7*tmp) + ncwrite(Rrs0,"HalfDegree.nc",varnames[n]) +end diff --git a/DataSet From DARWIN/Main_MakeGridDarwin2.jl b/DataSet From DARWIN/Main_MakeGridDarwin2.jl index 2c23e91..f37c905 100644 --- a/DataSet From DARWIN/Main_MakeGridDarwin2.jl +++ b/DataSet From DARWIN/Main_MakeGridDarwin2.jl @@ -3,10 +3,10 @@ include("MakeGridDarwin2.jl") -if ~isfile("MasterHalf.nc") +#if ~isfile("MasterHalf.nc") MakeMasterHalf() -end -realnames=["Latitude","Longitude","Month","THETA", "SALT", "Chl", "Rirr412", "Rirr443", "Rirr490", "Rirr510", "Rirr555","Rirr670","MXLDEPTH","Bottom Depth", "Wind Speed", "TKE","PAR","Euphotic Depth"] -for n in [.5 1 2 4] -MakeDarwinGrid(n,realnames) -end \ No newline at end of file +#end +# realnames=["Latitude","Longitude","Month","THETA", "SALT", "Chl", "Rirr412", "Rirr443", "Rirr490", "Rirr510", "Rirr555","Rirr670","MXLDEPTH","Bottom Depth", "Wind Speed", "TKE","PAR","Euphotic Depth"] +# for n in [.5 1 2 4] +# MakeDarwinGrid(n,realnames) +# end \ No newline at end of file diff --git a/DataSet From DARWIN/MakeGridDarwin2.jl b/DataSet From DARWIN/MakeGridDarwin2.jl index e14e0b5..f6ab097 100644 --- a/DataSet From DARWIN/MakeGridDarwin2.jl +++ b/DataSet From DARWIN/MakeGridDarwin2.jl @@ -9,6 +9,8 @@ function MakeMasterHalf() #This is the starting folder for the Climatology #name="C:\\Users\\folle\\Dropbox (MIT)\\201805-CBIOMES-climatology\\interp" name="C:\\Users\\folle\\Documents\\CBIOMES\\interp" + #name="/nfs/cnhlab001/cnh/cbiomes/201805-CBIOMES-climatology/interp/" + #This is the filename for the .5 degree monthly climatology filename="HalfDegree.nc" #These are the variable names which we would like on the grid @@ -121,7 +123,10 @@ function MakeMasterHalf() OCCIbands=["Rirr412", "Rirr443", "Rirr490", "Rirr510", "Rirr555", "Rirr670"] #export OCCIbands - if ~isfile("InterpedWavebands.nc") + if isfile("InterpedWavebands3.nc") + rm("InterpedWavebands3.nc") + end + if ~isfile("InterpedWavebands3.nc") for n in 1:6 nccreate("InterpedWavebands3.nc", OCCIbands[n],"Lon",collect(-179.75:.5:179.75),"Lat",collect(-89.75:.5:89.75),"Month",collect(1:12)) ncwrite(wbcci[:,:,:,n],"InterpedWavebands3.nc", OCCIbands[n]) @@ -191,7 +196,7 @@ function MakeMasterHalf() if isfile("EuphoticDepth.nc") rm("EuphoticDepth.nc") end - nccreate("EuphoticDepth.nc", "Euphotic Depth","Lon",collect(.5:.5:360),"Lat",collect(-89.75:.5:89.75),"Month",collect(1:12)) + nccreate("EuphoticDepth.nc", "Euphotic Depth","Lon",collect(-179.75:.5:179.75),"Lat",collect(-89.75:.5:89.75),"Month",collect(1:12)) ncwrite(Eudepth,"EuphoticDepth.nc", "Euphotic Depth") #lon lat depth t @@ -207,8 +212,8 @@ function MakeMasterHalf() #We have everything except PAR, Ekman pumping, Kd, and Euphotic Depth. vars=["THETA", "SALT", "Chl", "Rirr412", "Rirr443", "Rirr490", "Rirr510", "Rirr555", "Rirr670","MXLDEPTH","Bottom Depth", "EXFwspee", "GGL90TKE","PAR","Euphotic Depth"] - realnames=["THETA", "SALT", "Chl", "Rirr412", "Rirr443", "Rirr490", "Rirr510", "Rirr555", - "Rirr670","MXLDEPTH","Bottom Depth", "Wind Speed", "TKE","PAR","Euphotic Depth"] + realnames=["SST", "SALT", "Chl", "Rrs412", "Rrs443", "Rrs490", "Rrs510", "Rrs555", + "Rrs670","MLD","Bathymetry", "wind", "EKE","PAR","Euphotic Depth"] #export vars Tg=Array{Float64,4}(undef,720,360,12,length(vars)) for n in [1 2 3] @@ -255,7 +260,7 @@ function MakeDarwinGrid(deg,names) for n in 1:length(names) if sum(k.==names[n])>=1 vartmp=ncread("MasterHalf.nc",names[n]) - nccreate("DarwinGrid"*string(deg)*".nc", names[n],"Londim",collect(.5:deg:360),"Latdim",collect(-89.75:deg:89.75),"Monthdim",collect(1:12)) + nccreate("DarwinGrid"*string(deg)*".nc", names[n],"Londim",collect(-179.75:deg:179.75),"Latdim",collect(-89.75:deg:89.75),"Monthdim",collect(1:12)) #println(size(vartmp[1:Int(deg/.5):end,1:Int(deg/.5):end,:])) #println(size(vartmp)) tmp=vartmp[1:Int(deg/.5):end,1:Int(deg/.5):end,:]; @@ -268,7 +273,7 @@ function MakeDarwinGrid(deg,names) ncwrite(tmp,"DarwinGrid"*string(deg)*".nc", names[n]) end if sum(k.==names[n])==0 - tmp=zeros(length(collect(.5:deg:360)), length(collect(-89.75:deg:89.75)),12) + tmp=zeros(length(collect(-179.75:deg:179.75)), length(collect(-89.75:deg:89.75)),12) if names[n]=="Latitude" Lat=collect(-89.75:deg:89.75); for m in 1:length(Lat) @@ -276,7 +281,7 @@ function MakeDarwinGrid(deg,names) tmp2=fill(Lat[m],size(tmp[:,m,:])) tmp[:,m,:]=tmp2 end - nccreate("DarwinGrid"*string(deg)*".nc", names[n],"Londim",collect(.5:deg:360),"Latdim",collect(-89.75:deg:89.75),"Monthdim",collect(1:12)) + nccreate("DarwinGrid"*string(deg)*".nc", names[n],"Londim",collect(-179.75:deg:179.75),"Latdim",collect(-89.75:deg:89.75),"Monthdim",collect(1:12)) ncwrite(tmp,"DarwinGrid"*string(deg)*".nc", names[n]) sztmp=size(tmp) end diff --git a/DataSet From DARWIN/README.md b/DataSet From DARWIN/README.md index 74a8af3..ae8a294 100644 --- a/DataSet From DARWIN/README.md +++ b/DataSet From DARWIN/README.md @@ -14,19 +14,19 @@ This folder contains code written in Julia to generate gridded data for the foll - "Lat" (degrees) - "Lon" (degrees) - "Month" (month) -- "THETA" (degrees C) +- "SST" (degrees C) - "SALT" (psu) - "Chl" (mg/l) -- "Rirr412" (irradiance reflectance for waveband at wavelength) -- "Rirr443" (irradiance reflectance for waveband at wavelength) -- "Rirr490" (irradiance reflectance for waveband at wavelength) -- "Rirr510" (irradiance reflectance for waveband at wavelength) -- "Rirr555" (irradiance reflectance for waveband at wavelength) -- "Rirr670" (irradiance reflectance for waveband at wavelength) -- "MXLDEPTH" (depth, meters) -- "Bottom Depth" (depth, meters) -- "Wind Speed" (m/s) -- "TKE" (m^2/s^2) is GGL90 sub-grid turbulent kinetic energy +- "Rrs412" (irradiance reflectance for waveband at wavelength) +- "Rrs443" (irradiance reflectance for waveband at wavelength) +- "Rrs490" (irradiance reflectance for waveband at wavelength) +- "Rrs510" (irradiance reflectance for waveband at wavelength) +- "Rrs555" (irradiance reflectance for waveband at wavelength) +- "Rrs670" (irradiance reflectance for waveband at wavelength) +- "MLD" (depth, meters) +- "Bathymetry" (depth, meters) +- "wind" (m/s) +- "EKE" (m^2/s^2) is GGL90 sub-grid turbulent kinetic energy - "PAR" () - "Euphotic Depth" Depth (m) where light is 1% of the surface values diff --git a/Try Jackson/J17.nc b/Try Jackson/J17.nc new file mode 100644 index 0000000..c74da93 Binary files /dev/null and b/Try Jackson/J17.nc differ diff --git a/Try Jackson/JacksonModelRemote.jl b/Try Jackson/JacksonModelRemote.jl new file mode 100644 index 0000000..03e4698 --- /dev/null +++ b/Try Jackson/JacksonModelRemote.jl @@ -0,0 +1,92 @@ +#This script will attempt to apply the Jackson classifier to the DARWIN reflectance data + +using Plots, Distributions, NetCDF, NCDatasets, Interpolations + + +##This is Gael's Code, taken from 05-OceanColourClassifiers.ipynb## + +function fcm(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 + f +end + +function getJackson(file) + #Jackson et al 2017: + tmpM = ncread("J17.nc", "cluster_means") + tmpSinv = ncread("J17.nc", "inverse_covariance") + wv_cci=[412, 443, 490, 510, 555, 670] + 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 + + J17=Dict("M" => M, "Sinv" => Sinv, "S" => inv.(Sinv)) + #plot(wv_cci,M,w=3); xlabel!("nm"); ylabel!("Rrs") + + + ds = Dataset(file) + + Rrs_412=ds["Rrs412"] + Rrs_443=ds["Rrs443"] + Rrs_490=ds["Rrs490"] + Rrs_510=ds["Rrs510"] + Rrs_555=ds["Rrs555"] + Rrs_670=ds["Rrs670"] + + Rrs_412=Rrs_412[1:4:720,1:4:360,1:4:12] + Rrs_443=Rrs_443[1:4:720,1:4:360,1:4:12] + Rrs_490=Rrs_490[1:4:720,1:4:360,1:4:12] + Rrs_510=Rrs_510[1:4:720,1:4:360,1:4:12] + Rrs_555=Rrs_555[1:4:720,1:4:360,1:4:12] + Rrs_670=Rrs_670[1:4:720,1:4:360,1:4:12] + + tmp=fill(false,size(Rrs_412)) + for ii in eachindex(Rrs_412) + !ismissing(Rrs_412[ii]) ? tmp[ii]=!ismissing(Rrs_443[ii].*Rrs_490[ii].*Rrs_510[ii].*Rrs_555[ii].*Rrs_670[ii]) : nothing + end + ii=findall(tmp) + sz=size(Rrs_412) + #export ii, sz, Rrs_412, Rrs_443, Rrs_490, Rrs_510, Rrs_555, Rrs_670 + mbrshp=Array{Float64,4}(undef,(sz[1],sz[2],sz[3],14)); + for jj=1:length(ii); + kk=ii[jj] + Rrs_tmp=[Rrs_412[kk] Rrs_443[kk] Rrs_490[kk] Rrs_510[kk] Rrs_555[kk] Rrs_670[kk]] + mbrshp[kk[1],kk[2],kk[3],:]=fcm(J17["M"],J17["Sinv"],Rrs_tmp) + println(kk) + end + return mbrshp +end +#path="https://rsg.pml.ac.uk/shared_files/brj/CBIOMES_ecoregions/ver_0_2_5" +path="https://rsg.pml.ac.uk/shared_files/brj/CBIOMES_ecoregions/ver_0_2_5/" +filegeo="gridded_geospatial_montly_clim_360_720_ver_0_2.nc" +filedarwin="gridded_darwin_montly_clim_360_720_ver_0_2_2.nc" +#run(`wget $path/$file`) +if ~isfile(filegeo) + download(string(path, filegeo), filegeo) +end +if ~isfile(filedarwin) + download(string(path, filedarwin), filedarwin) +end + +#varnames=["Rirr412", "Rirr443", "Rirr490", "Rirr510", "Rirr555", "Rirr670"] +varnames=["Rrs412", "Rrs443", "Rrs490", "Rrs510", "Rrs555", "Rrs670"] + +for n in 1:length(varnames) + Rirr=ncread(filedarwin, varnames[n]) + tmp=Rirr/3 + Rrs0=(0.52*tmp)./(1.0 .-1.7*tmp) + ncwrite(Rrs0,filedarwin,varnames[n]) + #renameVar(filedarwin,varnames[n],newnames[n]) +end + +DarwinMembers=getJackson(filedarwin) +GeoMembers=getJackson(filegeo) +export DarwinMembers +export GeoMembers \ No newline at end of file diff --git a/Try Jackson/M09.jl b/Try Jackson/M09.jl new file mode 100644 index 0000000..a74178e --- /dev/null +++ b/Try Jackson/M09.jl @@ -0,0 +1,63 @@ + +#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]; + diff --git a/Try Jackson/README.md b/Try Jackson/README.md new file mode 100644 index 0000000..338ceaf --- /dev/null +++ b/Try Jackson/README.md @@ -0,0 +1,7 @@ +# TryJackson + +This folder contains code that will download the half degree data grids for DARWIN and GEO and attempt to apply the OC-CCI water class algorithm to it. + +The important code is JacksonModelRemote.jl. In the current version, the Rirr to Rrs correction is applied to the DARWIN dataset still. + + diff --git a/Try Jackson/cci_Rrs_490.bin b/Try Jackson/cci_Rrs_490.bin new file mode 100644 index 0000000..cd556c1 Binary files /dev/null and b/Try Jackson/cci_Rrs_490.bin differ diff --git a/Try Jackson/drwn3_Rirr.bin b/Try Jackson/drwn3_Rirr.bin new file mode 100644 index 0000000..14521c4 Binary files /dev/null and b/Try Jackson/drwn3_Rirr.bin differ diff --git a/Try Jackson/lat.bin b/Try Jackson/lat.bin new file mode 100644 index 0000000..6ad0c02 Binary files /dev/null and b/Try Jackson/lat.bin differ diff --git a/Try Jackson/lon.bin b/Try Jackson/lon.bin new file mode 100644 index 0000000..a0adbac Binary files /dev/null and b/Try Jackson/lon.bin differ diff --git a/Try Jackson/tryJackson.jl b/Try Jackson/tryJackson.jl new file mode 100644 index 0000000..44269a8 --- /dev/null +++ b/Try Jackson/tryJackson.jl @@ -0,0 +1,85 @@ +#This script will attempt to apply the Jackson classifier to the DARWIN reflectance data + +using Plots, Distributions, NetCDF, NCDatasets, Interpolations + + +##This is Gael's Code, taken from 05-OceanColourClassifiers.ipynb## + +function fcm(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 + f +end + +#Jackson et al 2017: +tmpM = ncread("J17.nc", "cluster_means") +tmpSinv = ncread("J17.nc", "inverse_covariance") +wv_cci=[412, 443, 490, 510, 555, 670] +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 + +J17=Dict("M" => M, "Sinv" => Sinv, "S" => inv.(Sinv)) +plot(wv_cci,M,w=3); xlabel!("nm"); ylabel!("Rrs") + +#path="https://rsg.pml.ac.uk/shared_files/brj/CBIOMES_ecoregions/ver_0_2_5" +path="https://https://rsg.pml.ac.uk/shared_files/brj/CBIOMES_ecoregions/ver_0_2_5/" +file="gridded_geospatial_montly_clim_360_720_ver_0_2.nc" +file="gridded_darwin_montly_clim_360_720_ver_0_2_2.nc" +#run(`wget $path/$file`) +download(string(path, file), file) +#Using NCDatasets +#Dataset(file) +#fil="C:\Users\folle\Dropbox (MIT)\CBIOMES_TransitionZone\Chris Follett\StatsGroup\ocean_clustering\DataSet From DARWIN\AddRrsCorrection\HalfDegree.nc" +#fil="C:\\Users\\folle\\Dropbox (MIT)\\CBIOMES_TransitionZone\\Chris Follett\\StatsGroup\\ocean_clustering\\DataSet From DARWIN\\AddRrsCorrection\\HalfDegree.nc" +#dir0="../samples/" +#fil="MasterHalf.nc" +## +varnames=["Rirr412", "Rirr443", "Rirr490", "Rirr510", "Rirr555", "Rirr670"] +for n in 1:length(varnames) + Rirr=ncread(filename, varnames[n]) + tmp=Rirr/3 + Rrs0=(0.52*tmp)./(1.0 .-1.7*tmp) + ncwrite(Rrs0,"HalfDegree.nc",varnames[n]) +end +## + +ds = Dataset(file) + +Rrs_412=ds["Rirr412"] +Rrs_443=ds["Rirr443"] +Rrs_490=ds["Rirr490"] +Rrs_510=ds["Rirr510"] +Rrs_555=ds["Rirr555"] +Rrs_670=ds["Rirr670"] + +Rrs_412=Rrs_412[1:4:720,1:4:360,1:4:12] +Rrs_443=Rrs_443[1:4:720,1:4:360,1:4:12] +Rrs_490=Rrs_490[1:4:720,1:4:360,1:4:12] +Rrs_510=Rrs_510[1:4:720,1:4:360,1:4:12] +Rrs_555=Rrs_555[1:4:720,1:4:360,1:4:12] +Rrs_670=Rrs_670[1:4:720,1:4:360,1:4:12] + +tmp=fill(false,size(Rrs_412)) +for ii in eachindex(Rrs_412) + !ismissing(Rrs_412[ii]) ? tmp[ii]=!ismissing(Rrs_443[ii].*Rrs_490[ii].*Rrs_510[ii].*Rrs_555[ii].*Rrs_670[ii]) : nothing +end +ii=findall(tmp) +sz=size(Rrs_412) +export ii, sz, Rrs_412, Rrs_443, Rrs_490, Rrs_510, Rrs_555, Rrs_670 +mbrshp=Array{Float64,4}(undef,(sz[1],sz[2],sz[3],14)); +for jj=1:length(ii); + kk=ii[jj] + Rrs_tmp=[Rrs_412[kk] Rrs_443[kk] Rrs_490[kk] Rrs_510[kk] Rrs_555[kk] Rrs_670[kk]] + mbrshp[kk[1],kk[2],kk[3],:]=fcm(J17["M"],J17["Sinv"],Rrs_tmp) + println(kk) +end +export mbrshp +heatmap(mbrshp[:,:,1,10]) \ No newline at end of file