
import numpy as np
import datetime as dt
import pandas as pd
import xarray as xr
import matplotlib as mpl  # Needed to access mpl.colormaps[cmapname]
import matplotlib.pyplot as plt
import cartopy.crs as ccrs


# Read in tracks for a year of LPT.
fn_lpt_systems = ('data/imerg/g50_72h/thresh12/systems/'
                  + 'lpt_systems_imerg_2022060100_2023063023.nc')
lpt_systems = xr.load_dataset(fn_lpt_systems)
lpt_systems

<xarray.Dataset>
Dimensions:                               (nlpt: 108, nstitch: 51707, nobj: 1615)
Coordinates:
  * nlpt                                  (nlpt) int32 0 1 2 3 ... 105 106 107
  * nstitch                               (nstitch) int32 0 1 2 ... 51705 51706
  * nobj                                  (nobj) int32 0 1 2 ... 1612 1613 1614
Data variables: (12/38)
    lptid                                 (nlpt) float64 0.1 1.1 ... 60.2 61.1
    lpt_begin_index                       (nlpt) int64 1 177 598 ... 51295 51526
    lpt_end_index                         (nlpt) int64 175 596 ... 51524 51705
    centroid_lon_start                    (nlpt) float64 264.9 298.4 ... 204.1
    centroid_lat_start                    (nlpt) float64 13.69 -12.56 ... -23.2
    centroid_lon_end                      (nlpt) float64 290.1 261.9 ... 206.4
    ...                                    ...
    min_filtered_running_field            (nstitch) float64 nan 12.0 ... nan
    max_inst_field                        (nstitch) float64 nan 40.38 ... nan
    max_running_field                     (nstitch) float64 nan 196.9 ... nan
    max_filtered_running_field            (nstitch) float64 nan 18.93 ... nan
    num_objects                           (nlpt) int32 175 420 723 ... 230 180
    objid                                 (nlpt, nobj) float64 2.022e+13 ... nan
Attributes:
    description:  LPT Systems NetCDF file.

array([  0,   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], dtype=int32)

array([    0,     1,     2, ..., 51704, 51705, 51706], dtype=int32)

array([   0,    1,    2, ..., 1612, 1613, 1614], dtype=int32)

array([ 0.1 ,  1.1 ,  2.1 ,  2.2 ,  2.3 ,  2.4 ,  2.5 ,  2.6 ,  2.7 ,
        2.8 ,  2.9 ,  3.1 ,  4.1 ,  4.2 ,  5.1 ,  6.1 ,  7.1 ,  8.1 ,
        8.2 ,  9.1 ,  9.2 , 10.1 , 11.1 , 12.1 , 13.1 , 14.1 , 15.1 ,
       16.1 , 17.1 , 18.1 , 19.1 , 20.1 , 21.1 , 21.2 , 21.3 , 21.4 ,
       21.5 , 21.6 , 21.7 , 21.8 , 21.9 , 22.1 , 23.1 , 24.1 , 24.2 ,
       25.1 , 26.1 , 27.1 , 28.1 , 29.1 , 29.2 , 29.3 , 29.4 , 30.1 ,
       31.1 , 32.1 , 33.1 , 34.1 , 35.1 , 36.1 , 37.01, 37.02, 37.03,
       37.04, 37.05, 37.06, 37.07, 37.08, 37.09, 37.1 , 37.11, 37.12,
       37.13, 37.14, 37.15, 37.16, 37.17, 37.18, 37.19, 38.1 , 38.2 ,
       39.1 , 40.1 , 41.1 , 42.1 , 43.1 , 44.1 , 44.2 , 44.3 , 45.1 ,
       46.1 , 47.1 , 48.1 , 49.1 , 49.2 , 50.1 , 51.1 , 52.1 , 53.1 ,
       54.1 , 55.1 , 56.1 , 57.1 , 58.1 , 59.1 , 60.1 , 60.2 , 61.1 ])

array([    1,   177,   598,  1166,  1822,  2407,  2867,  3415,  3892,
        4319,  4834,  5278,  5457,  6236,  6884,  7016,  7220,  7474,
        7666,  7940,  8301,  8710,  8908,  9081,  9346,  9536,  9716,
       10086, 10281, 10586, 10800, 11029, 11199, 11838, 12995, 13426,
       14375, 14903, 15172, 16223, 17792, 18061, 18367, 18576, 18821,
       19006, 19206, 19380, 19639, 20012, 20969, 22064, 22552, 23149,
       23339, 23553, 24393, 24814, 25129, 25361, 25831, 26166, 27417,
       28811, 29689, 30425, 31520, 32758, 33480, 34060, 34913, 35909,
       36389, 36727, 37608, 38632, 39140, 39506, 39822, 40281, 40787,
       40917, 41087, 41491, 41930, 42657, 43091, 43934, 44546, 45436,
       45607, 45950, 46198, 46788, 47233, 48122, 48444, 48653, 48921,
       49094, 49293, 49674, 50006, 50420, 50840, 51093, 51295, 51526])

array([  175,   596,  1164,  1820,  2405,  2865,  3413,  3890,  4317,
        4832,  5276,  5455,  6234,  6882,  7014,  7218,  7472,  7664,
        7938,  8299,  8708,  8906,  9079,  9344,  9534,  9714, 10084,
       10279, 10584, 10798, 11027, 11197, 11836, 12993, 13424, 14373,
       14901, 15170, 16221, 17790, 18059, 18365, 18574, 18819, 19004,
       19204, 19378, 19637, 20010, 20967, 22062, 22550, 23147, 23337,
       23551, 24391, 24812, 25127, 25359, 25829, 26164, 27415, 28809,
       29687, 30423, 31518, 32756, 33478, 34058, 34911, 35907, 36387,
       36725, 37606, 38630, 39138, 39504, 39820, 40279, 40785, 40915,
       41085, 41489, 41928, 42655, 43089, 43932, 44544, 45434, 45605,
       45948, 46196, 46786, 47231, 48120, 48442, 48651, 48919, 49092,
       49291, 49672, 50004, 50418, 50838, 51091, 51293, 51524, 51705])

array([264.87798633, 298.44829434, 133.54631309, 133.54631309,
       133.54631309,  85.91442916,  85.91442916,  85.91442916,
        63.41018398,  63.41018398,  63.41018398, 144.47854259,
        66.93761735, 142.74898238, 343.29780996, 272.49203179,
       291.86942424,  97.12735815, 137.29952693, 131.43546763,
        87.24813278, 160.00509571, 357.38203802, 261.05124077,
       132.21337565, 153.32930031, 113.9831724 , 298.38842988,
       264.86111209, 142.34445192, 305.1474526 ,  97.35626469,
       164.9026462 , 164.9026462 ,  93.56717789,  93.56717789,
       160.1898706 , 160.1898706 , 160.1898706 , 160.1898706 ,
       146.54337991, 286.71217638, 299.74176382, 306.41704515,
       286.98920271, 169.31832789, 219.51087736, 317.11635191,
       139.89183824, 114.8487026 , 114.8487026 , 114.8487026 ,
       114.8487026 , 304.36425844, 181.62523223, 305.59611129,
        27.93553668,  30.08748178, 216.36390185,  76.18666844,
       217.21045985, 145.49801374, 145.49801374, 145.49801374,
       145.49801374,  91.77879909,  91.77879909,  91.77879909,
        91.77879909, 202.77236803, 202.77236803, 202.77236803,
       202.77236803, 132.05140259, 132.05140259, 132.05140259,
       132.05140259, 191.87676388, 191.87676388, 318.15863692,
       318.15863692,  35.92924971, 292.43058189,  20.68682166,
       289.99423818,  41.1674298 ,  80.65962064,  80.65962064,
        80.65962064, 299.59199174, 294.1824324 ,   3.70651139,
        93.14198003, 154.36378926, 154.36378926, 205.24205012,
       123.69055503, 278.02056305, 290.34103108, 252.93249173,
        62.8686816 ,  89.4285249 , 138.57687965, 137.14376774,
       249.58336291, 111.37922297,  91.0178811 , 204.11091005])

array([ 13.68680658, -12.55953742,  -5.3026228 ,  -5.3026228 ,
        -5.3026228 ,  26.31486253,  26.31486253,  26.31486253,
         4.02750678,   4.02750678,   4.02750678, -22.09959478,
         2.71349452,  15.60333373,  10.17883172,   8.21850129,
         4.06212601,  15.59037854,   6.51385677,  14.51797171,
        -5.13582228,  22.54046838,  10.51538162,  11.56816525,
        20.62573742,  21.31750349,   1.43382765,  18.59131805,
        12.54004048,  20.22359661,  15.59161805,  -4.65019193,
         7.41392832,   7.41392832,   5.43949372,   5.43949372,
        13.23252293,  13.23252293,  13.23252293,  13.23252293,
        -1.00389433,  10.85520565,   8.81394351,  -9.39592149,
         3.67266405, -17.34906463, -19.58929795, -23.04381141,
       -13.35504259,   4.00657236,   4.00657236,   4.00657236,
         4.00657236,  27.178941  , -16.90566192, -13.42146215,
       -18.62442787, -18.70145262, -28.22281533,  -5.74273831,
       -30.85405915, -17.36389696, -17.36389696, -17.36389696,
       -17.36389696,  -5.06959774,  -5.06959774,  -5.06959774,
        -5.06959774, -25.31567829, -25.31567829, -25.31567829,
       -25.31567829, -21.02265887, -21.02265887, -21.02265887,
       -21.02265887, -12.98275541, -12.98275541, -22.15255864,
       -22.15255864, -12.54800793,  -5.83487848,  -4.15318361,
        -3.39001687, -19.46304076,  -5.66381583,  -5.66381583,
        -5.66381583, -22.67074187,  -3.14498111,   1.5604668 ,
        -1.49550722,   5.1667175 ,   5.1667175 , -17.87890786,
        34.10672115,   3.91485012,  23.14140291,   6.34554849,
         3.51376521,  10.03810035,   9.55456181,   1.87494723,
         8.58107359,  26.58626218,  20.83401845, -23.20271191])

array([290.11375956, 261.86379743, 164.01439161, 138.6722078 ,
        80.06594906, 164.01439161, 138.6722078 ,  80.06594906,
       164.01439161, 138.6722078 ,  80.06594906, 166.1449985 ,
        69.679826  ,  69.679826  , 335.49414543, 291.13106995,
       276.89695131,  90.6268527 ,  90.6268527 ,  96.55115465,
        96.55115465, 149.67817931, 353.39919368, 248.64958715,
       132.52708816, 128.6594842 ,  85.76963832, 310.68826849,
       261.54525087, 138.06665509, 297.51275771,  85.81538431,
       121.15634319, 121.27829238, 121.15634319, 121.27829238,
        94.04376021, 163.1983501 , 121.15634319, 121.27829238,
       121.27829238, 301.97472429, 268.74045578, 298.98689814,
       298.98689814, 188.69132818, 227.55234859, 303.29382759,
       134.08690148, 168.33548766, 125.31226093,  91.39157634,
       154.70482054, 335.22750566, 188.81300786, 291.19791247,
        20.46962725,  27.16128653, 189.44249049,  47.00376743,
       224.18966745, 111.69361299, 224.18966745, 165.9954038 ,
       101.71185027, 111.69361299, 224.18966745, 165.9954038 ,
       101.71185027, 111.69361299, 224.18966745, 165.9954038 ,
       101.71185027, 111.69361299, 224.18966745, 165.9954038 ,
       101.71185027, 111.69361299, 224.18966745, 296.31463644,
       322.08719331,  33.02979746, 293.22788169,  78.82917509,
       296.88016276,  28.38888061, 229.60342651, 126.86305395,
       166.93092903, 305.73335906, 306.94247067,   2.14528916,
        73.01983164, 179.67795088, 147.57375544, 227.0030676 ,
       161.68327739, 289.49584163, 289.22051747, 288.9970385 ,
        67.26417341,  91.4747527 , 131.09747592, 168.08070242,
       256.57964806,  74.30573446,  74.30573446, 206.38728477])


# Create a simple time-longitude plot covering the Indo-Pacific region
plt.figure()

# Markers with size according to by area
plt.scatter(lpt_systems['centroid_lon_stitched'],
         lpt_systems['timestamp_stitched'],
         s = 1e-5*lpt_systems['area_stitched'],
         c = 'none', edgecolors='k', zorder=200,
         linewidths=0.3)

# Draw a line for the centroid track
plt.plot(lpt_systems['centroid_lon_stitched'],
         lpt_systems['timestamp_stitched'], linewidth=1.3,
         color = 'darkorange', zorder=300)

# Set some plotting range limits
plt.gca().set_xlim([60, 240])
plt.gca().set_ylim([dt.datetime(2022,11,1,0,0,0),
                    dt.datetime(2023,5,1,0,0,0)])

# Some labeling
plt.gca().set_xlabel('Longitude')

# Add grid lines
plt.gca().grid(True, axis='both', linewidth=0.5,
               color='darkgrey', linestyle='--')
plt.title('Time-Longitude of LPT System Centroids')

# Save figure
fn_out = 'demo.time_lon_centroids.png'
print(f'--> {fn_out}')
plt.savefig(fn_out, dpi=100, bbox_inches='tight')

--> demo.time_lon_centroids.png


# Plot a map of LPT system tracks, colored by the LPT system ID.
# This can also be done using "scatter" and setting colors to the lptid_stitched
# But here I want to illustrate extracting
# the coordinates for individual LPT systems.
fig = plt.figure()
CMAPNAME='tab10'
cmap = mpl.colormaps[CMAPNAME]


# Set up map background
ax = fig.add_subplot(1,1,1,projection=ccrs.Mercator(central_longitude=180))
ax.coastlines()
ax.set_extent([60, 240, -40, 40], crs=ccrs.PlateCarree())

# Loop over the list of LPT IDs
count = 0
for n, this_lptid in enumerate(lpt_systems['lptid'].data):
    # Choose a color for this LPT system.
    # The colormap 'tab10' has ten colors
    # but we need to feed it a value from 0.0 - 1.0
    # Therefore, get the ramainder on integer division by 10
    # and divide that by 10.0.
    this_color = cmap((n % 10)/10.0)

    # Extract the longitude and latitude for this LPT system
    # Need to be careful since these IDs are decimal values!!!
    #
    # For this case, using the == method seems to work
    # because lptid and lptid_stitched were saved the same
    # way to the NetCDF file.
    #
    # However, this will NOT work when matching these MJO IDs
    # with the values in the MJO list files, because those are
    # saved as decimals in text files and read in using Pandas.
    #
    # USING THE == METHOD WITH LPTID FROM THE MJO LIST FILES 
    # WILL RESULT IN MISSED MJO CASES!!!
    #
    # Recommendation: Always use the <0.00001 method for generality.
    #
    # The == method (not recommended)
    # X = lpt_systems['centroid_lon_stitched'].data[
    #         lpt_systems['lptid_stitched'] == this_lptid]
    # Y = lpt_systems['centroid_lat_stitched'].data[
    #         lpt_systems['lptid_stitched'] == this_lptid]
    #
    # The <0.00001 method (recommended)
    X = lpt_systems['centroid_lon_stitched'].data[
            np.abs(lpt_systems['lptid_stitched'] - this_lptid) < 0.00001]
    Y = lpt_systems['centroid_lat_stitched'].data[
            np.abs(lpt_systems['lptid_stitched'] - this_lptid) < 0.00001]

    if len(X) > 0:
        count += 1

    # Plot the track
    ax.plot(X, Y, color=this_color, transform=ccrs.PlateCarree())

print(count, ' systems.')

# Save figure
fn_out = 'demo.map_centroids.png'
print(f'--> {fn_out}')
plt.savefig(fn_out, dpi=100, bbox_inches='tight')

108  systems.
--> demo.map_centroids.png


fn_mjo_list = ('data/imerg/g50_72h/thresh12/systems/'
               + 'mjo_lpt_list_imerg_2022060100_2023063023.txt')

mjo_list = pd.read_fwf(fn_mjo_list)
mjo_list


# Create a simple time-longitude plot covering the Indo-Pacific region
plt.figure()

# Loop over the list of LPT IDs
count = 0
for n, this_lptid in enumerate(lpt_systems['lptid'].data):

    # Extract the longitude and latitude for this LPT system
    # Need to be careful since these IDs are decimal values!!!
    #
    # For this case, using the == method seems to work
    # because lptid and lptid_stitched were saved the same
    # way to the NetCDF file.
    #
    # However, this will NOT work when matching these MJO IDs
    # with the values in the MJO list files, because those are
    # saved as decimals in text files and read in using Pandas.
    #
    # USING THE == METHOD WITH LPTID FROM THE MJO LIST FILES 
    # WILL RESULT IN MISSED MJO CASES!!!
    #
    # Recommendation: Always use the <0.00001 method for generality.
    #
    # The == method (not recommended)
    # X = lpt_systems['centroid_lon_stitched'].data[
    #         lpt_systems['lptid_stitched'] == this_lptid]
    # Y = lpt_systems['timestamp_stitched'].data[
    #         lpt_systems['lptid_stitched'] == this_lptid]
    #
    # The <0.00001 method (recommended)
    X = lpt_systems['centroid_lon_stitched'].data[
            np.abs(lpt_systems['lptid_stitched'] - this_lptid) < 0.00001]
    Y = lpt_systems['timestamp_stitched'].data[
            np.abs(lpt_systems['lptid_stitched'] - this_lptid) < 0.00001]

    # Check whether it is an MJO LPT system.
    is_mjo = False
    if np.nanmin(np.abs(np.array(mjo_list['lptid'].values) - this_lptid)) < 0.00001:
        is_mjo = True
        count += 1
        print(f'MJO LPT System: {this_lptid}.')
        # Choose a color for this LPT system.
        # The colormap 'tab10' has ten colors
        # but we need to feed it a value from 0.0 - 1.0
        # Therefore, get the ramainder on integer division by 10
        # and divide that by 10.0.
        this_color = cmap((count % 10)/10.0)

    # Plot the track
    if is_mjo:
        plt.gca().plot(X, Y, color=this_color, linewidth=2, zorder=200)
    else:
        plt.gca().plot(X, Y, color='darkgrey', linewidth=1, zorder=100)


# Set some plotting range limits
plt.gca().set_xlim([60, 240])
plt.gca().set_ylim([dt.datetime(2022,11,1,0,0,0),
                    dt.datetime(2023,5,1,0,0,0)])

# Some labeling
plt.gca().set_xlabel('Longitude')

# Add grid lines
plt.gca().grid(True, axis='both', linewidth=0.5,
               color='darkgrey', linestyle='--')
plt.title('Time-Longitude of LPT System Centroids With MJO Highlighted')

# Save figure
fn_out = 'demo.time_lon_centroids.mjo_lpt.png'
print(f'--> {fn_out}')
plt.savefig(fn_out, dpi=100, bbox_inches='tight')

MJO LPT System: 2.7.
MJO LPT System: 4.1.
MJO LPT System: 6.1.
MJO LPT System: 21.4.
MJO LPT System: 21.8.
MJO LPT System: 29.1.
MJO LPT System: 32.1.
MJO LPT System: 37.08.
MJO LPT System: 40.1.
MJO LPT System: 42.1.
MJO LPT System: 44.1.
MJO LPT System: 46.1.
MJO LPT System: 47.1.
MJO LPT System: 48.1.
MJO LPT System: 49.1.
MJO LPT System: 58.1.
--> demo.time_lon_centroids.mjo_lpt.png

	begin_tracking	end_tracking	lptid	lptgroup	duration	mean_zonal_spd	lpt_begin	lpt_end	eprop_begin_idx	eprop_end_idx	eprop_spd	eprop_dur	eprop_begin	eprop_end	eprop_lon_begin	eprop_lon_end
0	2022060100	2023063023	2.70	2	425.0	7.71	2022063017	2022071810	0	415	7.66	415.0	2022063017	2022071800	63.4	169.9
1	2022060100	2023063023	4.10	4	777.0	0.34	2022071714	2022081823	0	520	1.35	520.0	2022071714	2022080806	66.9	103.0
2	2022060100	2023063023	6.10	6	202.0	3.84	2022072110	2022072920	14	200	3.98	186.0	2022072200	2022072918	270.1	291.2
3	2022060100	2023063023	21.80	21	1567.0	-0.45	2022092022	2022112505	972	1495	1.33	523.0	2022103110	2022112205	88.2	127.3
4	2022060100	2023063023	21.40	21	947.0	0.52	2022101618	2022112505	7	308	3.17	301.0	2022101701	2022102914	93.2	124.7
5	2022060100	2023063023	21.80	21	1567.0	-0.45	2022092022	2022112505	358	614	6.58	256.0	2022100520	2022101612	95.5	139.6
6	2022060100	2023063023	29.10	29	955.0	2.54	2022120201	2023011020	140	955	2.93	815.0	2022120721	2023011020	89.2	168.3
7	2022060100	2023063023	32.10	32	884.0	-0.35	2022121616	2023012212	349	758	-0.07	409.0	2023010203	2023011904	290.8	301.3
8	2022060100	2023063023	37.08	37	720.0	4.22	2023012019	2023021919	0	610	4.25	610.0	2023012019	2023021505	91.8	186.0
9	2022060100	2023063023	40.10	40	416.0	0.36	2023021322	2023030306	14	243	2.36	229.0	2023021502	2023022415	292.1	314.1
10	2022060100	2023063023	42.10	42	725.0	0.43	2023030502	2023040407	72	290	3.67	218.0	2023030802	2023031704	285.1	311.3
11	2022060100	2023063023	44.10	44	841.0	5.03	2023032106	2023042507	31	841	5.06	810.0	2023032213	2023042507	77.9	229.6
12	2022060100	2023063023	46.10	46	341.0	0.27	2023041500	2023042905	0	341	0.27	341.0	2023041500	2023042905	294.2	306.9
13	2022060100	2023063023	47.10	47	246.0	1.00	2023042009	2023043015	0	170	6.24	170.0	2023042009	2023042711	3.7	34.4
14	2022060100	2023063023	48.10	48	588.0	0.15	2023042510	2023051922	174	454	2.11	280.0	2023050216	2023051408	71.2	90.2
15	2022060100	2023063023	49.10	49	465.0	0.97	2023042812	2023051721	0	443	0.97	465.0	2023042812	2023051721	154.4	179.7
16	2022060100	2023063023	58.10	58	418.0	2.37	2023061313	2023063023	186	374	8.77	188.0	2023062107	2023062903	127.1	172.1

Demo Python code for Large Scale Precipitation Tracking (LPT)¶

1. Reading in and plotting the centroid tracks¶

1.1. Time-longitude plot of LPT System centroids¶

1.2. Plot a map with individual LPTs highlighted in different colors¶

2. Extract MJO LPT systems¶

2.1. Time-Longitude plot highlighting the MJO LPTs¶