GC-Net-level-1-data-processing/jaws_tools.py at main · GEUS-Glaciology-and-Climate/GC-Net-level-1-data-processing · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
# -*- coding: utf-8 -*-
"""
Functions adapted from JAWS:
    https://github.com/jaws/jaws/blob/master/jaws/gcnet2nc.py

"""
import numpy as np
import pandas as pd
import xarray as xr
import sys

try:
    from metpy.units import units
    from metpy.calc import potential_temperature, density
    from metpy.calc import (
        mixing_ratio_from_relative_humidity,
        specific_humidity_from_mixing_ratio,
    )
except ImportError:
    print("ImportError: No module named metpy.units")
    print(
        "HINT: This is a known issue for a dependent package that occurs using pip installation in Python 2.7 \n"
        "To fix it, please perform either of below operations: \n"
        "1. Install jaws using conda as: conda install -c conda-forge jaws \n"
        "2. Upgrade to Python version >= 3.6 \n"
    )
    print(
        "If none of the above works for you, "
        "please inform the JAWS maintainers by opening an issue at https://github.com/jaws/jaws/issues."
    )

    sys.exit(1)

try:
    from jaws import common, sunposition, clearsky, tilt_angle, fsds_adjust
except ImportError:
    import common, sunposition, clearsky, tilt_angle, fsds_adjust


def gradient_fluxes(df):  # This method is very sensitive to input data quality
    """Returns Sensible Heat Flux and Latent Heat Flux based on Steffen & DeMaria (1996) method"""
    g = 9.81  # m/s**2
    cp = 1005  # J/kg/K
    k = 0.4  # von Karman
    Lv = 2.50e6  # J/kg

    # Average temp from both sensors for height1 and height2
    ta1 = df.loc[:, ("TA1", "TA3")]
    ta2 = df.loc[:, ("TA2", "TA4")]
    ta1 = ta1.mean(axis=1)
    ta2 = ta2.mean(axis=1)

    ht_low  = np.array(df["HW1"], dtype=float, copy=True)
    ht_high = np.array(df["HW2"], dtype=float, copy=True)
    ta_low  = np.array(ta1, dtype=float, copy=True)
    ta_high = np.array(ta2, dtype=float, copy=True)
    VW_low  = np.array(df["VW1"], dtype=float, copy=True)
    VW_high = np.array(df["VW2"], dtype=float, copy=True)
    RH_low  = np.array(df["RH1"], dtype=float, copy=True)
    RH_high = np.array(df["RH2"], dtype=float, copy=True)
    P       = np.array(df["P"], dtype=float, copy=True)

    hw1 = np.array(df["HW1"], dtype=float, copy=True)
    hw2 = np.array(df["HW2"], dtype=float, copy=True)
    vw1 = np.array(df["VW1"], dtype=float, copy=True)
    vw2 = np.array(df["VW2"], dtype=float, copy=True)
    rh1 = np.array(df["RH1"], dtype=float, copy=True)
    rh2 = np.array(df["RH2"], dtype=float, copy=True)
    ta1a = np.array(ta1, dtype=float, copy=True)
    ta2a = np.array(ta2, dtype=float, copy=True)

    ind = ht_high < ht_low

    ht_low[ind]  = hw2[ind]
    ht_high[ind] = hw1[ind]
    ta_low[ind]  = ta2a[ind]
    ta_high[ind] = ta1a[ind]
    VW_low[ind]  = vw2[ind]
    VW_high[ind] = vw1[ind]
    RH_low[ind]  = rh2[ind]
    RH_high[ind] = rh1[ind]

    # Potential Temperature
    pot_tmp_low = potential_temperature(
        units.Quantity(P, "hPa"), units.Quantity(ta_low, "degC")
    ).magnitude
    pot_tmp_high = potential_temperature(
        units.Quantity(P, "hPa"), units.Quantity(ta_high, "degC")
    ).magnitude
    pot_tmp_avg = (pot_tmp_low + pot_tmp_high) / 2
    ta_avg = (ta_low + ta_high) / 2

    # Ri
    du = VW_high - VW_low
    du[du == 0] = np.nan
    pot_tmp_avg[pot_tmp_avg == 0] = np.nan
    ri = g * (pot_tmp_high - pot_tmp_low) * (ht_high - ht_low) / (pot_tmp_avg * du)

    # Phi
    phi_m = (1 - 5.2 * ri) ** -1
    phi_h = phi_m
    phi_m[ri < -0.03] = (1 - 18 * ri[ri < -0.03]) ** -0.25
    phi_h[ri < -0.03] = phi_m[ri < -0.03] / 1.3
    phi_m[(-0.03 <= ri) & (ri < 0)] = (1 - 18 * ri[(-0.03 <= ri) & (ri < 0)]) ** -0.25
    phi_h[(-0.03 <= ri) & (ri < 0)] = phi_m[(-0.03 <= ri) & (ri < 0)]
    phi_e = phi_h

    # air density
    RHo = density(
        P * units.hectopascal, ta_avg * units.degC, 0
    ).magnitude  # Use average temperature

    # SH
    ht_low[ht_low == 0] = np.nan
    num = -RHo * cp * k ** 2 * (pot_tmp_high - pot_tmp_low) * (VW_high - VW_low)

    dnm = phi_h * phi_m * np.log(ht_high / ht_low) ** 2
    dnm[dnm == 0] = np.nan
    sh = num / dnm
    sh[np.abs(sh) >= 100] = np.nan

    # Specific Humidity
    mixing_ratio_low = mixing_ratio_from_relative_humidity(
        units.Quantity(P, "hPa"), units.Quantity(ta_low, "degC"), RH_low
    )
    mixing_ratio_high = mixing_ratio_from_relative_humidity(
        units.Quantity(P, "hPa"), units.Quantity(ta_high, "degC"), RH_high
    )
    q_low = specific_humidity_from_mixing_ratio(mixing_ratio_low).magnitude
    q_high = specific_humidity_from_mixing_ratio(mixing_ratio_high).magnitude
    q_low = q_low / 100  # Divide by 100 to make it in range [0,1]
    q_high = q_high / 100

    # LH
    num = -RHo * Lv * k ** 2 * (q_high - q_low) * (VW_high - VW_low)
    dnm = phi_e * phi_m * np.log(ht_high / ht_low) ** 2
    dnm[dnm == 0] = np.nan
    lh = num / dnm
    lh[np.abs(lh) >= 100] = np.nan
    return -sh, -lh