Swapped cut_through for lineout in cutthrough_timeseries

analysator/calculations/cutthrough.py

@@ -27,7 +27,7 @@
 import sys
 import logging
 
-def get_cellids_coordinates_distances( vlsvReader, xmax, xmin, xcells, ymax, ymin, ycells, zmax, zmin, zcells, cell_lengths, point1, point2 ):
+def get_cellids_coordinates_distances( vlsvReader, xmax, xmin, xcells, ymax, ymin, ycells, zmax, zmin, zcells, point1, point2 ):
    ''' Calculates coordinates to be used in the cut_through. The coordinates are calculated so that every cell gets picked in the coordinates.
        :param vlsvReader:       Some open VlsvReader
        :type vlsvReader:        :class:`vlsvfile.VlsvReader`
@@ -62,8 +62,9 @@ def get_cellids_coordinates_distances( vlsvReader, xmax, xmin, xcells, ymax, ymi
       if cellid == 0:
          raise ValueError("invalid cell id!")
       # Get the max and min boundaries:
-      min_bounds = vlsvReader.get_cell_coordinates(cellid) - 0.5 * cell_lengths
-      max_bounds = np.array([min_bounds[i] + cell_lengths[i] for i in range(0,3)])
+      cell_length = vlsvReader.get_cell_dx(cellid)
+      min_bounds = vlsvReader.get_cell_coordinates(cellid) - 0.5 * cell_length
+      max_bounds = np.array([min_bounds[i] + cell_length[i] for i in range(0,3)])
       # Check which boundary we hit first when we move in the unit_vector direction:
       coefficients_min = np.divide((min_bounds - iterator), unit_vector, out=np.full(min_bounds.shape, sys.float_info.max), where=np.logical_not(unit_vector==0.0))
       coefficients_max = np.divide((max_bounds - iterator), unit_vector, out=np.full(min_bounds.shape, sys.float_info.max), where=np.logical_not(unit_vector==0.0))
@@ -107,8 +108,8 @@ def get_cellids_coordinates_distances( vlsvReader, xmax, xmin, xcells, ymax, ymi
       if min((point2 - iterator)* unit_vector) < 0:
          break
    # Return the coordinates, cellids and distances for processing
-   from output import output_1d
-   return output_1d( [np.array(cellids, copy=False), np.array(distances, copy=False), np.array(coordinates, copy=False)], ["CellID", "distances", "coordinates"], ["", "m", "m"] )
+   from analysator.calculations.output import output_1d
+   return output_1d( [np.asarray(cellids), np.asarray(distances), np.asarray(coordinates)], ["CellID", "distances", "coordinates"], ["", "m", "m"] )
 
 
 def cut_through( vlsvReader, point1, point2 ):
@@ -159,7 +160,7 @@ def cut_through( vlsvReader, point1, point2 ):
    cell_lengths = np.array([(xmax - xmin)/(float)(xcells), (ymax - ymin)/(float)(ycells), (zmax - zmin)/(float)(zcells)])
 
    # Get list of coordinates for the cells:
-   return get_cellids_coordinates_distances( vlsvReader, xmax, xmin, xcells, ymax, ymin, ycells, zmax, zmin, zcells, cell_lengths, point1, point2 )
+   return get_cellids_coordinates_distances( vlsvReader, xmax, xmin, xcells, ymax, ymin, ycells, zmax, zmin, zcells, point1, point2 )
 
 
 def cut_through_swath(vlsvReader, point1, point2, width, normal):
@@ -259,8 +260,8 @@ def cut_through_step( vlsvReader, point1, point2 ):
          break
 
    # Return the coordinates, cellids and distances for processing
-   from output import output_1d
-   return output_1d( [np.array(cellids, copy=False), np.array(distances, copy=False), np.array(coordinates, copy=False)], ["CellID", "distances", "coordinates"], ["", "m", "m"] )
+   from analysator.calculations.output import output_1d
+   return output_1d( [np.asarray(cellids), np.asarray(distances), np.asarray(coordinates)], ["CellID", "distances", "coordinates"], ["", "m", "m"] )
 
 
 # Take a curve (list of 3-coords), return cells and distances along curve as with cut_through.
@@ -318,7 +319,7 @@ def cut_through_curve(vlsvReader, curve):
        raise ValueError('point1 in cut_through_curve out of bounds')
      if vlsvReader.get_cellid(point2) == 0:
        raise ValueError('point1 in cut_through_curve out of bounds')
-     cut = get_cellids_coordinates_distances( vlsvReader, xmax, xmin, xcells, ymax, ymin, ycells, zmax, zmin, zcells, cell_lengths, point1, point2)
+     cut = get_cellids_coordinates_distances( vlsvReader, xmax, xmin, xcells, ymax, ymin, ycells, zmax, zmin, zcells, point1, point2)
      ccid = cut[0].data
      cedges = cut[1].data
      try:
@@ -355,5 +356,5 @@ def cut_through_curve(vlsvReader, curve):
      rCellIds = cellIds
      rEdges = edges
      rCoords = coords
-   from output import output_1d
-   return output_1d( [np.array(rCellIds, copy=False), np.array(rEdges, copy=False), np.array(rCoords, copy=False)], ["CellID", "distances", "coordinates"], ["", "m", "m"] )
+   from analysator.calculations.output import output_1d
+   return output_1d( [np.asarray(rCellIds), np.asarray(rEdges), np.asarray(rCoords)], ["CellID", "distances", "coordinates"], ["", "m", "m"] )

analysator/calculations/lineout.py

@@ -28,7 +28,7 @@
 import logging
 
 def lineout( vlsvReader, point1, point2, variable, operator="pass",interpolation_order=1, points=100 ):
-   ''' Returns a line cut-through from a given VLSV file for distance, coordinates and variable values. The main difference between this and cut_through is that this function interpolates a given variable.
+   ''' Returns a line cut-through from a given VLSV file for distance, coordinates and variable values. The main difference between this and cut_through is that this function interpolates a given variable to equally-spaced sample points.
 
        :param vlsvReader:            Some open VlsvReader
        :type vlsvReader:             :class:`vlsvfile.VlsvReader`
@@ -56,14 +56,24 @@ def lineout( vlsvReader, point1, point2, variable, operator="pass",interpolation
    # Transform point1 and point2 into numpy array:
    point1 = np.array(point1)
    point2 = np.array(point2)
-   # Get parameters from the file to determine a good length between points (step length):
+
+   if (interpolation_order == 0):
+      logging.warn("Lineout called with interpolation order of 0. Consider using cut_through instead, as lineout cannot guarantee complete sampling of cells along the line.")
 
    # Make sure point1 and point2 are inside bounds
    if vlsvReader.get_cellid(point1) == 0:
       raise ValueError("point1 in lineout out of bounds!")
    if vlsvReader.get_cellid(point2) == 0:
       raise ValueError("point1 in lineout out of bounds!")
 
+   # Get parameters from the file to determine a sufficient number of points, warn if not enough:
+   cell_min_lengths = vlsvReader.get_spatial_mesh_min_cell_length()
+   line_length_in_cells = np.int64(np.abs((point2-point1)/cell_min_lengths))
+
+   if (points < np.max(line_length_in_cells)):
+      logging.warn("Lineout called with {:d} points over a span of (X:{:d},Y:{:d},Z:{:d}) cells at the highest resolution.\nSuggest using at least {:d} points to sample this span.".format(points, *line_length_in_cells, np.max(line_length_in_cells))
+      )
+
    value_len=len(np.atleast_1d(vlsvReader.read_interpolated_variable( variable, point1, operator)))
 
    if value_len==1:

analysator/calculations/output.py

@@ -47,7 +47,7 @@ def output_1d( arrays, names, units="" ):
    if( (len(arrays) != len(names)) or (len(arrays) != len(units)) ):
       raise ValueError("Bad array and name length in output_1d (pyCalculations/output.py)")
    new_format = []
-   from variable import VariableInfo
+   from analysator.calculations.variable import VariableInfo
    for i in range(len(arrays)):
       variable = VariableInfo(arrays[i])
       variable.name = names[i]

analysator/vlsvfile/vlsvcache.py

@@ -104,7 +104,7 @@ def __init__(self, reader) -> None:
          logging.info("Rtree found, but tread carefully - the file caching is somewhat unstable")
          __has_rtree = True
       except:
-         logging.info("No Rtree found")
+         logging.debug("No Rtree found")
          __has_rtree = False
       if __has_rtree:
          pass

analysator/vlsvfile/vlsvreader.py

@@ -3999,6 +3999,16 @@ def get_spatial_mesh_extent(self):
       '''
       return np.array([self.__xmin, self.__ymin, self.__zmin, self.__xmax, self.__ymax, self.__zmax])
 
+   def get_spatial_mesh_min_cell_length(self):
+      ''' Calculate the side length of the smallest spatial mesh cell
+      :returns: minimum existing SpatialMesh cell sise [dx_min,dy_min,dz_min]
+      '''
+      dxs = np.array([self.__dx,self.__dy,self.__dz])
+      amr = self.get_max_refinement_level()
+
+      return dxs/2**amr
+
+
    def get_fsgrid_mesh_size(self):
       ''' Read fsgrid mesh size
 

scripts/cutthrough_timeseries.py

@@ -5,6 +5,7 @@
 import matplotlib.pyplot as plt
 import argparse
 from analysator.calculations.cutthrough import cut_through
+from analysator.calculations.lineout import lineout
 
 r_e = 6.371e6
 
@@ -13,7 +14,7 @@
 
     This script takes 12 parameters, example usage:
 
-    python cutthrough_timeseries.py -var <var> -fnr <fnr1> <fnr2> -bulkpath <bulkpath> -bulkprefix <bulkprefix> -point <point1_x> <point1_y> <point1_z> <point2_x> <point2_y> <point2_z> -outputname <outputname> -outputdir <outputdir> -filt <filt> -op <op> -cmap <cmap>
+    python cutthrough_timeseries.py -var <var> -fnr <fnr1> <fnr2> -bulkpath <bulkpath> -bulkprefix <bulkprefix> -point <point1_x> <point1_y> <point1_z> <point2_x> <point2_y> <point2_z> -outputname <outputname> -outputdir <outputdir> -filt <filt> -op <op> -cmap <cmap> -re <re> -npoints <npoints> -interpolation_order <interpolation_order
 
 
     Parameter descriptions:
@@ -29,6 +30,9 @@
         filt: Filter out temporally slowly changing signal? (<=0: no filtering, >0: filter with specified window size), default=-1
         op: Variable operator, default="pass"
         cmap: Colormap, default="viridis"
+        re: Input points given in Earth radii? default=False
+        npoints: Number of points in line, default=100
+        interpolation_order: Order of interpolation (0 or 1), default=1
 """
 
 
@@ -45,6 +49,10 @@ def jplots(
     filt=-1,
     op="pass",
     cmap="viridis",
+    re=False,
+    npoints=100,
+    interpolation_order=1,
+    draw=True,
 ):
 
     fnr_arr = np.arange(fnr1, fnr2 + 0.1, 1, dtype=int)
@@ -56,9 +64,23 @@ def jplots(
     fobj = pt.vlsvfile.VlsvReader(
         bulkpath + bulkprefix + ".{}.vlsv".format(str(fnr1).zfill(7))
     )
-    cut = cut_through(fobj, point1, point2)
-    cellids = cut[0].data[:-1]
-    distances = cut[1].data[:-1]
+    if re:
+        point1 = [point1[0] * r_e, point1[1] * r_e, point1[2] * r_e]
+        point2 = [point2[0] * r_e, point2[1] * r_e, point2[2] * r_e]
+
+    lineout0 = lineout(
+        fobj,
+        point1,
+        point2,
+        variable=var,
+        operator=op,
+        points=npoints,
+        interpolation_order=interpolation_order,
+    )
+    distances, _, _ = lineout0
+
+    fobj.optimize_clear_fileindex_for_cellid()
+
     distances = np.array(distances) / r_e
 
     data_arr = np.zeros((fnr_arr.size, distances.size), dtype=float)
@@ -69,7 +91,18 @@ def jplots(
             bulkpath + bulkprefix + ".{}.vlsv".format(str(fnr).zfill(7))
         )
         t_arr[idx] = vlsvobj.read_parameter("time")
-        data_arr[idx, :] = vlsvobj.read_variable(var, operator=op, cellids=cellids)
+
+        linecut = lineout(
+            vlsvobj,
+            point1,
+            point2,
+            variable=var,
+            operator=op,
+            points=npoints,
+            interpolation_order=interpolation_order,
+        )
+        data_arr[idx, :] = linecut[2]
+        vlsvobj.optimize_clear_fileindex_for_cellid()
     if filt > 0:
         data_arr = data_arr - uniform_filter1d(data_arr, size=filt, axis=0)
 
@@ -94,17 +127,20 @@ def jplots(
 
     cb = fig.colorbar(im, ax=ax)
 
-    if not os.path.exists(outputdir):
-        try:
-            os.makedirs(outputdir)
-        except OSError:
-            pass
+    if not draw:
+        if not os.path.exists(outputdir):
+            try:
+                os.makedirs(outputdir)
+            except OSError:
+                pass
 
-    if outputdir[-1] != "/":
-        outputdir += "/"
+        if outputdir[-1] != "/":
+            outputdir += "/"
 
-    fig.savefig(outputdir + outputname, dpi=300)
-    plt.close(fig)
+        fig.savefig(outputdir + outputname, dpi=300)
+        plt.close(fig)
+    else:
+        return (fig,ax,XmeshXY,YmeshXY,data_arr)
 
 
 def main():
@@ -148,6 +184,18 @@ def main():
     parser.add_argument(
         "-cmap", help="Colormap to use for plot", type=str, default="viridis"
     )
+    parser.add_argument(
+        "-re", help="Are input points in Earth radii?", type=bool, default=False
+    )
+    parser.add_argument(
+        "-npoints", help="Number of points in line", type=int, default=100
+    )
+    parser.add_argument(
+        "-interpolation_order",
+        help="Order of interpolation (0 or 1)",
+        type=int,
+        default=1,
+    )
     args = parser.parse_args()
 
     jplots(
@@ -163,6 +211,10 @@ def main():
         filt=args.filt,
         op=args.op,
         cmap=args.cmap,
+        re=args.re,
+        npoints=args.npoints,
+        interpolation_order=args.interpolation_order,
+        draw=False
     )