diff --git a/analysator/calculations/cutthrough.py b/analysator/calculations/cutthrough.py index e0563d4e..b5070472 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)) @@ -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"] ) 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/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] 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 8ec4c68a..ff9c98ae 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 + python cutthrough_timeseries.py -var -fnr -bulkpath -bulkprefix -point -outputname -outputdir -filt -op -cmap -re -npoints -interpolation_order 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 )