[WIP] Implement Rec. ITU-R P.2001 propagation model

pycraf/itudata/p.2001-3/make_npy.py

@@ -0,0 +1,112 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+import numpy as np
+import os
+
+
+BASEPATH = 'R-REC-P.2001-3-201908'
+
+
+def main():
+
+    lons, lats = np.meshgrid(
+        np.linspace(0, 360, 241), np.linspace(90, -90, 121),
+        )
+
+    dn_median = np.genfromtxt(os.path.join(BASEPATH, 'DN_Median.txt'))
+    dn_supslope = np.genfromtxt(os.path.join(BASEPATH, 'DN_SupSlope.txt'))
+    dn_subslope = np.genfromtxt(os.path.join(BASEPATH, 'DN_SubSlope.txt'))
+    dn_dz = np.genfromtxt(os.path.join(BASEPATH, 'dndz_01.txt'))
+
+    surfwv_50 = np.genfromtxt(os.path.join(BASEPATH, 'surfwv_50_fixed.txt'))
+
+    foes_50 = np.genfromtxt(os.path.join(BASEPATH, 'FoEs50.txt'))
+    foes_10 = np.genfromtxt(os.path.join(BASEPATH, 'FoEs10.txt'))
+    foes_1 = np.genfromtxt(os.path.join(BASEPATH, 'FoEs01.txt'))
+    foes_01 = np.genfromtxt(os.path.join(BASEPATH, 'FoEs0.1.txt'))
+
+    h0 = np.genfromtxt(os.path.join(BASEPATH, 'h0.txt'))
+
+    # import matplotlib.pyplot as plt
+
+    # plt.contourf(lons, lats, dn_dz, 128)
+    # plt.show()
+
+    save_kwargs = {
+        'lons': lons.astype(np.float32),
+        'lats': lats.astype(np.float32),
+        'dn_median': dn_median.astype(np.float32),
+        'dn_supslope': dn_supslope.astype(np.float32),
+        'dn_subslope': dn_subslope.astype(np.float32),
+        'dn_dz': dn_dz.astype(np.float32),
+        }
+    np.savez('refract_map', **save_kwargs)
+
+    save_kwargs = {
+        'lons': lons.astype(np.float64),
+        'lats': lats.astype(np.float64),
+        'surfwv_50': surfwv_50.astype(np.float64),
+        }
+    np.savez('wv_map', **save_kwargs)
+
+    save_kwargs = {
+        'lons': lons.astype(np.float32),
+        'lats': lats.astype(np.float32),
+        'h0': h0.astype(np.float32),
+        }
+    np.savez('h0_map', **save_kwargs)
+
+    save_kwargs = {
+        'lons': lons.astype(np.float32),
+        'lats': lats.astype(np.float32),
+        'foes_50': foes_50.astype(np.float32),
+        'foes_10': foes_10.astype(np.float32),
+        'foes_1': foes_1.astype(np.float32),
+        'foes_01': foes_01.astype(np.float32),
+        }
+    np.savez('sporadic_e_map', **save_kwargs)
+
+    # read using:
+    # dat = np.load('refract_map.npz')
+
+    lons, lats = np.meshgrid(
+        np.linspace(0, 360, 321), np.linspace(90, -90, 161),
+        )
+
+    pr6 = np.genfromtxt(os.path.join(BASEPATH, 'Esarain_Pr6_v5.txt'))
+    mt = np.genfromtxt(os.path.join(BASEPATH, 'Esarain_Mt_v5.txt'))
+    beta = np.genfromtxt(os.path.join(BASEPATH, 'Esarain_Beta_v5.txt'))
+
+    # import matplotlib.pyplot as plt
+
+    # plt.contourf(lons, lats, pr6, 128)
+    # plt.show()
+
+    save_kwargs = {
+        'lons': lons.astype(np.float64),
+        'lats': lats.astype(np.float64),
+        'pr6': pr6.astype(np.float64),
+        'mt': mt.astype(np.float64),
+        'beta': beta.astype(np.float64),
+        }
+    np.savez('rain_map', **save_kwargs)
+
+    tropoclim = np.genfromtxt(
+        os.path.join(BASEPATH, 'TropoClim.txt'), dtype=np.int8
+        )
+
+    lons, lats = np.meshgrid(
+        np.linspace(-179.75, 179.75, 720), np.linspace(89.75, -89.75, 360),
+        )
+
+    save_kwargs = {
+        'lons': lons.astype(np.float32),
+        'lats': lats.astype(np.float32),
+        'tropoclim': tropoclim,  # use "nearest neighbor" for interpolation!
+        }
+    np.savez('tropoclim_map', **save_kwargs)
+
+
+if __name__ == '__main__':
+    main()

pycraf/itudata/p.2001-3/table_c.2.1.txt

@@ -0,0 +1,50 @@
+# index rel_height prob
+1 -2400 0.000555
+2 -2300 0.000802
+3 -2200 0.001139
+4 -2100 0.001594
+5 -2000 0.002196
+6 -1900 0.002978
+7 -1800 0.003976
+8 -1700 0.005227
+9 -1600 0.006764
+10 -1500 0.008617
+11 -1400 0.010808
+12 -1300 0.013346
+13 -1200 0.016225
+14 -1100 0.019419
+15 -1000 0.022881
+16 -900 0.026542
+17 -800 0.030312
+18 -700 0.034081
+19 -600 0.037724
+20 -500 0.041110
+21 -400 0.044104
+22 -300 0.046583
+23 -200 0.048439
+24 -100 0.049589
+25 0 0.049978
+26 100 0.049589
+27 200 0.048439
+28 300 0.046583
+29 400 0.044104
+30 500 0.041110
+31 600 0.037724
+32 700 0.034081
+33 800 0.030312
+34 900 0.026542
+35 1000 0.022881
+36 1100 0.019419
+37 1200 0.016225
+38 1300 0.013346
+39 1400 0.010808
+40 1500 0.008617
+41 1600 0.006764
+42 1700 0.005227
+43 1800 0.003976
+44 1900 0.002978
+45 2000 0.002196
+46 2100 0.001594
+47 2200 0.001139
+48 2300 0.000802
+49 2400 0.000555

pycraf/pathprof/helper.py

@@ -91,6 +91,103 @@ def __init__(self, *args, **kwargs):
     )
 
 
+_refract_data_p2001 = np.load(get_pkg_data_filename(
+    '../itudata/p.2001-3/refract_map.npz'
+    ))
+
+_dn_median_interpolator = RegularGridInterpolator(
+    (_refract_data_p2001['lons'][0], _refract_data_p2001['lats'][::-1, 0]),
+    _refract_data_p2001['dn_median'][::-1].T
+    )
+_dn_supslope_interpolator = RegularGridInterpolator(
+    (_refract_data_p2001['lons'][0], _refract_data_p2001['lats'][::-1, 0]),
+    _refract_data_p2001['dn_supslope'][::-1].T
+    )
+_dn_subslope_interpolator = RegularGridInterpolator(
+    (_refract_data_p2001['lons'][0], _refract_data_p2001['lats'][::-1, 0]),
+    _refract_data_p2001['dn_subslope'][::-1].T
+    )
+_dn_dz_interpolator = RegularGridInterpolator(
+    (_refract_data_p2001['lons'][0], _refract_data_p2001['lats'][::-1, 0]),
+    _refract_data_p2001['dn_dz'][::-1].T
+    )
+
+_wv_data_p2001 = np.load(get_pkg_data_filename(
+    '../itudata/p.2001-3/wv_map.npz'
+    ))
+
+_surfwv_50_interpolator = RegularGridInterpolator(
+    (_wv_data_p2001['lons'][0], _wv_data_p2001['lats'][::-1, 0]),
+    _wv_data_p2001['surfwv_50'][::-1].T
+    )
+
+_h0_data_p2001 = np.load(get_pkg_data_filename(
+    '../itudata/p.2001-3/h0_map.npz'
+    ))
+
+_h0_interpolator = RegularGridInterpolator(
+    (_h0_data_p2001['lons'][0], _h0_data_p2001['lats'][::-1, 0]),
+    _h0_data_p2001['h0'][::-1].T
+    )
+
+_spo_e_data_p2001 = np.load(get_pkg_data_filename(
+    '../itudata/p.2001-3/sporadic_e_map.npz'
+    ))
+
+_foes_50_interpolator = RegularGridInterpolator(
+    (_spo_e_data_p2001['lons'][0], _spo_e_data_p2001['lats'][::-1, 0]),
+    _spo_e_data_p2001['foes_50'][::-1].T
+    )
+_foes_10_interpolator = RegularGridInterpolator(
+    (_spo_e_data_p2001['lons'][0], _spo_e_data_p2001['lats'][::-1, 0]),
+    _spo_e_data_p2001['foes_10'][::-1].T
+    )
+_foes_1_interpolator = RegularGridInterpolator(
+    (_spo_e_data_p2001['lons'][0], _spo_e_data_p2001['lats'][::-1, 0]),
+    _spo_e_data_p2001['foes_1'][::-1].T
+    )
+_foes_01_interpolator = RegularGridInterpolator(
+    (_spo_e_data_p2001['lons'][0], _spo_e_data_p2001['lats'][::-1, 0]),
+    _spo_e_data_p2001['foes_01'][::-1].T
+    )
+
+_rain_data_p2001 = np.load(get_pkg_data_filename(
+    '../itudata/p.2001-3/rain_map.npz'
+    ))
+
+_pr6_interpolator = RegularGridInterpolator(
+    (_rain_data_p2001['lons'][0], _rain_data_p2001['lats'][::-1, 0]),
+    _rain_data_p2001['pr6'][::-1].T
+    )
+_mt_interpolator = RegularGridInterpolator(
+    (_rain_data_p2001['lons'][0], _rain_data_p2001['lats'][::-1, 0]),
+    _rain_data_p2001['mt'][::-1].T
+    )
+_beta_interpolator = RegularGridInterpolator(
+    (_rain_data_p2001['lons'][0], _rain_data_p2001['lats'][::-1, 0]),
+    _rain_data_p2001['beta'][::-1].T
+    )
+
+_tclim_data_p2001 = np.load(get_pkg_data_filename(
+    '../itudata/p.2001-3/tropoclim_map.npz'
+    ))
+
+# BEWARE: UNLIKE FOR THE OTHER INTERPOLATORS, LONS ARE IN [-180, 180]
+_tropoclim_interpolator = RegularGridInterpolator(
+    (_tclim_data_p2001['lons'][0], _tclim_data_p2001['lats'][::-1, 0]),
+    _tclim_data_p2001['tropoclim'][::-1].T, method='nearest',
+    )
+
+
+_rain_probs = np.genfromtxt(
+    get_pkg_data_filename(
+        '../itudata/p.2001-3/table_c.2.1.txt'
+        ),
+    dtype=(np.int8, np.float64, np.float64),
+    names=True,
+    )
+
+
 @utils.ranged_quantity_input(
     p_w=(0, 100, apu.percent),
     phi=(-90, 90, apu.deg),
@@ -150,6 +247,47 @@ def _DN_N0_from_map(lon, lat):
     return _DN, _N0
 
 
+def _DN_P2001_from_map(lon, lat):
+
+    dn_median = _dn_median_interpolator((lon % 360, lat))
+    dn_supslope = _dn_supslope_interpolator((lon % 360, lat))
+    dn_subslope = _dn_subslope_interpolator((lon % 360, lat))
+    dn_dz = _dn_dz_interpolator((lon % 360, lat))
+
+    return dn_median, dn_supslope, dn_subslope, dn_dz
+
+
+def _sporadic_E_P2001_from_map(lon, lat, p):
+
+    lon, lat, p = np.broadcast_arrays(lon, lat, p)
+    f_oes1 = np.empty(lon.shape, dtype=np.float64)
+    f_oes2 = np.empty(lon.shape, dtype=np.float64)
+    p1 = np.empty(lon.shape, dtype=np.float64)
+    p2 = np.empty(lon.shape, dtype=np.float64)
+
+    mask01 = p < 0.01
+    mask10 = p > 0.1
+    mask_m = (~mask01) & (~mask10)  # 1% <= p <= 10%
+    f_oes1[mask01] = _foes_01_interpolator((lon % 360, lat))
+    f_oes2[mask01] = _foes_1_interpolator((lon % 360, lat))
+    p1[mask01] = 0.001
+    p2[mask01] = 0.01
+
+    f_oes1[mask_m] = _foes_1_interpolator((lon % 360, lat))
+    f_oes2[mask_m] = _foes_10_interpolator((lon % 360, lat))
+    p1[mask_m] = 0.01
+    p2[mask_m] = 0.1
+
+    f_oes1[mask10] = _foes_10_interpolator((lon % 360, lat))
+    f_oes2[mask10] = _foes_50_interpolator((lon % 360, lat))
+    p1[mask10] = 0.1
+    p2[mask10] = 0.5
+
+    f_oes = f_oes1 + (f_oes2 - f_oes1) * np.log10(p / p1) / np.log10(p2 / p1)
+
+    return f_oes
+
+
 @utils.ranged_quantity_input(
     lon=(-180, 360, apu.deg),
     lat=(-90, 90, apu.deg),