From fc611c6de9446b41460c432cdd6a6e53ff98674d Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Martin=20Schr=C3=B6n?= <martin.schroen@ufz.de>
Date: Mon, 27 Nov 2023 16:41:00 +0100
Subject: [PATCH] Feature: muons processing and plotting
---
VERSION | 2 +-
corny/version.py | 2 +-
default.cfg | 14 ++++++++
instantPASDy.py | 88 ++++++++++++++++++++++++++++++++++++++++++++++++
4 files changed, 104 insertions(+), 2 deletions(-)
diff --git a/VERSION b/VERSION
index 972ef76..100435b 100644
--- a/VERSION
+++ b/VERSION
@@ -1 +1 @@
-0.7.9
+0.8.2
diff --git a/corny/version.py b/corny/version.py
index 07b67b5..40c8236 100644
--- a/corny/version.py
+++ b/corny/version.py
@@ -1,6 +1,6 @@
# Corny module: version.py
# Version management
-VERSION = "0.8.1"
+VERSION = "0.8.2"
from .io import app_file_path
diff --git a/default.cfg b/default.cfg
index c1518bd..d3812e0 100644
--- a/default.cfg
+++ b/default.cfg
@@ -122,6 +122,10 @@ latitude_column =
longitude_column =
altitude_column =
+## Muons
+# Column for muon counts
+muon_column =
+
[clean]
# Cut the data, find outliers, fill gaps, and smooth out noise.
@@ -280,6 +284,10 @@ temperature_range = -60, 60
# Units: m³/m³
sm_range =
+## Muon count range
+# Units: cph
+muon_range =
+
## Split tracks
# Identify single tracks between locations A and B (yes/no)
@@ -373,6 +381,12 @@ NM_station = JUNG
# - 1day # 1 day temporal resolution
NM_resolution = 1min
+## Muon counts
+# Muon correction parameters based on Stevanato et al. (2022)
+muon_beta = 625
+muon_alpha = 0.0021
+muon_T_ref = 0
+
[conversion]
# New column name for the soil moisture product
diff --git a/instantPASDy.py b/instantPASDy.py
index 0c65782..f795e90 100644
--- a/instantPASDy.py
+++ b/instantPASDy.py
@@ -575,6 +575,8 @@ def main(configfile=None, cfg_prefix=None):
(Ta,Tb) = xsplit(config['clean']['temperature_range'], range=2, type=float)
if 'timeres_range' in config['clean']:
(tra,trb)= xsplit(config['clean']['timeres_range'], range=2, type=float)
+
+
# Define new column N+clean
data[N+clean] = data[N]
@@ -710,6 +712,31 @@ def main(configfile=None, cfg_prefix=None):
#data[['GpsUTC', 'LatDec', 'LongDec', 'Alt']].to_csv('testgps.csv', index=False)
+ ####################################################
+ ###### Muons #######################################
+ ####################################################
+
+ Muons = False
+ if cc(config,'input','muon_column'):
+ Muons = True
+ mu = gc(config,'input','muon_column')
+ (Ma,Mb) = gc(config,'clean','muon_range', range=2, dtype=list, types=float)
+ data[mu+"_cph"] = Unit_convert(data[mu], data[tres], 'cph')
+ data[mu+clean] = find_outlier(data, mu+"_cph", lower=Ma, upper=Mb)
+
+ beta_M = gc(config,'correction','muon_beta', dtype=float)
+ alpha_M = gc(config,'correction','muon_alpha', dtype=float)
+ Tref_M = gc(config,'correction','muon_T_ref', dtype=float)
+
+ data[mu+clean+"_p"] = data[mu+clean] * np.exp((data[p] - 1013) / beta_M)
+ data[mu+clean+"_ph"] = data[mu+clean+"_p"] * (1 + alpha_M * (data[T] - Tref_M))
+
+ mufinal = mu+clean+"_ph"
+
+ J.report("muons",
+ "Muons have been dropped beyond below {:.0f} and above {:.0f} cph. " \
+ +"Corrections were applied based on Stevanato et al. (2022) using beta={:.0f}, alpha={:.0f}, and T_ref={:.0f}.",
+ Ma, Mb, beta_M, alpha_M, Tref_M)
####################################################
###### Soil #######################################
@@ -1638,6 +1665,13 @@ def main(configfile=None, cfg_prefix=None):
J.report("neutronsth",
"Temporal smoothing was applied over {:.0f} time steps.",
smooth_window)
+
+ if Muons:
+ data[mufinal+mov] = data[mufinal].rolling(smooth_window, center=True).mean()
+ mufinal += mov
+ J.report("muons",
+ "Temporal smoothing was applied over {:.0f} time steps.",
+ smooth_window)
if cc(config,'clean','smooth_radius'):
@@ -2791,6 +2825,60 @@ def main(configfile=None, cfg_prefix=None):
buffer = F.buff
R.add_image(buffer)
+ if Muons:
+ Pr.update("Muons ")
+ R.add_page()
+ R.add_title("Muons")
+ R.add_par(J.post("muons"))
+
+ if not data[mu].isnull().all():
+ with Figure(size=(10,11), layout=(3,1),
+ **kw_fig_buff, **kw_fig_time,
+ x_ticks_last=True) as F:
+
+ ax = F.axes[0]
+ # ax.scatter(data_orig.index, data_orig[N],
+ # label="$N$ raw",
+ # marker='o', s=5, lw=0, color='#DDDDDD')
+ ax.scatter(data.index, data[mu+"_cph"],
+ label="Muons raw",
+ marker='o', s=10, lw=0, color='#CCCCCC')
+ data_outliers = data[data[mu+clean] != data[mu+"_cph"]]
+ ax.scatter(data_outliers.index, data_outliers[mu+"_cph"],
+ label="Muons raw (dropped)",
+ marker='+', s=40, lw=2, color='red')
+
+ ax = F.axes[1]
+ ax.scatter(data.index, data[mu+clean],
+ label="Muons cleaned",
+ marker='o', s=10, lw=0, color='#CCCCCC')
+ ax.scatter(data.index, data[mu+clean+"_ph"],
+ label="$Muons cleaned, atm. corrections",
+ marker='o', s=10, lw=0, color='#999999')
+
+ ax = F.axes[2]
+ ax.plot(data.index, (data[mufinal] - data[mufinal].mean())/data[mufinal].std(ddof=0),
+ label="Muons, normalized",
+ lw=1, color='C0')
+ ax.plot(data.index, (data["NM"] - data["NM"].mean())/data["NM"].std(ddof=0),
+ label="Neutron Monitor, normalized",
+ lw=1, color='C3')
+ ax.set_ylabel("z score (--)")
+
+ for ax in F.axes[0:2]:
+ ax.set_ylabel("Muons (cph)")
+ for ax in F.axes:
+ ax.set_xlim(date1, date2)
+ ax.grid(color="k", alpha=0.1)
+ ax.legend(loc="upper left", bbox_to_anchor=(0, 1.12),
+ ncol=3, fancybox=False, framealpha=1, edgecolor="k",
+ shadow=True)
+ if ax != F.axes[-1]:
+ ax.set_xticklabels([])
+
+ buffer = F.buff
+ R.add_image(buffer)
+
Pr.update("Atmosphere")
R.add_page()
R.add_title("Atmosphere")
--
GitLab