From cab03e569e8c80b001e7686546ee714e3a8db6a7 Mon Sep 17 00:00:00 2001 From: JonasSuni Date: Thu, 5 Mar 2026 11:21:31 +0200 Subject: [PATCH 01/12] Change jplot to use interpolated variable reading --- scripts/cutthrough_timeseries.py | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/scripts/cutthrough_timeseries.py b/scripts/cutthrough_timeseries.py index 8ec4c68a..691c4ebc 100644 --- a/scripts/cutthrough_timeseries.py +++ b/scripts/cutthrough_timeseries.py @@ -61,6 +61,10 @@ def jplots( distances = cut[1].data[:-1] distances = np.array(distances) / r_e + uvector_p1_p2 = np.array( + [point2[0] - point1[0], point2[1] - point1[1], point2[2] - point1[2]] + ) / (distances[-1] * r_e) + data_arr = np.zeros((fnr_arr.size, distances.size), dtype=float) for idx in range(fnr_arr.size): @@ -69,7 +73,11 @@ 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) + for idx2 in range(distances.size): + coords = np.array(point1) + uvector_p1_p2 * distances[idx2] * r_e + data_arr[idx, idx2] = vlsvobj.read_interpolated_variable( + var, coords, operator=op + ) if filt > 0: data_arr = data_arr - uniform_filter1d(data_arr, size=filt, axis=0) From 962118a1c8d852d1bf6601d481e82ec6e953fdc4 Mon Sep 17 00:00:00 2001 From: JonasSuni Date: Thu, 5 Mar 2026 11:29:38 +0200 Subject: [PATCH 02/12] Add re option --- scripts/cutthrough_timeseries.py | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/scripts/cutthrough_timeseries.py b/scripts/cutthrough_timeseries.py index 691c4ebc..176604cb 100644 --- a/scripts/cutthrough_timeseries.py +++ b/scripts/cutthrough_timeseries.py @@ -13,7 +13,7 @@ This script takes 12 parameters, example usage: - python cutthrough_timeseries.py -var -fnr -bulkpath -bulkprefix -point -outputname -outputdir -filt -op -cmap + python cutthrough_timeseries.py -var -fnr -bulkpath -bulkprefix -point -outputname -outputdir -filt -op -cmap -re Parameter descriptions: @@ -29,6 +29,7 @@ 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 """ @@ -45,6 +46,7 @@ def jplots( filt=-1, op="pass", cmap="viridis", + re=False, ): fnr_arr = np.arange(fnr1, fnr2 + 0.1, 1, dtype=int) @@ -56,6 +58,9 @@ def jplots( fobj = pt.vlsvfile.VlsvReader( bulkpath + bulkprefix + ".{}.vlsv".format(str(fnr1).zfill(7)) ) + if re: + point1 = point1 * r_e + point2 = point2 * r_e cut = cut_through(fobj, point1, point2) cellids = cut[0].data[:-1] distances = cut[1].data[:-1] @@ -156,6 +161,9 @@ 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 + ) args = parser.parse_args() jplots( @@ -171,6 +179,7 @@ def main(): filt=args.filt, op=args.op, cmap=args.cmap, + re=args.re, ) From c3d79070b6595cd3eb5c3b3c46b126866744b1f2 Mon Sep 17 00:00:00 2001 From: JonasSuni Date: Thu, 5 Mar 2026 11:33:49 +0200 Subject: [PATCH 03/12] Fix --- scripts/cutthrough_timeseries.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/scripts/cutthrough_timeseries.py b/scripts/cutthrough_timeseries.py index 176604cb..05eec96b 100644 --- a/scripts/cutthrough_timeseries.py +++ b/scripts/cutthrough_timeseries.py @@ -59,8 +59,8 @@ def jplots( bulkpath + bulkprefix + ".{}.vlsv".format(str(fnr1).zfill(7)) ) if re: - point1 = point1 * r_e - point2 = point2 * r_e + 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] cut = cut_through(fobj, point1, point2) cellids = cut[0].data[:-1] distances = cut[1].data[:-1] From b5843e43cc40624301a45a8edc752b0b74cbafac Mon Sep 17 00:00:00 2001 From: JonasSuni Date: Thu, 5 Mar 2026 11:37:31 +0200 Subject: [PATCH 04/12] Test --- analysator/calculations/cutthrough.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/analysator/calculations/cutthrough.py b/analysator/calculations/cutthrough.py index e0563d4e..68f783d4 100644 --- a/analysator/calculations/cutthrough.py +++ b/analysator/calculations/cutthrough.py @@ -259,7 +259,7 @@ def cut_through_step( vlsvReader, point1, point2 ): break # Return the coordinates, cellids and distances for processing - from output import output_1d + from analysator.calculations.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 5244ecd0ae55a8282f8c880234bb0cb8b27772da Mon Sep 17 00:00:00 2001 From: JonasSuni Date: Thu, 5 Mar 2026 11:38:42 +0200 Subject: [PATCH 05/12] Test --- analysator/calculations/cutthrough.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/analysator/calculations/cutthrough.py b/analysator/calculations/cutthrough.py index 68f783d4..8ac109a6 100644 --- a/analysator/calculations/cutthrough.py +++ b/analysator/calculations/cutthrough.py @@ -107,7 +107,7 @@ 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 + from analysator.calculations.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 37ab1a17364de1d6f6e8d3349cbd724bb8d8b742 Mon Sep 17 00:00:00 2001 From: JonasSuni Date: Thu, 5 Mar 2026 11:41:11 +0200 Subject: [PATCH 06/12] Fix old crud in cutthrough.py --- analysator/calculations/cutthrough.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/analysator/calculations/cutthrough.py b/analysator/calculations/cutthrough.py index 8ac109a6..2c725823 100644 --- a/analysator/calculations/cutthrough.py +++ b/analysator/calculations/cutthrough.py @@ -108,7 +108,7 @@ def get_cellids_coordinates_distances( vlsvReader, xmax, xmin, xcells, ymax, ymi break # Return the coordinates, cellids and distances for processing from analysator.calculations.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"] ) + return output_1d( [np.asarray(cellids), np.asarray(distances), np.asarray(coordinates)], ["CellID", "distances", "coordinates"], ["", "m", "m"] ) def cut_through( vlsvReader, point1, point2 ): @@ -260,7 +260,7 @@ def cut_through_step( vlsvReader, point1, point2 ): # Return the coordinates, cellids and distances for processing from analysator.calculations.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"] ) + 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. @@ -355,5 +355,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"] ) From eff896d44d05bc532d1626d8325edaef4df1f14b Mon Sep 17 00:00:00 2001 From: JonasSuni Date: Thu, 5 Mar 2026 11:42:14 +0200 Subject: [PATCH 07/12] Fix old crud in output.py --- analysator/calculations/output.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/analysator/calculations/output.py b/analysator/calculations/output.py index 4c0e03cf..389476cf 100644 --- a/analysator/calculations/output.py +++ b/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] From 7f7f63f71c8983e3eebbddb3f5300b230b1310d0 Mon Sep 17 00:00:00 2001 From: Markku Alho Date: Thu, 5 Mar 2026 12:26:45 +0200 Subject: [PATCH 08/12] Initial commit, cell_lengths deprecation began --- analysator/calculations/cutthrough.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/analysator/calculations/cutthrough.py b/analysator/calculations/cutthrough.py index e0563d4e..10cd8270 100644 --- a/analysator/calculations/cutthrough.py +++ b/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)) @@ -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): @@ -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: From c8eb0c49c6dc0a189e1d6258a3430598c71e0528 Mon Sep 17 00:00:00 2001 From: JonasSuni Date: Thu, 5 Mar 2026 13:48:03 +0200 Subject: [PATCH 09/12] Swap cut_through for lineout in cutthrough_timeseries --- scripts/cutthrough_timeseries.py | 51 ++++++++++++++++++++++++-------- 1 file changed, 38 insertions(+), 13 deletions(-) diff --git a/scripts/cutthrough_timeseries.py b/scripts/cutthrough_timeseries.py index 05eec96b..d5502884 100644 --- a/scripts/cutthrough_timeseries.py +++ b/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 -fnr -bulkpath -bulkprefix -point -outputname -outputdir -filt -op -cmap -re + python cutthrough_timeseries.py -var -fnr -bulkpath -bulkprefix -point -outputname -outputdir -filt -op -cmap -re -npoints -interpolation_order 0: data_arr = data_arr - uniform_filter1d(data_arr, size=filt, axis=0) @@ -164,6 +178,15 @@ def main(): 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( @@ -180,6 +203,8 @@ def main(): op=args.op, cmap=args.cmap, re=args.re, + npoints=args.npoints, + interpolation_order=args.interpolation_order, ) From eb775d1d2944975fabaafcb8dea5715c16891d27 Mon Sep 17 00:00:00 2001 From: Markku Alho Date: Thu, 5 Mar 2026 16:36:11 +0200 Subject: [PATCH 10/12] Optimize vlsvreader cleaning for time-stepping --- scripts/cutthrough_timeseries.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/scripts/cutthrough_timeseries.py b/scripts/cutthrough_timeseries.py index d5502884..f63e88c5 100644 --- a/scripts/cutthrough_timeseries.py +++ b/scripts/cutthrough_timeseries.py @@ -76,6 +76,8 @@ def jplots( points=npoints, interpolation_order=interpolation_order, ) + fobj.optimize_clear_fileindex_for_cellid() + distances, coords, values = lineout0 distances = np.array(distances) / r_e @@ -96,6 +98,7 @@ def jplots( points=npoints, interpolation_order=interpolation_order, ) + vlsvobj.optimize_clear_fileindex_for_cellid() data_arr[idx, :] = linecut[2] if filt > 0: data_arr = data_arr - uniform_filter1d(data_arr, size=filt, axis=0) From e8c790771b23f1b577355b0fa93038354e8ac959 Mon Sep 17 00:00:00 2001 From: Markku Alho Date: Fri, 6 Mar 2026 11:47:50 +0200 Subject: [PATCH 11/12] Add sampling check to lineout to warn about undersampling, harvest unused variables to dummy. Add helper function to find smallest SpatialGrid dx/cell lengths. --- analysator/calculations/lineout.py | 14 ++++++++++++-- analysator/vlsvfile/vlsvcache.py | 2 +- analysator/vlsvfile/vlsvreader.py | 10 ++++++++++ scripts/cutthrough_timeseries.py | 6 ++++-- 4 files changed, 27 insertions(+), 5 deletions(-) diff --git a/analysator/calculations/lineout.py b/analysator/calculations/lineout.py index e34eb6c1..d32a390f 100644 --- a/analysator/calculations/lineout.py +++ b/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,7 +56,9 @@ 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: @@ -64,6 +66,14 @@ def lineout( vlsvReader, point1, point2, variable, operator="pass",interpolation 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: diff --git a/analysator/vlsvfile/vlsvcache.py b/analysator/vlsvfile/vlsvcache.py index b4b2835f..d70d7665 100644 --- a/analysator/vlsvfile/vlsvcache.py +++ b/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 diff --git a/analysator/vlsvfile/vlsvreader.py b/analysator/vlsvfile/vlsvreader.py index 717ea934..3f2f46a2 100644 --- a/analysator/vlsvfile/vlsvreader.py +++ b/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 diff --git a/scripts/cutthrough_timeseries.py b/scripts/cutthrough_timeseries.py index f63e88c5..0920d871 100644 --- a/scripts/cutthrough_timeseries.py +++ b/scripts/cutthrough_timeseries.py @@ -76,9 +76,10 @@ def jplots( points=npoints, interpolation_order=interpolation_order, ) + distances, _, _ = lineout0 + fobj.optimize_clear_fileindex_for_cellid() - distances, coords, values = lineout0 distances = np.array(distances) / r_e data_arr = np.zeros((fnr_arr.size, distances.size), dtype=float) @@ -89,6 +90,7 @@ def jplots( bulkpath + bulkprefix + ".{}.vlsv".format(str(fnr).zfill(7)) ) t_arr[idx] = vlsvobj.read_parameter("time") + linecut = lineout( vlsvobj, point1, @@ -98,8 +100,8 @@ def jplots( points=npoints, interpolation_order=interpolation_order, ) - vlsvobj.optimize_clear_fileindex_for_cellid() 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) From 5f4d3c7eaead785eb431b44097ca10af3f503ef7 Mon Sep 17 00:00:00 2001 From: Jonas Suni Date: Mon, 9 Mar 2026 09:46:39 +0200 Subject: [PATCH 12/12] Add option to jplots to return figure, axes, coordinate mesh and data instead of saving an image --- scripts/cutthrough_timeseries.py | 23 ++++++++++++++--------- 1 file changed, 14 insertions(+), 9 deletions(-) diff --git a/scripts/cutthrough_timeseries.py b/scripts/cutthrough_timeseries.py index 0920d871..ff9c98ae 100644 --- a/scripts/cutthrough_timeseries.py +++ b/scripts/cutthrough_timeseries.py @@ -52,6 +52,7 @@ def jplots( re=False, npoints=100, interpolation_order=1, + draw=True, ): fnr_arr = np.arange(fnr1, fnr2 + 0.1, 1, dtype=int) @@ -126,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(): @@ -210,6 +214,7 @@ def main(): re=args.re, npoints=args.npoints, interpolation_order=args.interpolation_order, + draw=False )