|
| 1 | +# # Examples from the ESDL paper |
| 2 | +# ## Earth Syst. Dynam., 11, 201–234, 2020 [doi](https://doi.org/10.5194/esd-11-201-2020) |
| 3 | + |
| 4 | +# **NOTE:** This section is based on the case studies from the paper "Earth system data cubes unravel global multivariate dynamics" by Mahecha, Gans et al. (2019). Original scripts are available at https://github.com/esa-esdl/ESDLPaperCode.jl. |
| 5 | +# - We have slightly adjusted the scripts. A few differences are that these new scripts are updated to Julia 1.9, and the YAXArrays.jl package is used. |
| 6 | +# - The dataset has been updated but it has less available variables. Therefore the results might differ. |
| 7 | +# - The calculations are performed with a very coarse spatial (2.5°) and temporal resolution (monthly). |
| 8 | +# - These are examples for illustrative purposes of the packages and do not intend any deeper scientific interpretation. For scientific analysis use the higher spatio-temporal resolution datasets. |
| 9 | + |
| 10 | +# ## Case study 1: Seasonal dynamics on the land surface |
| 11 | +# ### Based on simple seasonal statistics |
| 12 | + |
| 13 | +# * The code is written based on Julia 1.9 |
| 14 | + |
| 15 | +# * Normal text are explanations referring to notation and equations in the paper |
| 16 | + |
| 17 | +# * `# comments in the code are itended explain specific aspects of the coding` |
| 18 | + |
| 19 | +# * **New steps in workflows are introduced with bold headers** |
| 20 | + |
| 21 | +# Load requiered packages |
| 22 | + |
| 23 | +using Pkg |
| 24 | + |
| 25 | +## for plotting later on (need to be loaded first, to avoid conflicts) |
| 26 | +using PyCall, PyPlot, PlotUtils |
| 27 | + |
| 28 | +## for operating data cubes |
| 29 | +using Zarr, YAXArrays |
| 30 | +using EarthDataLab |
| 31 | + |
| 32 | +## for data analysis |
| 33 | +using Statistics, Dates, SkipNan |
| 34 | + |
| 35 | + |
| 36 | +# Next we get a handle to the Earth System Data Cube we want to use, which provides a description of the cube: |
| 37 | +cube_handle = esdc(res="tiny") |
| 38 | + |
| 39 | +## if we want the names of the variables: |
| 40 | +println(getAxis("Var", cube_handle).values) |
| 41 | + |
| 42 | +# We decide which variables to plot: |
| 43 | +vars = ["gross_primary_productivity", "air_temperature_2m", "surface_moisture"] |
| 44 | + |
| 45 | +## time based on data span availability |
| 46 | +time_overlap = Date("2001-01-01")..Date("2020-12-31") |
| 47 | + |
| 48 | +# So we "virtually get" the cube data virtually: |
| 49 | +cube_subset = subsetcube(cube_handle, variable=vars, time=time_overlap) |
| 50 | + |
| 51 | +# 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: |
| 52 | + # $$ |
| 53 | + # f_{\{time\}}^{\{doy\}} : \mathcal{C}(\{lat, lon, time, var\}) \rightarrow \mathcal{C}(\{lat, lon, doy, var\}) |
| 54 | + # $$ |
| 55 | + |
| 56 | +# median seasonal cycle built-in function |
| 57 | +cube_msc = getMedSC(cube_subset) |
| 58 | + |
| 59 | +# The resulting cube `cube_msc` has is of the form $\mathcal{C}(\{lat, lon, doy, var\})$. On this cube we want to apply function `nan_med` (see below) to estimate latitudinal averages for all variables. The atomic function (Eq. 10) needs to |
| 60 | +# have the form, i.e. expecting a longitude and returning a scalar: |
| 61 | +# $$ |
| 62 | +# f_{\{lon\}}^{\{\}} : \mathcal{C}(\{lat, lon, doy, var\}) \rightarrow \mathcal{C}(\{lat, doy, var\}) |
| 63 | +# $$ |
| 64 | + |
| 65 | +import Statistics.median |
| 66 | +cube_msc_lat = mapslices(median ∘ skipmissing, cube_msc, dims = "Lon") |
| 67 | + |
| 68 | +# 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. |
| 69 | + |
| 70 | +# ### Visualization |
| 71 | + |
| 72 | +# At this point we leave the `Earth Data Lab` and and go for visualizations. Using PyPlot we can generate fig. 3 of the paper, the data can be exatracted from the resutls cube via array-access `A[:, :]`. |
| 73 | + |
| 74 | + |
| 75 | +## Create a plot and make it a polar plot |
| 76 | +function zonal_polar_plot(d_msc_lat, sbp, it, vari, lab) |
| 77 | + |
| 78 | + ax = subplot(sbp, polar = "true") |
| 79 | + ourcmap = ColorMap(get_cmap("Spectral_r", 72)) |
| 80 | + |
| 81 | + ## set polar ticks and labels |
| 82 | + month_ang = 0:11 |
| 83 | + month_ang = month_ang .* (360/12) |
| 84 | + month_lab = ["Jan"; "Feb"; "Mar"; "Apr"; "May"; "Jun"; "Jul"; "Aug"; "Sep"; "Oct"; "Nov"; "Dec"] |
| 85 | + |
| 86 | + ## tuning |
| 87 | + ax.set_yticklabels([]) |
| 88 | + ax.set_thetagrids(angles = month_ang, |
| 89 | + labels = month_lab, |
| 90 | + #rotation = month_ang |
| 91 | + ) |
| 92 | + |
| 93 | + ## set Jan to the top of the plot |
| 94 | + ax.set_theta_zero_location("N") |
| 95 | + |
| 96 | + ## switch to clockwise |
| 97 | + ax.set_theta_direction(-1) |
| 98 | + |
| 99 | + ## color setup |
| 100 | + if isequal(vari, "gross_primary_productivity") |
| 101 | + N_min = 0 |
| 102 | + N_max = 8 |
| 103 | + N_var = 8 |
| 104 | + elseif isequal(vari, "air_temperature_2m") |
| 105 | + N_min = -35 |
| 106 | + N_max = 35 |
| 107 | + N_var = 15 |
| 108 | + elseif isequal(vari, "surface_moisture") |
| 109 | + N_min = 0 |
| 110 | + N_max = 50 |
| 111 | + N_var = 10 |
| 112 | + end |
| 113 | + |
| 114 | + ## background setup |
| 115 | + |
| 116 | + ## define time (all could be done more elegantly in julia I guess... ) |
| 117 | + time_ang = range(0, stop = 12, step = 0.25) |
| 118 | + N_time = length(time_ang) |
| 119 | + time_ang = time_ang .* (360/12) |
| 120 | + time_rad = time_ang .* (pi / 180) |
| 121 | + |
| 122 | + ## create a continous var for the background |
| 123 | + N_div = N_time*10 |
| 124 | + y = range(N_min, stop = N_max, length = N_div) |
| 125 | + |
| 126 | + ## grid to fill the polar plot |
| 127 | + xgrid = repeat(time_rad', N_div, 1) |
| 128 | + ygrid = repeat(y, 1, N_time) |
| 129 | + |
| 130 | + ## a grid of NaNs to make an extra colorbar later |
| 131 | + nangrid = zeros(size(ygrid)).*NaN |
| 132 | + levels = range(-90, stop = 90, step = 10) |
| 133 | + ticks = range(-80, stop = 80, step = 20) |
| 134 | + |
| 135 | + axsurf2 = ax.contourf(xgrid, ygrid, ygrid.*NaN, N_max, |
| 136 | + cmap = ourcmap, |
| 137 | + levels = levels, |
| 138 | + ticks = ticks) |
| 139 | + |
| 140 | + ## background to the range of values |
| 141 | + if isequal(vari, "gross_primary_productivity") |
| 142 | + axsurf = ax.contourf(xgrid, ygrid, ygrid, N_max, |
| 143 | + cmap = ColorMap("gray_r"), |
| 144 | + extend = "max") |
| 145 | + elseif isequal(vari, "air_temperature_2m") |
| 146 | + axsurf = ax.contourf(xgrid, ygrid, ygrid, N_max, |
| 147 | + cmap = ColorMap("gray_r"), |
| 148 | + extend = "both") |
| 149 | + elseif isequal(vari, "surface_moisture") |
| 150 | + axsurf = ax.contourf(xgrid, ygrid, ygrid, N_max, |
| 151 | + cmap = ColorMap("gray_r"), |
| 152 | + extend = "max") |
| 153 | + end |
| 154 | + |
| 155 | + ## colorbar setup |
| 156 | + if isodd(sbp) |
| 157 | + ## add forground colorbar |
| 158 | + cban = colorbar(axsurf2, fraction = 0.05, shrink = 0.5, pad = 0.18) |
| 159 | + cban.ax.set_title(label = "Latitude") |
| 160 | + cban.set_ticks([-80, -60, -40, -20, 0, 20, 40, 60, 80]) |
| 161 | + cban.set_ticklabels(["80°N", "60°N", "40°N", "20°N", "0°", "20°S", "40°S", "60°S", "80°S"]) |
| 162 | + cban.ax.invert_yaxis() |
| 163 | + else |
| 164 | + ## add background colorbar |
| 165 | + cbsurf = colorbar(axsurf, fraction = 0.05, shrink = 0.5, pad = 0.18) |
| 166 | + if vari == "gross_primary_productivity" |
| 167 | + cbsurf.ax.set_title(label = "GPP") |
| 168 | + cbsurf.set_label(label = "g C / (m2 d)", rotation=270, labelpad=+20) |
| 169 | + elseif vari == "air_temperature_2m" |
| 170 | + cbsurf.ax.set_title(label = "Tair") |
| 171 | + cbsurf.set_label(label = "°C", rotation = 270, labelpad=+20) |
| 172 | + elseif vari == "surface_moisture" |
| 173 | + cbsurf.ax.set_title(label = "Surf. Moist.") |
| 174 | + cbsurf.set_label(label = "[]", rotation = 270, labelpad=+20) |
| 175 | + cbsurf.set_ticks([0, 10, 20, 30, 40, 50]) |
| 176 | + cbsurf.set_ticklabels(["0", "0.1", "0.2", "0.3", "0.4", "0.5"]) |
| 177 | + end |
| 178 | + end |
| 179 | + |
| 180 | + ## forground setup |
| 181 | + |
| 182 | + ## plot the real data |
| 183 | + N_msc = size(d_msc_lat)[1] |
| 184 | + time_ang_dat = range(1/N_msc, step = 1/N_msc, length = N_msc) |
| 185 | + time_ang_dat = time_ang_dat .* (360) |
| 186 | + time_rad_dat = time_ang_dat .* (pi / 180) |
| 187 | + time_rad_dat = [time_rad_dat; time_rad_dat[1]] |
| 188 | + |
| 189 | + ## add your data |
| 190 | + for j = it |
| 191 | + jj = convert(Int, j) |
| 192 | + try |
| 193 | + var_idx = findall(vari .== getAxis("Variable", d_msc_lat).values)[1] |
| 194 | + ts = d_msc_lat[:, jj, var_idx] |
| 195 | + va = [ts; ts[1]] |
| 196 | + ## correction for temperature |
| 197 | + if isequal(vari, "air_temperature_2m") |
| 198 | + va = va.-273.15 |
| 199 | + elseif isequal(vari, "surface_moisture") |
| 200 | + va = va.*100 |
| 201 | + end |
| 202 | + p = ax.plot(time_rad_dat, va, |
| 203 | + color = ourcmap(jj), |
| 204 | + linewidth = 0.8) |
| 205 | + catch |
| 206 | + end |
| 207 | + end |
| 208 | + |
| 209 | +end |
| 210 | + |
| 211 | +# create a new figure |
| 212 | +f1 = figure("polar_lineplot_new", figsize = (10, 15)) |
| 213 | + |
| 214 | +# get the latitude values for which we have data |
| 215 | +L = collect(getAxis("lat", caxes(cube_msc_lat)).values) |
| 216 | + |
| 217 | +sbps = 321:2:332 |
| 218 | +labtoshow = ["a)", "b)", "c)", "d)", "e)", "f)"] |
| 219 | +vari = getAxis("Variable", caxes(cube_msc_lat)).values |
| 220 | + |
| 221 | +for (sbp, lab, vari) in zip(sbps, labtoshow, vari) |
| 222 | + it1 = range(72/2, stop = 1, step = -2) |
| 223 | + it2 = range(72/2+1, stop = 72, step = 2) |
| 224 | + zonal_polar_plot(cube_msc_lat, sbp, it1, vari, lab) |
| 225 | + zonal_polar_plot(cube_msc_lat, sbp+1, it2, vari, lab) |
| 226 | +end |
| 227 | + |
| 228 | +f1 |
0 commit comments