Corrections

Project.toml

@@ -3,6 +3,7 @@ uuid = "2502ad13-e485-5df4-bc88-d00dc12670bb"
 
 [deps]
 AWSCore = "4f1ea46c-232b-54a6-9b17-cc2d0f3e6598"
+Blosc = "a74b3585-a348-5f62-a45c-50e91977d574"
 Conda = "8f4d0f93-b110-5947-807f-2305c1781a2d"
 Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
@@ -16,5 +17,7 @@ Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0"
 PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
+Shapefile = "8e980c4a-a4fe-5da2-b3a7-4b4b0353a2f4"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+WeightedOnlineStats = "bbac0a1f-7c9d-5672-960b-c6ca726e5d5d"
 Zarr = "0a941bbe-ad1d-11e8-39d9-ab76183a1d99"

docs/make.jl

@@ -1,17 +1,16 @@
-using Pkg
-Pkg.activate(joinpath(@__DIR__,".."))
 using Documenter
 
 makedocs(
   format = Documenter.HTML(),
-  sitename = "ESDLExampleCode.jl",
+  sitename = "ESDLPaperCode.jl",
   authors = "Miguel Mahecha, Fabian Gans",
   pages = [
     "Home" => "index.md",
     "Case studies" => [
       "Seasonality" => "ESDL case study 1 seasonality.md",
-      "Intrinsic Dimensions" => "ESDL case study 2 intrinsic dimension.md",
+      "Intrinsic Dimensions" => "ESDL case study 2 Intrinsic dimension.md",
       "Q10" => "ESDL case study 3 q10.md",
+      "Polygons" => "ESDL case study 4 Polygons.md",
       ]
   ],
   workdir = @__DIR__

docs/src/ESDL case study 1 seasonality.md

@@ -1,5 +1,5 @@
 ```@meta
-EditURL = "@__REPO_ROOT_URL__/"
+EditURL = "<unknown>/src/ESDL case study 1 seasonality.jl"
 ```
 
 ## Case study 1: Seasonal dynamics on the land surface
@@ -27,15 +27,11 @@ using PyCall, PyPlot, PlotUtils
 
 # for operating the Earth system data lab
 using ESDL
-
-# for parallel computing
-using Distributed
 ```
 
 Next we get a handle to the Earth System Data Cube we want to use, which provides a description of the cube:
 
 ```@example ESDL case study 1 seasonality
-cd(@__DIR__)
 cube_handle = Cube("../data/subcube")
 ```
 
@@ -62,11 +58,6 @@ So we "virtually get" the cube data virually:
 cube_subset = subsetcube(cube_handle, variable = vars)
 ```
 
-```@example ESDL case study 1 seasonality
-# The next step requires a bit of CPU -> add some parallel processors:
-addprocs(4)
-```
-
 The next function estimates the median seasonal cycle. This changes the dimension of the cube, as the time domain is replaced by day of year (doy); Eq. 9 in the manuscript:
 ```math
     f_{\{time\}}^{\{doy\}} : \mathcal{C}(\{lat, lon, time, var\}) \rightarrow \mathcal{C}(\{lat, lon, doy, var\})
@@ -83,18 +74,12 @@ The resulting cube `cube_msc` has is of the form $\mathcal{C}(\{lat, lon, doy, v
 ```
 
 ```@example ESDL case study 1 seasonality
-# The median function; the @everywhere brings the contents to each core
-@everywhere import Statistics.median
+import Statistics.median
 
 # Applied to the dimension "Lon"
 cube_msc_lat = mapslices(median ∘ skipmissing, cube_msc, dims = "Lon")
 ```
 
-```@example ESDL case study 1 seasonality
-# Now the hard work is done and we can remove the workers
-rmprocs(workers())
-```
-
 The result of each operation on a data cube is a data cube. Here the resulting cube has the form $\mathcal{C}(\{doy, lat, var\})$
 as expected but in different order, which is, irrelevant as axes have no natural order.
 
@@ -267,10 +252,13 @@ for (sbp, lab, vari) in zip(sbps,labtoshow,getAxis(VariableAxis, caxes(cube_msc_
 
 end
 
-
+mkpath("../figures")
 savefig("../figures/zonalmeans.png",
         bbox_inches = "tight");
+nothing #hide
 ```
 
+---
+
 *This page was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*
 

docs/src/ESDL case study 2 Intrinsic dimension.md

@@ -1,5 +1,5 @@
 ```@meta
-EditURL = "@__REPO_ROOT_URL__/"
+EditURL = "<unknown>/src/ESDL case study 2 Intrinsic dimension.jl"
 ```
 
 ## Case study 2: Intrinsic dimensions of ecosystem dynamics
@@ -27,9 +27,8 @@ using PyCall, PyPlot
 
 # for operating the Earth system data lab
 using ESDL
-
-# for parallel computing
-using Distributed
+import Blosc
+Blosc.set_num_threads(Threads.nthreads())
 ```
 
 In this study we investigate the redundancy the different variables in each pixel. Therefore we calculate a linear dimensionality reduction (PCA) and check how many dimensions are needed to explain 90% of the variance of a cube that contained originally 6 variables.  First we check out the variables from the cube and add some processors, because we want to do a global study
@@ -39,7 +38,6 @@ In this study we investigate the redundancy the different variables in each pixe
 We need to choose a cube and here select a 8-dayly, 0.25° resolution global cube. The cube name suggests it is chunked such that we have one time chunk and 720x1440 spatial chunks
 
 ```@example ESDL case study 2 Intrinsic dimension
-cd(@__DIR__)
 cube_handle = Cube("../data/subcube")
 ```
 
@@ -87,51 +85,47 @@ vars = ["evaporative_stress",
 timespan = Date("2003-01-01")..Date("2011-12-31")
 
 # subset the grand cube and get the cube we will analyse here
-cube_subset = subsetcube(cube_handle, time = timespan, variable = vars)
+cube_subset = subsetcube(cube_handle, time = timespan, variable = vars, region="Europe")
 ```
 
 An important preprocessing step is gapfilling. We do not want to enter the debate on the optimal gapfilling method. What we do here is gapfilling first with the mean seasonal cycle (where it can be estimated), and interpolating long-recurrent gaps (typically in winter seasons).
 
 ```@example ESDL case study 2 Intrinsic dimension
-# gapfilling this requires a bit of CPU -> add some parallel processors:
-addprocs(4)
-
 # use the ESDL buit-in function
 cube_fill = gapFillMSC(cube_subset)
 ```
 
-The interpolation of wintergpas needs a function that we code here an call `LinInterp`.
+The interpolation of wintergaps needs a function that we code here an call `LinInterp`.
 
 ```@example ESDL case study 2 Intrinsic dimension
-# Function LinInterp should be available on every core, i.e. @everywhere
-@everywhere begin
-
-    # package on each core
-    using Interpolations
-
-    function LinInterp(y)
-        try
-            # find the values we need to input
-            idx_nan = findall(ismissing, y)
-            idx_ok  = findall(!ismissing, y)
-
-            # make sure to have a homogenous input array
-            y2 = Float32[y[i] for i in idx_ok]
-
-            # generate an interpolation object based on the good data
-            itp = interpolate((idx_ok,), y2, Gridded(Linear()))
-
-            # fill the missing values based on a linter interpolation
-            y[idx_nan] = itp(idx_nan)
-            return y
-        catch
-            return y
-        end
-    end
+using Interpolations
+
+function LinInterp(y)
+
+  try
+    # find the values we need to input
+    idx_nan = findall(ismissing, y)
+    idx_ok  = findall(!ismissing, y)
+
+    # make sure to have a homogenous input array
+    y2 = Float64[y[i] for i in idx_ok]
+
+    # generate an interpolation object based on the good data
+    itp = extrapolate(interpolate((idx_ok,), y2, Gridded(Linear())),Flat())
+
+    # fill the missing values based on a linter interpolation
+    y[idx_nan] = itp(idx_nan)
+    return y
+  catch
+    idx_nan = findall(ismissing, y)
+    y[idx_nan] .= mean(skipmissing(y))
+    return y
+  end
 end
 
+
 # short test
-x = [2.5,missing,3.8,missing,8.9]
+x = [2.5,missing,3.8,missing,8.9,missing]
 LinInterp(x)
 ```
 
@@ -153,11 +147,13 @@ As we describe in the paper, we estimate the intrinsic dimensions from the raw,
   f_{\{time\}}^{\{time, freq\}} : \mathcal{C}(\{lat, lon, time, var\}) \rightarrow \mathcal{C}(\{lat, lon, time, var, freq\}).
 \end{equation}
 
-which can be done using a pre-implemented ESDL function:
+which can be done using a pre-implemented ESDL function. Note that this step
+will use a lot of computing time.
 
 ```@example ESDL case study 2 Intrinsic dimension
 import Zarr
-cube_decomp = filterTSFFT(cube_fill_itp, compressor=Zarr.BloscCompressor(clevel=1))
+ESDL.ESDLDefaults.compressor[] = Zarr.BloscCompressor(clevel=1)
+cube_decomp = filterTSFFT(cube_fill_itp)
 ```
 
 ### Estimate intrinic dimension via PCA
@@ -185,22 +181,18 @@ We can now apply this to the cube: The latter was the operation described in the
 \end{equation}
 
 ```@example ESDL case study 2 Intrinsic dimension
-# Function sufficient_dimensions should be available on every core, i.e. @everywhere
-@everywhere begin
-
-    # packages needed on each core
-    using MultivariateStats, Statistics
-
-    function sufficient_dimensions(xin::AbstractArray, expl_var::Float64 = 0.95)
-
-        any(ismissing,xin) && return NaN
-        npoint, nvar = size(xin)
-        means = mean(xin, dims = 1)
-        stds  = std(xin,  dims = 1)
-        xin   = broadcast((y,m,s) -> s>0.0 ? (y-m)/s : one(y), xin, means, stds)
-        pca = fit(PCA, xin', pratio = 0.999, method = :svd)
-        return findfirst(cumsum(principalvars(pca)) / tprincipalvar(pca) .> expl_var)
-    end
+# packages needed on each core
+using MultivariateStats, Statistics
+
+function sufficient_dimensions(xin::AbstractArray, expl_var::Float64 = 0.95)
+
+  any(ismissing,xin) && return NaN
+  npoint, nvar = size(xin)
+  means = mean(xin, dims = 1)
+  stds  = std(xin,  dims = 1)
+  xin   = broadcast((y,m,s) -> s>0.0 ? (y-m)/s : one(y), xin, means, stds)
+  pca = fit(PCA, xin', pratio = 0.999, method = :svd)
+  return findfirst(cumsum(principalvars(pca)) / tprincipalvar(pca) .> expl_var)
 end
 ```
 
@@ -213,7 +205,7 @@ cube_int_dim = mapslices(sufficient_dimensions, cube_fill_itp, 0.95, dims = ("Ti
 Saving intermediate results can save CPU later, not needed to guarantee reproducability tough
 
 ```@example ESDL case study 2 Intrinsic dimension
-saveCube(cube_int_dim, "../data/IntDim")
+saveCube(cube_int_dim, "../data/IntDim", overwrite=true)
 ```
 
 Now we apply the same function
@@ -233,7 +225,7 @@ cube_int_dim_dec = mapslices(sufficient_dimensions, cube_decomp, 0.95, dims = ("
 ```
 
 ```@example ESDL case study 2 Intrinsic dimension
-saveCube(cube_int_dim_dec, "../data/IntDimDec")
+saveCube(cube_int_dim_dec, "../data/IntDimDec", overwrite=true)
 ```
 
 ### Visualizing results is not part of the ESDL package.
@@ -325,7 +317,7 @@ for i in 1:4
 
     name = prelab[i]*scale_name[i]
     fig = plot_robin(prelab[i]*scale_name[i], DAT, "Intrinsic DImensions")
-
+    mkpath("../figures")
     savefig("../figures/IntDim_" * save_name[i] * ".pdf",
         orientation = "landscape",
         bbox_inches = "tight")
@@ -396,5 +388,7 @@ xlabel("Intrinsic dimension", fontsize = 14)
 savefig("../figures/IntDim_Hist.pdf", bbox_inches = "tight")
 ```
 
+---
+
 *This page was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*
 

docs/src/ESDL case study 3 q10.md

@@ -1,5 +1,5 @@
 ```@meta
-EditURL = "@__REPO_ROOT_URL__/"
+EditURL = "<unknown>/src/ESDL case study 3 q10.jl"
 ```
 
 ## Case study 3: Model-parameter estimation in the ESDL
@@ -35,7 +35,6 @@ using Statistics
 We need to choose a cube and here select a 8-dayily, 0.25° resolution global cube. The cube name suggests it is chunked such that we have one time chunk and 720x1440 spatial chunks
 
 ```@example ESDL case study 3 q10
-cd(@__DIR__)
 cube_handle = Cube("../data/subcube")
 ```
 
@@ -207,6 +206,7 @@ end
 ```
 
 ```@example ESDL case study 3 q10
+mkpath("../figures")
 p1 = plot_robin("a) Confounded Parameter Estimation", q10_direct[:,:], "Q10",dmin=1, dmax=2.5);
 savefig("../figures/q10_confounded.pdf",orientation="landscape",bbox_inches="tight")
 p2 = plot_robin("b) Scale Dependent Parameter Estimation", q10_scape[:,:], "Q10",dmin=1,dmax=2.5);
@@ -226,7 +226,10 @@ q10_diff = q10_direct - q10_scape
 ```@example ESDL case study 3 q10
 plot_robin("Correlation Tair and Rb",cor_tair_rb[:,:],"Coefficient",dmin=-1.0, dmax=1.0)
 plot_robin("Ratio of Q10conv and Q10Scape",q10_diff[:,:],"Ratio",dmin=-1.0,dmax=1.0);
+nothing #hide
 ```
 
+---
+
 *This page was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*