Add code changes from revision

Makefile

@@ -2,11 +2,14 @@
 NSAMPLES := 1000
 
 # How many ensemble members in the sampled ensemble?
-SAMPLE_SIZE := 75
+SAMPLE_SIZE := 100
 
 # Threshold for the correlation distribution width below which datapoints are
 # considered to be statistically significant
-SIGNIFICANCE_LEVEL := 0.1
+SIGNIFICANCE_LEVEL := 0.09
+
+# How many forecast members per cluster?
+NCLUSTER := 15
 
 # Determine filename suffix for C-extensions to installed Python
 PYEXT := $(shell python3-config --extension-suffix)
@@ -18,14 +21,15 @@ FIGURES := \
 		figures/event24-plume.pdf \
 		figures/event24-budget.pdf \
 		figures/event24-evaluation.pdf \
+		figures/idealized.pdf \
 		figures/event24-maps+hovmoeller.pdf \
 		figures/event24-cluster+scatter.pdf \
 		figures/event24-maps-separate.pdf \
 		figures/event24-cluster-f2.pdf \
 		figures/event24-nh18fig4.pdf
 
 
-.PHONY: all clean local
+.PHONY: all clean
 
 all: $(FIGURES)
 
@@ -38,6 +42,28 @@ clean:
 	rm -f scripts/hn2016_falwa/*.so
 
 
+# Simulations with the traffic jam model
+
+SIMULATIONS := \
+	data/idealized-0.40-60-2.0.nc \
+	data/idealized-0.40-55-2.0.nc \
+	data/idealized-0.40-65-2.0.nc \
+	data/idealized-0.40-60-1.8.nc \
+	data/idealized-0.40-60-2.2.nc
+
+figures/idealized.pdf: $(SIMULATIONS) scripts/plot-idealized.py
+	python3 -m scripts.plot-idealized $(SIMULATIONS) --output $@ --data-output "data/idealized.json"
+
+data/idealized-%-60-2.0.nc: scripts/run-idealized.py
+	python3 -m scripts.run-idealized $@ --alpha="$*" --urefcg="60" --forcing-peak="-2.0"
+
+data/idealized-0.40-%-2.0.nc: scripts/run-idealized.py
+	python3 -m scripts.run-idealized $@ --alpha="0.40" --urefcg="$*" --forcing-peak="-2.0"
+
+data/idealized-0.40-60-%.nc: scripts/run-idealized.py
+	python3 -m scripts.run-idealized $@ --alpha="0.40" --urefcg="60" --forcing-peak="-$*"
+
+
 # Python scipts for data processing and plotting, C-extensions
 
 scripts/calculate-lwa.py: scripts/common/regrid.py scripts/hn2016_falwa/oopinterface.py
@@ -86,7 +112,7 @@ scripts/hn2016_falwa/%$(PYEXT): scripts/hn2016_falwa/%.f90
 
 # Write processed data to netcdf files:
 event24_ana := data/ERA-2016-12-10-to-2016-12-24-lwa-qg.nc
-event24_ens := data/ENS-2016-12-10T00Z-lwa-qg.nc data/ENS-2016-12-10T12Z-lwa-qg.nc data/ENS-2016-12-11T00Z-lwa-qg.nc
+event24_ens := data/ENS-2016-12-09T12Z-lwa-qg.nc data/ENS-2016-12-10T00Z-lwa-qg.nc data/ENS-2016-12-10T12Z-lwa-qg.nc data/ENS-2016-12-11T00Z-lwa-qg.nc
 event24_eval := data/EVAL-2016-12-18T00Z-lwa-qg.nc
 
 # Don't remove these intermediate files
@@ -102,7 +128,7 @@ data/ENS-%-lwa-qg.nc: scripts/calculate-lwa.py data/IFS-ENS/ENS-%/ENS-*-[uvt].nc
 	python3 -m scripts.calculate-lwa --half-resolution $(filter-out $<,$^)
 
 # Forecast evaluation data for 18 Dec 00 UTC
-data/EVAL-2016-12-18T00Z-lwa-qg.nc: scripts/calculate-lwa.py data/IFS-ENS/EVAL/ENS-DEC18-EVAL-*.nc | data
+data/EVAL-2016-12-18T00Z-lwa-qg.nc: scripts/calculate-lwa.py data/IFS-ENS/EVAL/ENS-DEC18-EVAL-*.nc
 	python3 -m scripts.calculate-lwa --half-resolution --eval data/IFS-ENS/EVAL/ENS-DEC18-EVAL-*.nc
 
 
@@ -114,40 +140,51 @@ figures/forecast-combination.pdf: figures/forecast-combination/forecast-combinat
 	cd $(dir $<) && lualatex $(notdir $<)
 	mv $(<:.tex=.pdf) $@
 
+
 # Figure 2: Reanalysis overview
+figures/event24-reanalysis+nh18fig5.pdf: scripts/plot.py $(event24_ana) $(event24_ens)
+	python3 -m scripts.plot reanalysis+nh18fig5 2016-12-18T00,65,10,35,-30 $(event24_ens) \
+			--reanalysis $(event24_ana) \
+			--output $@ \
+			--hovmoeller-extent "target" \
+			--end 2016-12-24T00
+
+
 # Figure 3: Ensemble target metric overview
-# Figure 6: ESA maps and Hovmöller with full flux
-# Figure 7: Cluster and scatter analysis of upstream region
-figures/event24-%.pdf: scripts/plot.py $(event24_ana) $(event24_ens)
+# Figure 7: ESA maps and Hovmöller with full flux
+# Figure 8: Cluster and scatter analysis of upstream region (needs data from idealized experiment)
+figures/event24-%.pdf: scripts/plot.py figures/idealized.pdf $(event24_ana) $(event24_ens)
 	python3 -m scripts.plot $* 2016-12-18T00,65,10,35,-30 $(event24_ens) \
 			--reanalysis $(event24_ana) \
 			--output $@ \
 			--nsamples $(NSAMPLES) \
 			--sample-size $(SAMPLE_SIZE) \
 			--significance-level $(SIGNIFICANCE_LEVEL) \
+			--ncluster $(NCLUSTER) \
 			--delta-hours "0,-24,-48,-72,-96" \
 			--source "f1+f2+f3" \
 			--scatter 2016-12-16T00,65,-90,35,-15 \
+			--scatter-idealized "data/idealized.json" \
 			--hovmoeller-extent "target" \
-			--end 2016-12-24T00
+			--end 2016-12-21T00
 
 
 # Figure 4: Integrated budget in 2.5 days prior to onset
-figures/event24-budget.pdf: scripts/plot-budget.py $(event24_ana) | figures
+figures/event24-budget.pdf: scripts/plot-budget.py $(event24_ana)
 	python3 -m scripts.plot-budget $(event24_ana) 2016-12-15T12 2016-12-18T00,65,10,35,-30 \
 			--statistics-output "data/budget-statistics.json" \
 			--output $@
 
 
 # Figure 5: Target forecast evaluation
-figures/event24-evaluation.pdf: scripts/plot-evaluation.py $(event24_eval) $(event24_ana) | figures
+figures/event24-evaluation.pdf: scripts/plot-evaluation.py $(event24_eval) $(event24_ana)
 	python3 -m scripts.plot-evaluation 65,10,35,-30 $(event24_eval) \
 		--reanalysis $(event24_ana) \
-		--highlight 2016-12-10T00 2016-12-10T12 2016-12-11T00 \
+		--highlight 2016-12-09T12 2016-12-10T00 2016-12-10T12 2016-12-11T00 \
 		--output $@
 
 
-# Figure 8: ESA maps with flux terms
+# Figure 9: ESA maps with flux terms
 figures/event24-maps-separate.pdf: figures/event24-maps-separate/event24-maps-separate.tex \
 		figures/event24-maps-separate/f1f3.pdf figures/event24-maps-separate/f2.pdf
 	cd $(dir $<) && lualatex $(notdir $<)
@@ -160,6 +197,7 @@ figures/event24-maps-separate/f1f3.pdf: scripts/plot.py $(event24_ana) $(event24
 			--nsamples $(NSAMPLES) \
 			--sample-size $(SAMPLE_SIZE) \
 			--significance-level $(SIGNIFICANCE_LEVEL) \
+			--ncluster $(NCLUSTER) \
 			--delta-hours "24, 0,-24, -48" \
 			--source "f1+f3" \
 			--panel-letters "abcd"
@@ -171,27 +209,31 @@ figures/event24-maps-separate/f2.pdf: scripts/plot.py $(event24_ana) $(event24_e
 			--nsamples $(NSAMPLES) \
 			--sample-size $(SAMPLE_SIZE) \
 			--significance-level $(SIGNIFICANCE_LEVEL) \
+			--ncluster $(NCLUSTER) \
 			--delta-hours "24, 0,-24, -48" \
 			--source "f2" \
 			--panel-letters "efgh"
 
-# Figure 9: Cluster analysis of F2 in block
+
+# Figure 10: Cluster analysis of F2 in block
 figures/event24-cluster-f2.pdf: scripts/plot.py $(event24_ana) $(event24_ens)
 	python3 -m scripts.plot cluster 2016-12-18T00,65,10,35,-30 $(event24_ens) \
 			--reanalysis $(event24_ana) \
 			--output $@ \
 			--nsamples $(NSAMPLES) \
 			--sample-size $(SAMPLE_SIZE) \
 			--significance-level $(SIGNIFICANCE_LEVEL) \
+			--ncluster $(NCLUSTER) \
 			--source "f2"
 
-# Figure 10: Flux-LWA relationship in block (similar to NH18's Fig. 4)
+
+# Figure 11: Flux-LWA relationship in block (similar to NH18's Fig. 4)
 figures/event24-nh18fig4.pdf: scripts/plot.py $(event24_ana) $(event24_ens)
 	python3 -m scripts.plot nh18fig4 2016-12-18T00,65,10,35,-30 $(event24_ens) \
 			--reanalysis $(event24_ana) \
 			--output $@ \
 			--nsamples $(NSAMPLES) \
 			--sample-size $(SAMPLE_SIZE) \
 			--significance-level $(SIGNIFICANCE_LEVEL) \
+			--ncluster $(NCLUSTER) \
 			--scatter 2016-12-18T00,50,-5,40,-15
-

README.md

@@ -1,6 +1,6 @@
 # Block-Dec2016-LWA-ESA
 
-Code to produce all figures of the article **The Onset of a Blocked Flow Event as a "Traffic Jam": Characterization with Ensemble Sensitivity Analysis** (in revision at JAS).
+Code to produce all figures of the article **The Onset of a Blocking Event as a "Traffic Jam": Characterization with Ensemble Sensitivity Analysis** (in revision at JAS).
 
 - Research by [Christopher Polster](https://dynmet.ipa.uni-mainz.de/christopher-polster/) and [Volkmar Wirth](https://dynmet.ipa.uni-mainz.de/volkmar-wirth/).
 - Software by Christopher Polster.
@@ -23,6 +23,10 @@ Use the scripts and instructions in the [`data`](data) directory to obtain the r
 
     data
     + IFS-ENS
+    | + ENS-2016-12-09T12Z
+    | | + ENS-2016-12-09T12Z-t.nc
+    | | + ENS-2016-12-09T12Z-u.nc
+    | | + ENS-2016-12-09T12Z-v.nc
     | + ENS-2016-12-10T00Z
     | | + ENS-2016-12-10T00Z-t.nc
     | | + ENS-2016-12-10T00Z-u.nc
@@ -78,11 +82,12 @@ Figures will appear in the `figures` directory:
 - Fig. 3: `figures/event24-plume.pdf`
 - Fig. 4: `figures/event24-budget.pdf`
 - Fig. 5: `figures/event24-evaluation.pdf`
-- Fig. 6: `figures/event24-maps+hovmoeller.pdf`
-- Fig. 7: `figures/event24-cluster+scatter.pdf`
-- Fig. 8: `figures/event24-maps-separate.pdf`
-- Fig. 9: `figures/event24-cluster-f2.pdf`
-- Fig. 10: `figures/event24-nh18fig4.pdf`, does not contain the fits of [Nakamura and Huang (2018)](https://doi.org/10.1126/science.aat0721) due to licensing restrictions
+- Fig. 6: `figures/idealized.pdf`
+- Fig. 7: `figures/event24-maps+hovmoeller.pdf`
+- Fig. 8: `figures/event24-cluster+scatter.pdf`
+- Fig. 9: `figures/event24-maps-separate.pdf`
+- Fig. 10: `figures/event24-cluster-f2.pdf`
+- Fig. 11: `figures/event24-nh18fig4.pdf`, does not contain the fits of [Nakamura and Huang (2018)](https://doi.org/10.1126/science.aat0721) due to licensing restrictions
 
 Because the statistical significance test is based on a randomized procedure, results may vary slightly each time the plots are created.
 

data/README.md

@@ -4,6 +4,7 @@
 
 IFS data can be downloaded from the ECMWF MARS archive. If you have the required credentials, log in to ecgate and use the `job.sh` script provided in the `IFS-ENS` folder to obtain the files
 
+    $ sbatch job.sh 2016-12-09 12
     $ sbatch job.sh 2016-12-10 00
     $ sbatch job.sh 2016-12-10 12
     $ sbatch job.sh 2016-12-11 00

figures/forecast-combination/forecast-combination.tex

@@ -6,12 +6,12 @@
 
 \begin{document}
     \begin{tikzpicture}[xscale=0.45,yscale=0.5]
-        % No-data zones
-        \fill[black!15!white] (4,0) rectangle (22,4);
+        % Data zone
+        \fill[black!15!white] (4,0) rectangle (20,4);
         % Forecasts
-        \foreach \x in {1,2,3} {
-            \draw[thick,-latex] (\x-3,\x) -- (\x+22,\x);
-            \fill[black] (\x-3,\x) circle (0.15) node[anchor=east] {ENS\x\;};
+        \foreach \x in {1,2,3,4} {
+            \draw[thick,-latex] (\x-4,\x) -- (\x+21,\x);
+            \fill[black] (\x-4,\x) circle (0.15) node[anchor=east] {ENS\x\;};
         }
         % Last initial and target evaluation time
         \draw[thick,dashed] (0,0) -- (0,4);

scripts/common/plotting.py

@@ -64,6 +64,7 @@ def format_date_full(x):
 
 # Axes limits as a multiple of some value
 def mult_limit(values, mult=10, tol=0, minval=None, maxval=None):
+    values = np.asarray(values).flatten()
     # Find next multiple below minium
     ymin = int(min(values) / mult) * mult
     while min(values) < ymin:
@@ -151,7 +152,7 @@ def plot_box(ax, box, **kwargs):
 
 # Plot presets
 
-def plume(ax, x, ys, clusters, reanalysis=None, percentile=None, linewidth=1.0, **kwargs):
+def plume(ax, x, ys, clusters, reanalysis=None, ncluster=None, linewidth=1.0, **kwargs):
     # Unpack clusters
     bot_cluster, not_cluster, top_cluster = clusters
     # Number of members in ensemble to determine line alpha
@@ -164,8 +165,8 @@ def plume(ax, x, ys, clusters, reanalysis=None, percentile=None, linewidth=1.0,
     legend_handles = [line_not[0], line_top[0], line_bot[0]]
     legend_labels = [
         "ENS member",
-        "Top" if percentile is None else "{}% top".format(percentile),
-        "Bottom" if percentile is None else "{}% bottom".format(percentile),
+        "Top" + ("" if ncluster is None else " {}".format(ncluster)),
+        "Bottom" + ("" if ncluster is None else " {}".format(ncluster)),
     ]
     # Add the reanalysis timeseries to the plot and legend
     if reanalysis is not None:

scripts/plot-evaluation.py

@@ -90,13 +90,17 @@ def parse_datehour(arg):
         half = len(fx) // 2
         ax.plot(fx[:half], fy[:half], linewidth=0, color=plotting.esa_blue, marker="x", markersize=4, zorder=10)
         ax.plot(fx[half:], fy[half:], linewidth=0, color=plotting.esa_red, marker="x", markersize=4, zorder=10)
+        # Print some statistics about the inividual forecasts
+        print("INIT {:%Y-%m-%d %H}Z: {:6.2f} max  {:6.2f} min".format(init, max(values), min(values)))
 
     # Add marker and horizontal line for reanalysis value if available
     if args.reanalysis is not None:
         assert init == target_valid
         ax.plot([i], [target_ana], marker="*", markersize=10, color="#000", linewidth=1, zorder=99, label="Reanalysis")
         ax.axhline(target_ana, linewidth=1, color="#000", zorder=0)
         ax.legend(loc="center right", fontsize="small")
+        # Print the reanalysis target value too
+        print("REANALYSIS: {:6.2f}".format(target_ana.values))
 
     # Translate x-ticks to dates and add nice labels
     ticks = np.arange(0, len(target_ens)+1, 4)