merged master branch commits to this branch

Merge branch 'master' into mergeday

.github/workflows/test_python.yml

@@ -61,7 +61,7 @@ jobs:
   lint:
 
     runs-on: ubuntu-latest
-    continue-on-error: true
+    #continue-on-error: true
     strategy:
       fail-fast: false
 
@@ -78,6 +78,6 @@ jobs:
     - name: Lint with flake8
       run: |
         # stop the build if there are Python syntax errors or undefined names
-        flake8 . --count --select=E9,F63,F7,F82 --show-source --statistics
+        flake8 . --count --select=E9,F63,F7,F82 --show-source --statistics --extend-exclude=trash_can,pyMayaVi
         # exit-zero treats all errors as warnings. The GitHub editor is 127 chars wide
-        flake8 . --count --exit-zero --max-complexity=10 --max-line-length=127 --statistics
+        flake8 . --count --exit-zero --max-complexity=10 --max-line-length=127 --statistics --extend-exclude=trash_can,pyMayaVi

examples/generate_panel.py

@@ -188,8 +188,8 @@ def cavitoncontours(ax, XmeshXY,YmeshXY, extmaps, requestvariables=False):
 # it has access to the name of the .vlsv file, and then calls plot_vdf to
 # create insets on top of the variable map.
 def insetVDF(ax, XmeshXY,YmeshXY, pass_maps):
-    if requestvariables==True:
-        return []    
+    # if requestvariables==True: # Never triggers, copy-paste leftovers?
+    #     return []    
     # pass_maps is a list of numpy arrays, not used here.
     from mpl_toolkits.axes_grid1.inset_locator import inset_axes
 
@@ -365,7 +365,7 @@ def insetVDF(ax, XmeshXY,YmeshXY, pass_maps):
                           vmin=1e-2, 
                           vmax=1e0, 
                           pass_vars=['E','B','V'], 
-                          expression=expr_Slippage)
+                          expression=pt.plot.plot_helpers.expr_Slippage)
 
 
     # Plot beam number density with inset VDF

examples/multi_panel_plot.py

@@ -7,7 +7,7 @@
 
 
 # Avoids opening a figure window
-if str(matplotlib.get_backend()) is not 'Agg':
+if str(matplotlib.get_backend()) != 'Agg':
     plt.switch_backend('Agg') 
 
 fileLocation="/proj/vlasov/2D/BCH/bulk/"

pyCalculations/cut3d.py

@@ -62,8 +62,8 @@ def cut3d( vlsvReader, xmin, xmax, ymin, ymax, zmin, zmax, variable, operator="p
    # Read the cell lengths:
    ##################################################
    # Get xmax, xmin and xcells_ini
-   mesh_limits = f.get_spatial_mesh_extent()
-   mesh_size = f.get_spatial_mesh_size()
+   mesh_limits = vlsvReader.get_spatial_mesh_extent()
+   mesh_size = vlsvReader.get_spatial_mesh_size()
    xmax = mesh_limits[3]
    xmin = mesh_limits[0]
    xcells = mesh_size[0]

pyCalculations/find_x_and_o.py

@@ -62,7 +62,7 @@ def findIntersection(v1,v2):
 elif len(sys.argv)==3:
    indexes = np.arange(int(sys.argv[1]),int(sys.argv[2]))
 elif len(sys.argv)==4:
-   indexes = np.arange(int(sys.argv[1]),int(sys.argv[2]),,int(sys.argv[3]))
+   indexes = np.arange(int(sys.argv[1]),int(sys.argv[2]),int(sys.argv[3]))
 else:
    sys.stderr.write("Syntax: size.py <index>\n")
    sys.stderr.write("or: size.py <starting_index> <final_index+1>\n")

pyCalculations/themis_observation.py

@@ -104,7 +104,7 @@ def simulation_to_observation_frame(x_axis,y_axis):
     z_axis = np.cross(x_axis,y_axis)
     return np.array([x_axis,y_axis,z_axis])
 
-def themis_plot_detector(vlsvReader, cellID, detector_axis=np.array([0,1,0])):
+def themis_plot_detector(vlsvReader, cellID, detector_axis=np.array([0,1,0]), pop="proton"):
     ''' Plots a view of the detector countrates using matplotlib
     :param vlsvReader:        Some VlsvReader class with a file open
     :type vlsvReader:         :class:`vlsvfile.VlsvReader`
@@ -133,7 +133,7 @@ def themis_plot_detector(vlsvReader, cellID, detector_axis=np.array([0,1,0])):
     fig.colorbar(cax)
     pl.show()
 
-def themis_plot_phasespace_contour(vlsvReader, cellID, plane_x=np.array([1.,0,0]), plane_y=np.array([0,0,1.]), smooth=False, xlabel="Vx", ylabel="Vy"):
+def themis_plot_phasespace_contour(vlsvReader, cellID, plane_x=np.array([1.,0,0]), plane_y=np.array([0,0,1.]), smooth=False, xlabel="Vx", ylabel="Vy", pop="proton"):
     ''' Plots a contour view of phasespace, as seen by a themis detector, at the given cellID
     :param vlsvReader:        Some VlsvReader class with a file open
     :type vlsvReader:         :class:`vlsvfile.VlsvReader`

pyPlots/plot_helpers.py

@@ -25,6 +25,7 @@
 import numpy as np
 import sys
 from rotation import rotateTensorToVector
+import warnings
 
 PLANE = 'XY'
 # or alternatively, 'XZ'
@@ -485,7 +486,9 @@ def expr_timeavg(pass_maps, requestvariables=False):
         var = next(listofkeys)
         if var!="dstep": break
     ntimes = len(pass_maps)
-    thismap = thesemaps[var]
+
+    warnings.warn("expr_timeavg cleaned to not produce undefined variable errors, see commit e1d2dd8ecaa7a0444ce56215f795a5237f792b1e for applied changes and check the output!")
+    thismap = pass_maps[var]
     avgmap = np.zeros(np.array(thismap.shape))
     for i in range(ntimes):
         avgmap = np.add(avgmap, pass_maps[i][var])
@@ -1134,7 +1137,7 @@ def expr_ja(pass_maps, requestvariables=False):
 
 def expr_dLstardt(pass_maps, requestvariables=False):
     if requestvariables==True:
-        return ['B','E']
+        return ['B','E','V']
 
     if type(pass_maps) is not list:
         # Not a list of time steps, calculating this value does not make sense.
@@ -1148,9 +1151,10 @@ def expr_dLstardt(pass_maps, requestvariables=False):
     thesemaps = pass_maps[curri]
     pastmaps = pass_maps[previ]
 
-    thisB = TransposeVectorArray(thesemaps['B'])
-    pastB = TransposeVectorArray(pastmaps['B'])
-    Vddt = (thisV-pastV)/DT
+    warnings.warn("expr_dLstardt cleaned to not produce undefined variable errors, see commit e1d2dd8ecaa7a0444ce56215f795a5237f792b1e for applied changes and check the output!")
+    thisV = TransposeVectorArray(thesemaps['V'])
+    pastV = TransposeVectorArray(pastmaps['V'])
+    dVdt = (thisV-pastV)/DT
 
     Bmap = TransposeVectorArray(thesemaps['B'])
     upBmag2 = np.linalg.norm(Bmap,axis=-1)**(-2)
@@ -1184,7 +1188,8 @@ def overplotvectors(ax, XmeshXY,YmeshXY, pass_maps):
     step = int(np.sqrt(colors.shape[0] * colors.shape[1]/100.))
 
     # inplane unit length vectors
-    vectmap = pt.plot.plot_helpers.inplanevec(vectmap)
+    warnings.warn("usage of inplanevec(vf) is unverified! Used to be inplanevec(vectmap), with vectmap undefined. See changes in commit e1d2dd8ecaa7a0444ce56215f795a5237f792b1e and check if results are as expected!")
+    vectmap = pt.plot.plot_helpers.inplanevec(vf)
     vectmap = vectmap / np.linalg.norm(vectmap, axis=-1)[:,:,np.newaxis]
 
     X = XmeshXY[::step,::step]
@@ -1197,8 +1202,11 @@ def overplotvectors(ax, XmeshXY,YmeshXY, pass_maps):
     elif PLANE=="YZ":
         V = vectmap[::step,::step,0]
     C = colors[::step,::step]
-    ax.quiver(X,Y,U,V,C, cmap='gray', units='dots', scale=0.03/scale, headlength=2, headwidth=2,
-                       headaxislength=2, scale_units='dots', pivot='middle')
+
+    ax.quiver(X,Y,U,V,C, cmap='gray', units='dots',
+              #scale=0.03/scale, # scale not defined - if you need to use this, adjust as necessary!
+              headlength=2, headwidth=2,
+              headaxislength=2, scale_units='dots', pivot='middle')
 
 
 def overplotstreamlines(ax, XmeshXY,YmeshXY, pass_maps):

pyPlots/plot_isosurface.py

@@ -589,14 +589,14 @@ def plot_isosurface(filename=None,
             highresscale = float(highres)
             if np.isclose(highresscale, 1.0):
                 highresscale = 2
-        figsize= [x * highresscale for x in figsize]
+        # figsize= [x * highresscale for x in figsize] # figsize not defined, neither used in this function
         fontsize=fontsize*highresscale
         fontsize2=fontsize2*highresscale
         fontsize3=fontsize3*highresscale
         scale=scale*highresscale
         thick=thick*highresscale
-        streamlinethick=streamlinethick*highresscale
-        vectorsize=vectorsize*highresscale
+        # streamlinethick=streamlinethick*highresscale # streamlinethick not defined, neither used in this function
+        # vectorsize=vectorsize*highresscale # vectorsize not defined, neither used in this function
 
     # Generate virtual bounding box to get equal aspect
     maxrange = np.array([boxcoords[1]-boxcoords[0], boxcoords[3]-boxcoords[2], boxcoords[5]-boxcoords[4]]).max() / 2.0

pyPlots/plot_vdf.py

@@ -658,7 +658,7 @@ def plot_vdf(filename=None,
                 pass
 
         if not os.access(savefigdir, os.W_OK):
-            print("No write access for directory "+outputdir2+"! Exiting.")
+            print("No write access for directory "+savefigdir+"! Exiting.")
             return
 
 

pyPlots/plot_vdf_profiles.py

@@ -434,7 +434,7 @@ def plot_vdf_profiles(filename=None,
                 pass
 
         if not os.access(savefigdir, os.W_OK):
-            print("No write access for directory "+outputdir2+"! Exiting.")
+            print("No write access for directory "+savefigdir+"! Exiting.")
             return
 
 

pyPlots/plot_vdfdiff.py

@@ -288,7 +288,7 @@ def plot_vdfdiff(filename1=None, filename2=None,
                 pass
 
         if not os.access(savefigdir, os.W_OK):
-            print("No write access for directory "+outputdir2+"! Exiting.")
+            print("No write access for directory "+savefigdir+"! Exiting.")
             return
 
 

pyVlsv/reduction.py

@@ -636,15 +636,16 @@ def thermalvelocity( variables ):
    return thermalvelocity
 
 def Vstream( variables ):
-   rhoVstream = variables[0]
-   rhostream = variables[1]
-   rhoVNonBackstream = variables[2]
-   rhoNonBackstream = variables[3]
-   # get velocity of both populations:
-   vBackstream = v( [rhoVBackstream, rhoBackstream] )
-   vNonBackstream = v( [rhoVNonBackstream, rhoNonBackstream] )
-   vBeam = vBackstream - vNonBackstream
-   return vBeam # <- is a vector quantity
+   raise NotImplementedError("rhoVBackstream, rhoBackstream not defined here. Check implementaiton if required!")
+   # rhoVstream = variables[0]
+   # rhostream = variables[1]
+   # rhoVNonBackstream = variables[2]
+   # rhoNonBackstream = variables[3]
+   # # get velocity of both populations:
+   # vBackstream = v( [rhoVBackstream, rhoBackstream] )
+   # vNonBackstream = v( [rhoVNonBackstream, rhoNonBackstream] )
+   # vBeam = vBackstream - vNonBackstream
+   # return vBeam # <- is a vector quantity
 
 def v_beam( variables ):
    vBackstream = variables[0]

pyVlsv/vlsvreader.py

@@ -3018,8 +3018,8 @@ def read_velocity_cells(self, cellid, pop="proton"):
             self.__set_cell_offset_and_blocks(pop) 
          # Check that cells has vspace
          if not cellid in self.__fileindex_for_cellid_blocks[pop]:
-            print("Cell does not have velocity distribution")
-            return []
+            warnings.warn("Cell(s) does not have velocity distribution")
+            return {}
          # Navigate to the correct position:
          offset = self.__fileindex_for_cellid_blocks[pop][cellid][0]
          num_of_blocks = self.__fileindex_for_cellid_blocks[pop][cellid][1]
@@ -3031,8 +3031,8 @@ def read_velocity_cells(self, cellid, pop="proton"):
          try:
             cells_with_blocks_index = self.__order_for_cellid_blocks[pop][cellid]
          except:
-            print("Cell does not have velocity distribution")
-            return []
+            warnings.warn("Cell(s) does not have velocity distribution")
+            return {}
          # Navigate to the correct position:
          offset = self.__blocks_per_cell_offsets[pop][cells_with_blocks_index]
          num_of_blocks = self.__blocks_per_cell[pop][cells_with_blocks_index]
@@ -3075,8 +3075,7 @@ def read_velocity_cells(self, cellid, pop="proton"):
             elif datatype == "uint" and element_size == 8:
                data_block_ids = np.fromfile(fptr, dtype = np.uint64, count = vector_size*num_of_blocks)
             else:
-               print("Error! Bad data type in blocks!")
-               return
+               raise TypeError("Error! Bad data type in blocks! datatype found was "+datatype)
 
          if (pop=="avgs") and (child.tag == "BLOCKIDS"): # Old avgs files did not have the name set for BLOCKIDS
             vector_size = ast.literal_eval(child.attrib["vectorsize"])
@@ -3092,16 +3091,15 @@ def read_velocity_cells(self, cellid, pop="proton"):
             elif datatype == "uint" and element_size == 8:
                data_block_ids = np.fromfile(fptr, dtype = np.uint64, count = vector_size*num_of_blocks)
             else:
-               print("Error! Bad data type in blocks!")
-               return
+               raise TypeError("Error! Bad data type in blocks! datatype found was "+datatype)
 
             data_block_ids = data_block_ids.reshape(num_of_blocks, vector_size)
 
       fptr.close()
 
       # Check to make sure the sizes match (just some extra debugging)
       if len(data_avgs) != len(data_block_ids):
-         print("BAD DATA SIZES")
+         raise ValueError("BAD DATA SIZES")
       # Make a dictionary (hash map) out of velocity cell ids and avgs:
       velocity_cells = {}
       array_size = len(data_avgs)

pyVlsv/vlsvwriter.py

@@ -209,7 +209,7 @@ def copy_variables_list( self, vlsvReader, vars ):
             if 'mesh' in child.attrib:
                 mesh = child.attrib['mesh']
             else:
-               if tag in ['VARIABLE']:
+               if child.tag in ['VARIABLE']:
                   print('MESH required')
                   return
                mesh = None

scripts/create_time_energy_spectrogram.py

@@ -91,6 +91,11 @@
 # bin edges of kinetic energy in electron volts (energies below and above the last and first and )
 EkinBinEdges = np.logspace(np.log10(100),np.log10(80e3),66)
 
+# Initialize as none
+xReq = None
+yReq = None
+zReq = None
+
 # give a list of cids
 cids = (4502051,4951951,5551701)
 
@@ -152,7 +157,7 @@
 cidsTemp = []
 if 'cids' not in locals():
  print('Finding nearest cells with vspace from given coordinates')
- if ('xReq' not in locals()) or ('yReq' not in locals()) or ('zReq' not in locals()):
+ if (xReq is None) or (yReq is None) or (zReq is None):
   print('ERROR: cids or (xReq,yReq,zReq) coordinates must be given')
   quit()
  if xReq.shape == yReq.shape == zReq.shape:
@@ -257,7 +262,7 @@ def doSpectra(vlsvFile):
  locs = vlsvReader.get_cellid_locations()
  cellids = list(locs.keys())
  # sort variable array according to cell ids
- locs_sorted = sorted(locs.items()), key=oper.itemgetter(0))
+ locs_sorted = sorted(locs.items(), key=oper.itemgetter(0))
  fileNameStr = os.path.basename(vlsvFile)
  spectraStr = [] # spectra file contents
  bulkStr = [] # bulk parameter file contents

scripts/magnetopause3d.py

@@ -584,10 +584,11 @@ def main():
 
 
     if plotting[2]: #plot Poynting vetors (not working!)
-        coords = np.array(coords)
-        poynting = np.array(poynting)
+        raise NotImplementedError("Poynting vector plotting not implemented!")
+        # coords = np.array(coords)
+        # poynting = np.array(poynting)
 
-        ax.quiver(coords[:,0], coords[:,1], coords[:,2], P[:,0], P[:,1], P[:,2])
+        # ax.quiver(coords[:,0], coords[:,1], coords[:,2], P[:,0], P[:,1], P[:,2])
 
 
     if plotting[3]: #plot surface

scripts/plot_jet_criteria.py

@@ -144,7 +144,7 @@ def jetcontours(ax, XmeshXY,YmeshXY, pass_maps):
                           vmin=0.8, vmax=5, 
                           external=jetcontours,
                           boxre=[8,16,-6,6],
-                          cbtitle='$n_\mathrm{p}$ [cm$^{-3}$]', 
+                          cbtitle=r'$n_\mathrm{p}$ [cm$^{-3}$]', 
                           title='', usesci=0, thick=1.2)
 
 

scripts/ulf_waves_filter.py

@@ -100,25 +100,25 @@ def ulf_filter(
         data_arr,
         "{:s}_move_ave_{:d}".format(var_to_filter, windowlength),
         units="T",
-        latex="$B_{\mathrm{ave}}$",
-        latexunits="$\mathrm{T}$",
+        latex=r"$B_{\mathrm{ave}}$",
+        latexunits=r"$\mathrm{T}$",
     )
     writer.write_variable_info(varinfo, "SpatialGrid", unitConversion=1)
     varinfo = pt.calculations.VariableInfo(
         delta_var,
         "{:s}_move_ave_{:d}_delta".format(var_to_filter, windowlength),
         units="T",
-        latex="$\delta{}B_{\mathrm{ave}}$",
-        latexunits="$\mathrm{T}$",
+        latex=r"$\delta{}B_{\mathrm{ave}}$",
+        latexunits=r"$\mathrm{T}$",
     )
     writer.write_variable_info(varinfo, "SpatialGrid", unitConversion=1)
     # window_pad =50 see above this refers to the target file (centered file)
     varinfo = pt.calculations.VariableInfo(
         filtered_dataPc2[window_pad, rev_sorti, :],
         "vg_b_vol_{:s}_{:d}_xyz".format(target_wave, windowlength),
         units="T^2Hz^-1",
-        latex="$\delta{}B_{\mathrm{ave}}$",
-        latexunits="$\mathrm{T}^2\,\mathrm{Hz}^{-1}\mathrm{orsomething}$",
+        latex=r"$\delta{}B_{\mathrm{ave}}$",
+        latexunits=r"$\mathrm{T}^2\,\mathrm{Hz}^{-1}\mathrm{orsomething}$",
     )
     writer.write_variable_info(varinfo, "SpatialGrid", unitConversion=1)
     del reader

tools/vlsvintpol.py

@@ -37,9 +37,9 @@ def extract_file(filename):
         f.optimize_open_file()
         t=f.read_parameter("time")
         if t == None:
-	    t=f.read_parameter("t")
-            if t == None:	    
-		print("Unknown time format in file " + filename)
+            t=f.read_parameter("t")
+        if t == None:	    
+            print("Unknown time format in file " + filename)
 
         for coord in coords:
             if(args.re):