Merge pull request #34 from wx4stg/nldn

Add NLDN/ENTLN plotting

Author: deeplycloudy

examples/network_samples/gld360enldnns_20231224_daily_v1_lit.raw

@@ -0,0 +1,3 @@
+12/24/23 00:57:01.123456789  33.582 -101.881  +12.0 kA  0 0.4  0.2  1.8   75  0.8  4 C
+12/24/23 00:57:31.987654321  33.590 -102.032  -9.0 kA  0 0.2  0.2  1.0    3  1.0  5 G
+12/24/23 00:57:58.135792468  33.585 -101.872  -15.0 kA  0 2.2  1.5  1.5   26  2.1  6 C

examples/network_samples/lxarchive_pulse20231224.csv

@@ -0,0 +1,4 @@
+type,timestamp,latitude,longitude,peakcurrent,icheight,numbersensors,majoraxis,minoraxis,bearing
+ 1, 2023-12-24T00:57:04.123456789,  33.581914, -101.880986,   12345, 3014,   8, 371.00, 58.00, 104.00  
+ 0, 2023-12-24T00:57:26.987654321,  33.590077, -102.032033,   9876, 6028,   7, 693.00, 108.00, 99.00  
+ 1, 2023-12-24T00:57:47.246813579,  33.584480, -101.871498,   15790, 13591,  11, 178.00, 62.00, 102.00  

pyxlma/lmalib/io/read.py

@@ -220,6 +220,116 @@ def to_dataset(lma_file, event_id_start=0):
     return ds
 
 
+def nldn(filenames):
+    """
+    Read Viasala NLDN data
+     
+    Reads in one or multiple NLDN files and and returns a pandas dataframe with appropriate column names
+
+
+    Parameters
+    ----------
+    filenames : str or list of str
+        The file or files to read in
+    
+    
+    Returns
+    -------
+    full_df : `pandas.DataFrame`
+        A pandas dataframe of entln data, with columns:
+        'latitude' - the latitude of the event
+        'longitude' - the longitude of the event
+        'peak_current_kA' - the peak current in kA
+        'multiplicity' - the number of strokes in the event
+        'semimajor' - the semimajor axis length in km of the 50% confidence ellipse
+        'semiminor' - the semiminor axis length in km of the 50% confidence ellipse
+        'ellipseangle' - the angle of the 50% confidence ellipse
+        'chi2' - the reduced chi-squared value of the event
+        'num_stations' - the number of stations contributing to the event
+        'type' - 'IC' or 'CG' for intracloud or cloud-to-ground
+        'datetime' - the time of the event
+
+        
+    Notes
+    -----
+    This function is designed to read NLDN files in the format that matches the format of the file located
+    in examples/network_samples/gld360enldnns_20231224_daily_v1_lit.raw (This is not real NLDN data, but provides a correct sample of the format)
+
+    Other file formats may exist and may not be read in correctly. If you have a file that is not read in correctly,
+    please open an issue on the pyxlma GitHub page.
+    """
+    if type(filenames) is str:
+        filenames = [filenames]
+    full_df = pd.DataFrame({})
+    for filename in filenames:
+        this_file = pd.read_csv(filename, delim_whitespace=True, header=None, 
+                    names=[
+                        'date', 'time', 'latitude', 'longitude', 'peak_current_kA', 'curr_unit', 'multiplicity', 'semimajor',
+                        'semiminor', 'majorminorratio', 'ellipseangle', 'chi2', 'num_stations', 'type'
+                        ])
+        if len(this_file['curr_unit'].drop_duplicates()) == 1:
+            this_file.drop(columns=['curr_unit'], inplace=True)
+        else:
+            raise ValueError('Multiple current units in file')
+        this_file['datetime'] = pd.to_datetime(this_file['date']+' '+this_file['time'], format='%m/%d/%y %H:%M:%S.%f')
+        this_file.drop(columns=['date','time'], inplace=True)
+        this_file['type'] = this_file['type'].map({'G':'CG','C':'IC'})
+        full_df = pd.concat([full_df, this_file])
+    return full_df
+
+
+def entln(filenames):
+    """
+    Read Earth Networks Total Lightning Network data
+     
+    Reads in one or multiple ENTLN files and and returns a pandas dataframe with appropriate column names
+
+
+    Parameters
+    ----------
+    filenames : str or list of str
+        The file or files to read in
+    
+    
+    Returns
+    -------
+    full_df : `pandas.DataFrame`
+        A pandas dataframe of entln data, with columns:
+        'type' - 'IC' or 'CG' for intracloud or cloud-to-ground
+        'datetime' - the time of the event
+        'latitude' - the latitude of the event
+        'longitude' - the longitude of the event
+        'peak_current_kA' - the peak current in kA
+        'icheight' - the height of the IC event in meters
+        'num_stations' - the number of stations contributing to the event
+        'ellipseangle' - the angle of the 50% confidence ellipse
+        'semimajor' - the semimajor axis length in km of the 50% confidence ellipse
+        'semiminor' - the semiminor axis length in km of the 50% confidence ellipse
+
+    Notes
+    -----
+    This function is designed to read ENTLN files in the format that matches the format of the file located
+    in examples/network_samples/lxarchive_pulse20231224.csv (This is not real ENTLN data, but provides a correct sample of the format)
+
+    Other file formats may exist and may not be read in correctly. If you have a file that is not read in correctly,
+    please open an issue on the pyxlma GitHub page.
+    """
+    if type(filenames) is str:
+        filenames = [filenames]
+    full_df = pd.DataFrame({})
+    for filename in filenames:
+        this_file = pd.read_csv(filename, parse_dates=['timestamp'])
+        this_file['peakcurrent'] = this_file['peakcurrent']/1000
+        this_file['type'] = this_file['type'].map({0:'CG',1:'IC'})
+        this_file['semimajor'] = this_file['majoraxis']/2000
+        this_file['semiminor'] = this_file['minoraxis']/2000
+        this_file.drop(columns=['majoraxis','minoraxis'], inplace=True)
+        rename = {'timestamp' : 'datetime', 'peakcurrent' : 'peak_current_kA', 'numbersensors' : 'num_stations', 'bearing' : 'ellipseangle'}
+        this_file.rename(columns=rename, inplace=True)
+        full_df = pd.concat([full_df, this_file])
+    return full_df
+
+    
 class lmafile(object):
     def __init__(self,filename):
         """

pyxlma/plot/xlma_plot_feature.py

@@ -63,7 +63,7 @@ def setup_hist(lon_data, lat_data, alt_data, time_data,
 
 def plot_points(bk_plot, lon_data, lat_data, alt_data, time_data,
                   plot_cmap=None, plot_s=None, plot_vmin=None, plot_vmax=None, plot_c=None, edge_color='face',
-                  edge_width=0, add_to_histogram=True, **kwargs):
+                  edge_width=0, add_to_histogram=True, marker='o', **kwargs):
     """
     Plot scatter points on an existing bk_plot object given x,y,z,t for each
     and defined plotting colormaps and ranges
@@ -88,16 +88,16 @@ def plot_points(bk_plot, lon_data, lat_data, alt_data, time_data,
 
     art_plan = bk_plot.ax_plan.scatter(lon_data, lat_data,
                             c=plot_c,vmin=plot_vmin, vmax=plot_vmax, cmap=plot_cmap,
-                            s=plot_s,marker='o', linewidths=edge_width, edgecolors=edge_color, **kwargs)
+                            s=plot_s, marker=marker, linewidths=edge_width, edgecolors=edge_color, **kwargs)
     art_th = bk_plot.ax_th.scatter(time_data, alt_data,
                           c=plot_c,vmin=plot_vmin, vmax=plot_vmax, cmap=plot_cmap,
-                          s=plot_s,marker='o', linewidths=edge_width, edgecolors=edge_color, **kwargs)
+                          s=plot_s, marker=marker, linewidths=edge_width, edgecolors=edge_color, **kwargs)
     art_lon = bk_plot.ax_lon.scatter(lon_data, alt_data,
                           c=plot_c,vmin=plot_vmin, vmax=plot_vmax, cmap=plot_cmap,
-                          s=plot_s,marker='o',  linewidths=edge_width, edgecolors=edge_color, **kwargs)
+                          s=plot_s, marker=marker, linewidths=edge_width, edgecolors=edge_color, **kwargs)
     art_lat = bk_plot.ax_lat.scatter(alt_data, lat_data,
                           c=plot_c,vmin=plot_vmin, vmax=plot_vmax, cmap=plot_cmap,
-                          s=plot_s,marker='o', linewidths=edge_width, edgecolors=edge_color, **kwargs)
+                          s=plot_s, marker=marker, linewidths=edge_width, edgecolors=edge_color, **kwargs)
     art_out = [art_plan, art_th, art_lon, art_lat]
 
     if add_to_histogram:
@@ -111,6 +111,84 @@ def plot_points(bk_plot, lon_data, lat_data, alt_data, time_data,
         art_out.append(art_hist)
     return art_out
 
+
+def plot_2d_network_points(bk_plot, netw_data, actual_height=None, fake_ic_height=18, fake_cg_height=1,
+                        color_by='time', pos_color='red', neg_color='blue', **kwargs):
+    """
+    Plot points from a 2D lightning mapping neworks (ie, NLDN, ENTLN, etc)
+
+    Parameters
+    ----------
+    bk_plot : `pyxlma.plot.xlma_base_plot.BlankPlot`
+        A BlankPlot object to plot the data on
+    netw_data : `pandas.DataFrame` or `xarray.Dataset`
+        data object with columns/variables 'longitude', 'latitude', 'type' (CG/IC), and 'datetime'
+    actual_height : `numpy.ndarray` or `pandas.Series` or `xarray.DataArray`
+        the hieghts of the events to be plotted (default None, fake_ic_height and fake_cg_height used)
+    fake_ic_height : float
+        the altitude to plot IC points (default 18 km)
+    fake_cg_height : float
+        the altitude to plot CG points (default 1 km)
+    color_by : ['time', 'polarity']
+        Whether to color the points by time or polarity. Default 'time'. Ignored if **kwargs contains 'c'.
+    pos_color : str
+        color for positive points (default 'red') if color_by='polarity'
+    neg_color : str
+        color for negative points (default 'blue') if color_by='polarity'
+    **kwargs
+        additional keyword arguments to pass to plt.scatter
+
+    Returns
+    -------
+    art_out : list
+        nested lists of artists created by plot_points (first list CG, second list IC)
+
+    """
+
+    plot_c = kwargs.pop('c', None)
+    vmin = kwargs.pop('vmin', None)
+    vmax = kwargs.pop('vmax', None)
+    marker = kwargs.pop('marker', '^')
+    if actual_height is not None:
+        netw_data['height'] = actual_height
+
+    if plot_c is not None:
+        netw_data['plot_c'] = plot_c
+    elif color_by == 'time':
+        netw_data['plot_c'] = color_by_time(netw_data.datetime, bk_plot.tlim)[2]
+    elif color_by == 'polarity':
+        pass
+    else:
+        raise ValueError("color_by must be 'time' or 'polarity'")
+    
+    cgs = netw_data[netw_data['type']=='CG']
+    ics = netw_data[netw_data['type']=='IC']
+
+    if actual_height is None:
+        cgs['height'] = np.full_like(cgs.longitude, fake_cg_height)
+        ics['height'] = np.full_like(ics.longitude, fake_ic_height)
+    art_out = []
+    if color_by == 'polarity':
+        cgpos = cgs[cgs.peak_current_kA>0]
+        cgneg = cgs[cgs.peak_current_kA<0]
+        icpos = ics[ics.peak_current_kA>0]
+        icneg = ics[ics.peak_current_kA<0]
+        art_out.append(plot_points(bk_plot, cgneg.longitude, cgneg.latitude, cgneg.height,
+                    cgneg.datetime, c=neg_color, marker=marker, add_to_histogram=False, **kwargs))
+        art_out.append(plot_points(bk_plot, cgpos.longitude, cgpos.latitude, cgpos.height,
+                    cgpos.datetime, c=pos_color, marker=marker, add_to_histogram=False, **kwargs))
+        art_out.append(plot_points(bk_plot, icneg.longitude, icneg.latitude, icneg.height,
+                    icneg.datetime, c=neg_color, marker=marker, add_to_histogram=False, **kwargs))
+        art_out.append(plot_points(bk_plot, icpos.longitude, icpos.latitude, icpos.height,
+                    icpos.datetime, c=pos_color, marker=marker, add_to_histogram=False, **kwargs))
+    else:
+        art_out.append(plot_points(bk_plot, cgs.longitude, cgs.latitude, cgs.height,
+                    cgs.datetime, c=cgs.plot_c, vmin=vmin, vmax=vmax, marker=marker, add_to_histogram=False, **kwargs))
+        art_out.append(plot_points(bk_plot, ics.longitude, ics.latitude, ics.height,
+                    ics.datetime, c=ics.plot_c, vmin=vmin, vmax=vmax, marker=marker, add_to_histogram=False, **kwargs))
+    return art_out
+
+
 def plot_3d_grid(bk_plot, xedges, yedges, zedges, tedges,
                 alt_lon, alt_lat, alt_time, lat_lon,
                 alt_data, plot_cmap=None, **kwargs):

tests/test_base_plot.py

@@ -3,6 +3,7 @@
 from pyxlma.plot.xlma_base_plot import *
 from pyxlma.plot.xlma_plot_feature import *
 from pyxlma.lmalib.grid import *
+from pyxlma.lmalib.io import read as lma_read
 import datetime as dt
 import pandas as pd
 import matplotlib.dates as md
@@ -135,4 +136,53 @@ def test_plot_feature_inset_view():
                 c=colors, edgecolors='k', linewidths=0.1, add_to_histogram=True)
     inset_view(bk_plot, dataset.event_longitude.data, dataset.event_latitude.data,
                [-102.75, -102.25], [32, 32.5], .01, .01)
-    return bk_plot.fig
+    return bk_plot.fig
+
+
+@pytest.mark.mpl_image_compare
+def test_plot_feature_ndln_time():
+    start_time = dt.datetime(2023, 12, 24, 0, 57, 0, 0)
+    end_time = start_time + dt.timedelta(seconds=60)
+    bk_plot = BlankPlot(start_time, bkgmap=True, xlim=[-103.5, -99.5], ylim=[31.5, 35.5], zlim=[0, 20], tlim=[start_time, end_time], title='XLMA Test Plot')
+    nldn_data = lma_read.nldn('examples/network_samples/gld360enldnns_20231224_daily_v1_lit.raw')
+    plot_2d_network_points(bk_plot, nldn_data)
+    return bk_plot.fig
+
+
+def test_plot_feature_ndln_bad_arg():
+    start_time = dt.datetime(2023, 12, 24, 0, 57, 0, 0)
+    end_time = start_time + dt.timedelta(seconds=60)
+    bk_plot = BlankPlot(start_time, bkgmap=True, xlim=[-103.5, -99.5], ylim=[31.5, 35.5], zlim=[0, 20], tlim=[start_time, end_time], title='XLMA Test Plot')
+    nldn_data = lma_read.nldn('examples/network_samples/gld360enldnns_20231224_daily_v1_lit.raw')
+    with pytest.raises(ValueError, match="color_by must be 'time' or 'polarity'"):
+        plot_2d_network_points(bk_plot, nldn_data, color_by='bad_arg')
+
+
+@pytest.mark.mpl_image_compare
+def test_plot_feature_ndln_custom_colors():
+    start_time = dt.datetime(2023, 12, 24, 0, 57, 0, 0)
+    end_time = start_time + dt.timedelta(seconds=60)
+    bk_plot = BlankPlot(start_time, bkgmap=True, xlim=[-103.5, -99.5], ylim=[31.5, 35.5], zlim=[0, 20], tlim=[start_time, end_time], title='XLMA Test Plot')
+    nldn_data = lma_read.nldn('examples/network_samples/gld360enldnns_20231224_daily_v1_lit.raw')
+    plot_2d_network_points(bk_plot, nldn_data, c=[0, 10, 5], cmap='plasma', vmin=0, vmax=10)
+    return bk_plot.fig
+
+
+@pytest.mark.mpl_image_compare
+def test_plot_feature_entln_real_height():
+    start_time = dt.datetime(2023, 12, 24, 0, 57, 0, 0)
+    end_time = start_time + dt.timedelta(seconds=60)
+    bk_plot = BlankPlot(start_time, bkgmap=True, xlim=[-103.5, -99.5], ylim=[31.5, 35.5], zlim=[0, 20], tlim=[start_time, end_time], title='XLMA Test Plot')
+    entln_data = lma_read.entln('examples/network_samples/lxarchive_pulse20231224.csv')
+    plot_2d_network_points(bk_plot, entln_data, actual_height=entln_data['icheight'])
+    return bk_plot.fig
+
+
+@pytest.mark.mpl_image_compare
+def test_plot_feature_ndln_polarity():
+    start_time = dt.datetime(2023, 12, 24, 0, 57, 0, 0)
+    end_time = start_time + dt.timedelta(seconds=60)
+    bk_plot = BlankPlot(start_time, bkgmap=True, xlim=[-103.5, -99.5], ylim=[31.5, 35.5], zlim=[0, 20], tlim=[start_time, end_time], title='XLMA Test Plot')
+    nldn_data = lma_read.nldn('examples/network_samples/gld360enldnns_20231224_daily_v1_lit.raw')
+    plot_2d_network_points(bk_plot, nldn_data, color_by='polarity')
+    return bk_plot.fig

tests/test_io.py

@@ -3,6 +3,8 @@
 from os import listdir
 from datetime import datetime as dt
 import xarray as xr
+import pandas as pd
+import numpy as np
 from pyxlma.lmalib.io import read as lma_read
 from tests.test_flash import compare_dataarrays
 
@@ -22,3 +24,39 @@ def test_read_onefile_dataset():
     truth = xr.open_dataset('tests/truth/lma_netcdf/small_lma.nc')
     for var in truth.data_vars:
         compare_dataarrays(dataset, truth, var)
+
+
+def test_read_nldn():
+    dataset = lma_read.nldn('examples/network_samples/gld360enldnns_20231224_daily_v1_lit.raw')
+    truth = pd.DataFrame({
+        'latitude' : [33.582, 33.590, 33.585],
+        'longitude' : [-101.881, -102.032, -101.872],
+        'peak_current_kA' : [12.0, -9.0, -15.0],
+        'multiplicity' : [0, 0, 0],
+        'semimajor': [0.4, 0.2, 2.2],
+        'semiminor': [0.2, 0.2, 1.5],
+        'majorminorratio': [1.8, 1.0, 1.5],
+        'ellipseangle' : [75, 3, 26],
+        'chi2' : [0.8, 1.0, 2.1],
+        'num_stations': [4, 5, 6],
+        'type': ['IC', 'CG', 'IC'],
+        'datetime' : [np.datetime64('2023-12-24T00:57:01.123456789'), np.datetime64('2023-12-24T00:57:31.987654321'), np.datetime64('2023-12-24T00:57:58.135792468')]
+    })
+    assert dataset.equals(truth)
+
+
+def test_read_entln():
+    dataset = lma_read.entln('examples/network_samples/lxarchive_pulse20231224.csv')
+    truth = pd.DataFrame({
+        'type': ['IC', 'CG', 'IC'],
+        'datetime' : [np.datetime64('2023-12-24T00:57:04.123456789'), np.datetime64('2023-12-24T00:57:26.987654321'), np.datetime64('2023-12-24T00:57:47.246813579')],
+        'latitude' : [33.581914, 33.590077, 33.584480],
+        'longitude' : [-101.880986, -102.032033, -101.871498],
+        'peak_current_kA' : [12.345, 9.876, 15.79],
+        'icheight' : [3014, 6028, 13591],
+        'num_stations': [8, 7, 11],
+        'ellipseangle' : [104., 99., 102.],
+        'semimajor': [0.1855, 0.3465, 0.089],
+        'semiminor': [0.029, 0.054, 0.031]
+    })
+    assert dataset.equals(truth)