Merge branch 'master' of https://github.com/deeplycloudy/xlma-python into io-improv

Author: wx4stg

.gitignore

@@ -50,6 +50,7 @@ coverage.xml
 *.cover
 .hypothesis/
 .pytest_cache/
+tests/mpl-results/
 
 # Translations
 *.mo

examples/glueviz/config.py

@@ -21,7 +21,7 @@
 #             return line_no
 #
 #     def readfile(self):
-#         lmad = pd.read_csv(self.file,delim_whitespace=True,header=None,
+#         lmad = pd.read_csv(self.file,sep='\\s+',header=None,
 #                                           compression='gzip',skiprows=self.datastart()+1)
 #         columns = ['time','lat','lon','alt','chi','p','mask']
 #         lmad.columns = columns

pyxlma/coords.py

@@ -27,7 +27,7 @@ class CoordinateSystem(object):
         transformations to/from an ECEF cartesian system, e.g.
         >>> self.ERSxyz = proj4.Proj(proj='geocent', ellps='WGS84', datum='WGS84')
         >>> self.ERSlla = proj4.Proj(proj='latlong', ellps='WGS84', datum='WGS84')
-        >>> projectedData = proj4.transform(self.ERSlla, self.ERSxyz, lat, lon, alt )
+        >>> projectedData = proj4.Transformer.from_crs(self.ERSlla.crs, self.ERSxyz.crs).transform(lon, lat, alt)
     The ECEF system has its origin at the center of the earth, with the +Z toward the north pole,
         +X toward (lat=0, lon=0), and +Y right-handed orthogonal to +X, +Z
 
@@ -80,7 +80,7 @@ def toECEF(self, lon, lat, alt):
         lon = atleast_1d(lon)
         alt = atleast_1d(alt)
         if (lat.shape[0] == 0): return lon, lat, alt # proj doesn't like empties
-        projectedData = array(proj4.transform(self.ERSlla, self.ERSxyz, lon, lat, alt ))
+        projectedData = array(proj4.Transformer.from_crs(self.ERSlla.crs, self.ERSxyz.crs).transform(lon, lat, alt))
         if len(projectedData.shape) == 1:
             return projectedData[0], projectedData[1], projectedData[2]
         else:
@@ -91,7 +91,7 @@ def fromECEF(self, x, y, z):
         y = atleast_1d(y)
         z = atleast_1d(z)
         if (x.shape[0] == 0): return x, y, z # proj doesn't like empties
-        projectedData = array(proj4.transform(self.ERSxyz, self.ERSlla, x, y, z ))
+        projectedData = array(proj4.Transformer.from_crs(self.ERSxyz.crs, self.ERSlla.crs).transform(x, y, z))
         if len(projectedData.shape) == 1:
             return projectedData[0], projectedData[1], projectedData[2]
         else:
@@ -125,15 +125,15 @@ def toECEF(self, x, y, z):
         x += self.cx
         y += self.cy
         z += self.cz
-        projectedData = array(proj4.transform(self.projection, self.ERSxyz, x, y, z ))
+        projectedData = array(proj4.Transformer.from_crs(self.projection.crs, self.ERSxyz.crs).transform(x, y, z))
         if len(projectedData.shape) == 1:
             px, py, pz = projectedData[0], projectedData[1], projectedData[2]
         else:
             px, py, pz = projectedData[0,:], projectedData[1,:], projectedData[2,:]
         return px, py, pz
 
     def fromECEF(self, x, y, z):
-        projectedData = array(proj4.transform(self.ERSxyz, self.projection, x, y, z ))
+        projectedData = array(proj4.Transformer.from_crs(self.ERSxyz.crs, self.projection.crs).transform(x, y, z))
         if len(projectedData.shape) == 1:
             px, py, pz = projectedData[0], projectedData[1], projectedData[2]
         else:
@@ -221,10 +221,10 @@ def __init__(self, subsat_lon=0.0, subsat_lat=0.0, sweep_axis='y',
 
     def toECEF(self, x, y, z):
         X, Y, Z = x*self.h, y*self.h, z*self.h
-        return proj4.transform(self.fixedgrid, self.ECEFxyz, X, Y, Z)
+        return proj4.Transformer.from_crs(self.fixedgrid.crs, self.ECEFxyz.crs).transform(X, Y, Z)
 
     def fromECEF(self, x, y, z):
-        X, Y, Z = proj4.transform(self.ECEFxyz, self.fixedgrid, x, y, z)
+        X, Y, Z = proj4.Transformer.from_crs(self.ECEFxyz.crs, self.fixedgrid.crs).transform(x, y, z)
         return X/self.h, Y/self.h, Z/self.h
 
 # class AltitudePreservingMapProjection(MapProjection):
@@ -253,8 +253,8 @@ def __init__(self, ctrLat, ctrLon, ctrAlt, datum='WGS84', ellps='WGS84', effecti
         self.lla = proj4.Proj(proj='latlong', ellps=self.ellps, datum=self.datum)
         self.xyz = proj4.Proj(proj='geocent', ellps=self.ellps, datum=self.datum)
 
-        self.Requator, foo1, foo2 = proj4.transform(self.lla,self.xyz,0,0,0) # Equatorial radius  - WGS-84 value = 6378137.0
-        foo1, foo2, self.Rpolar = proj4.transform(self.lla,self.xyz,0,90,0) # Polar radius  - WGS-84 value = 6356752.314
+        self.Requator, _, _ = proj4.Transformer.from_crs(self.lla.crs, self.xyz.crs).transform(0,0,0) # Equatorial radius  - WGS-84 value = 6378137.0
+        _, _, self.Rpolar = proj4.Transformer.from_crs(self.lla.crs, self.xyz.crs).transform(0,90,0) # Polar radius  - WGS-84 value = 6356752.314
         self.flattening = (self.Requator-self.Rpolar)/self.Requator
 
         self.eccen = (2.0-self.flattening)*self.flattening   # First eccentricity squared - WGS-84 value = 0.00669437999013
@@ -374,10 +374,10 @@ def __init__(self, ctrLat=0.0, ctrLon=0.0, ctrAlt=0.0):
 
         ERSlla = proj4.Proj(proj='latlong', ellps='WGS84', datum='WGS84')
         ERSxyz = proj4.Proj(proj='geocent',  ellps='WGS84', datum='WGS84')
-        self.centerECEF = array(proj4.transform(ERSlla, ERSxyz, ctrLon, ctrLat, ctrAlt))
+        self.centerECEF = array(proj4.Transformer.from_crs(ERSlla.crs, ERSxyz.crs).transform(ctrLon, ctrLat, ctrAlt))
 
         #location of point directly above local center
-        aboveCenterECEF = array(proj4.transform(ERSlla, ERSxyz, ctrLon, ctrLat, self.ctrAlt+1e3))
+        aboveCenterECEF = array(proj4.Transformer.from_crs(ERSlla.crs, ERSxyz.crs).transform(ctrLon, ctrLat, self.ctrAlt+1e3))
 
         #normal vector to earth's surface at the center is the local z direction
         n = aboveCenterECEF - self.centerECEF
@@ -396,7 +396,7 @@ def __init__(self, ctrLat=0.0, ctrLon=0.0, ctrAlt=0.0):
         # Point just to the north of the center on earth's surface, projected onto the tangent plane
         # This calculation seems like it should only be done with latitude/north since the local x
         #   direction curves away along a non-straight line when projected onto the plane
-        northCenterECEF = array(proj4.transform(ERSlla, ERSxyz, self.ctrLon, self.ctrLat+1.01, self.ctrAlt))
+        northCenterECEF = array(proj4.Transformer.from_crs(ERSlla.crs, ERSxyz.crs).transform(self.ctrLon, self.ctrLat+1.01, self.ctrAlt))
         localy = dot(P, northCenterECEF[:,None] )
         localy = localy / norm(localy)
 

pyxlma/lmalib/flash/properties.py

@@ -1,8 +1,7 @@
 import numpy as np
 import xarray as xr
-from scipy.spatial import Delaunay, ConvexHull
+from scipy.spatial import Delaunay, ConvexHull, QhullError
 from scipy.special import factorial
-from scipy.spatial.qhull import QhullError
 from pyxlma.lmalib.traversal import OneToManyTraversal
 
 def local_cartesian(lon, lat, alt, lonctr, latctr, altctr):
@@ -135,7 +134,7 @@ def event_discharge_energy(z,area):
     eta_c      = 0.004 #Scale the energy to depict the fraction of energy neutralized by
                        #each flash in the capacitor model
     #Capacitor model:
-    w = 4 * ((sigma_crit**2. * d * area.iloc[0])/(2* e)) #The quantity for appears when considering image charges (2*sigma)^2=4sigma^2
+    w = 4 * ((sigma_crit**2. * d * area)/(2* e)) #The quantity for appears when considering image charges (2*sigma)^2=4sigma^2
     return(w*eta_c)
 
 
@@ -240,12 +239,12 @@ def flash_stats(ds, area_func=None, volume_func=None):
     # could only run on those where point counts meet threshold, instead
     # testing inside the function
     # Index of event_area and event_volume are event_parent_flash_id
-    event_area = fl_gb.apply(area_func)
-    event_volume = fl_gb.apply(volume_func)
+    event_area = fl_gb.apply(area_func, include_groups=False)
+    event_volume = fl_gb.apply(volume_func, include_groups=False)
 
     #Compute flash discharge energy using parallel plate capacitor
-    event_energy = fl_gb.apply(lambda df: event_discharge_energy(df['event_z'],
-                                                                 event_area[df['event_parent_flash_id']]))
+    event_energy = fl_gb.apply(lambda df1: event_discharge_energy(df1['event_z'],
+                                                                 event_area[df1.name]), include_groups=False)
 
 
     # set the index for the original dataset's flash dimension to the flash_id

pyxlma/lmalib/grid.py

@@ -132,7 +132,7 @@ def assign_regular_bins(dsg, ds, var_to_grid_map, pixel_id_var='pixel_id',
             var_dim = ds[var_name].dims[0]
         else:
             assert var_dim == ds[var_name].dims[0]
-        if ds.dims[var_dim] < 1:
+        if ds.sizes[var_dim] < 1:
             # No data on the dim, set in_range to nothing and stop checking other vars
             in_range = []
             have_data = False

pyxlma/lmalib/io/read.py

@@ -413,7 +413,7 @@ def nldn(filenames):
         filenames = [filenames]
     full_df = pd.DataFrame({})
     for filename in filenames:
-        this_file = pd.read_csv(filename, delim_whitespace=True, header=None, 
+        this_file = pd.read_csv(filename, sep='\\s+', header=None, 
                     names=[
                         'date', 'time', 'latitude', 'longitude', 'peak_current_kA', 'curr_unit', 'multiplicity', 'semimajor',
                         'semiminor', 'majorminorratio', 'ellipseangle', 'chi2', 'num_stations', 'type'
@@ -629,7 +629,7 @@ def readfile(self):
                 comp = 'gzip'
             else:
                 comp = None
-            lmad = pd.read_csv(self.file,compression=comp,delim_whitespace=True,
+            lmad = pd.read_csv(self.file,compression=comp,sep='\\s+',
                                header=None,skiprows=self.data_starts+1,on_bad_lines='skip')
             lmad.columns = self.names
         except pd.errors.EmptyDataError:

pyxlma/lmalib/traversal.py

@@ -77,7 +77,7 @@ def __init__(self, dataset, entity_id_vars, parent_id_vars):
         self.child_to_parent = collections.OrderedDict()
         self.parent_to_child = collections.OrderedDict()
         for (entity_var, parent_var) in self._descend():
-            if dataset.dims[dataset[entity_var].dims[0]] == 0:
+            if dataset.sizes[dataset[entity_var].dims[0]] == 0:
                 # No data, so groupby will fail in xarray > 0.13
                 entity_grouper = None
             else:
@@ -86,7 +86,7 @@ def __init__(self, dataset, entity_id_vars, parent_id_vars):
             if parent_var is None:
                 parent_grouper = None
             else:
-                if dataset.dims[dataset[parent_var].dims[0]] == 0:
+                if dataset.sizes[dataset[parent_var].dims[0]] == 0:
                     # No data, so groupby will fail in xarray > 0.13
                     parent_grouper = None
                 else:

pyxlma/plot/xlma_plot_feature.py

@@ -161,8 +161,8 @@ def plot_2d_network_points(bk_plot, netw_data, actual_height=None, fake_ic_heigh
     else:
         raise ValueError("color_by must be 'time' or 'polarity'")
 
-    cgs = netw_data[netw_data['type']=='CG']
-    ics = netw_data[netw_data['type']=='IC']
+    cgs = netw_data[netw_data['type']=='CG'].copy()
+    ics = netw_data[netw_data['type']=='IC'].copy()
 
     if actual_height is None:
         cgs['height'] = np.full_like(cgs.longitude, fake_cg_height)

tests/test_coords.py

@@ -20,7 +20,6 @@ def test_geographic():
     geosys = GeographicSystem()
     ecef_coords = geosys.toECEF(test_lons, test_lats, test_alts)
     lons, lats, alts = geosys.fromECEF(*ecef_coords)
-
     assert np.allclose(ecef_coords[0], test_ecef_X)
     assert np.allclose(ecef_coords[1], test_ecef_Y)
     assert np.allclose(ecef_coords[2], test_ecef_Z)
@@ -31,9 +30,7 @@ def test_geographic():
 def test_geographic_one_point():
     geosys = GeographicSystem()
     ecef_coords = geosys.toECEF(np.atleast_1d(test_lons[-1]), np.atleast_1d(test_lats[-1]), np.atleast_1d(test_alts[-1]))
-    print(len(np.atleast_1d(test_lons[-1]).shape))
     lons, lats, alts = geosys.fromECEF(*ecef_coords)
-
     assert np.allclose(ecef_coords[0], test_ecef_X[-1])
     assert np.allclose(ecef_coords[1], test_ecef_Y[-1])
     assert np.allclose(ecef_coords[2], test_ecef_Z[-1])
@@ -119,7 +116,6 @@ def test_radar_system_height():
 def test_radar_system_elevation():
     ADRAD_rcs = RadarCoordinateSystem(30.6177, -96.3365, 114)
     tornado_slant_range, radar_elevation = ADRAD_rcs.getSlantRangeElevation(17144.013390611748, 550.2784673999995)
-    print(tornado_slant_range, radar_elevation)
     assert np.allclose(tornado_slant_range, 17150)
     assert np.allclose(radar_elevation, 1.4)
 

tests/test_grid.py

@@ -10,8 +10,8 @@ def test_create_regular_grid():
     grid_h_res = 0.1
     grid_height = 20
     grid_v_res = 1
-    lon_range = (dataset.network_center_longitude - grid_range, dataset.network_center_longitude + grid_range, grid_h_res)
-    lat_range = (dataset.network_center_latitude - grid_range, dataset.network_center_latitude + grid_range, grid_h_res)
+    lon_range = (dataset.network_center_longitude.data.item() - grid_range, dataset.network_center_longitude.data.item() + grid_range, grid_h_res)
+    lat_range = (dataset.network_center_latitude.data.item() - grid_range, dataset.network_center_latitude.data.item() + grid_range, grid_h_res)
     alt_range = (0, grid_height, grid_v_res)
     time_range = (dataset.event_time.data.min(), dataset.event_time.data.max(), np.timedelta64(1, 'm'))
     grid_edge_ranges ={
@@ -66,4 +66,4 @@ def test_events_to_grid():
     truth = xr.open_dataset('tests/truth/lma_netcdf/gridded_lma.nc')
 
     for var in truth.data_vars:
-        compare_dataarrays(gridded_lma, truth, var)
+        compare_dataarrays(gridded_lma, truth, var)