Improvements to region finding

analysator/calculations/magnetopause_sw_streamline_3d.py

@@ -245,7 +245,7 @@ def make_streamlines(vlsvfile, streamline_seeds=None, seeds_n=25, seeds_x0=20*63
     f = pt.vlsvfile.VlsvReader(file_name=vlsvfile)
 
     # Create streamline starting points if needed
-    if streamline_seeds == None:
+    if not np.any(streamline_seeds):
         streamline_seeds = np.zeros([seeds_n**2, 3])
 
         t = np.linspace(seeds_range[0], seeds_range[1], seeds_n)
@@ -301,7 +301,7 @@ def make_magnetopause(streams, end_x=-15*6371000, x_point_n=50, sector_n=36, ign
         subsolar_x = np.partition(x_axis_points[:,0], ignore)[ignore] # take the nth point as subsolar point
 
     # divide the x point numbers between x > 0 (radial) an x < 0 (yz-planes) by ratio
-    dayside_x_point_n =  int((subsolar_x/np.abs(end_x))*x_point_n)
+    dayside_x_point_n =  int((subsolar_x/np.abs(subsolar_x - end_x))*x_point_n)
 
     ### dayside magnetopause ###
     # for x > 0, look for magnetopause radially
@@ -355,11 +355,13 @@ def grid_mid_point(theta_idx, phi_idx):
 
         dayside_magnetopause[ring_idx] = ring_points
 
-
     ### x < 0 magnetopause ###
     # rest: look for magnetopause in yz-planes
     ## define points in the x axis where to find magnetopause points on the yz-plane
-    x_points = np.linspace(0.0, end_x, x_point_n-dayside_x_point_n)
+    if end_x < 0:
+        x_points = np.linspace(0.0, end_x, int((np.abs(end_x)/np.abs(subsolar_x - end_x))*x_point_n))
+    else:
+        x_points = np.array([],dtype=float)
 
     ## interpolate more exact points for streamlines at exery x_point
     new_streampoints = np.zeros((len(x_points), len(streams), 2)) # new array for keeping interpolated streamlines in form new_streampoints[x_point, streamline, y and z -coordinates] 

analysator/vlsvfile/vlsvreader.py

@@ -2491,6 +2491,19 @@ def read_variable_to_cache(self, name, operator="pass"):
       # self.__read_fileindex_for_cellid()
       return self.__variable_cache[(name,operator)]
 
+   def add_cached_variable(self, data, name, operator="pass"):
+      ''' Add a variable to cache, for the whole grid and after applying
+          operator. Works by injecting a (variable,operator) key and the data to cache.
+
+          :param data: data array - make sure it is of the correct shape and size.
+          :param name: Name of the variable (or datareducer)
+          :param operator: Datareduction operator. "pass" does no operation on data.
+
+      '''
+
+      # add data to dict, use a tuple of (name,operator) as the key [tuples are immutable and hashable]
+      self.__variable_cache[(name,operator)] = data
+
    def read_variable(self, name, cellids=-1,operator="pass"):
       ''' Read variables from the open vlsv file.
       Arguments:

analysator/vlsvfile/vlsvvtkinterface.py

@@ -78,13 +78,21 @@ def __init__(self):
 
    def SetFileName(self, filename):
       if filename != self.__FileName:
-         self.Modified()
-         self.__FileName = filename
-         if self.__FileName is not None:
-            self.__reader = pt.vlsvfile.VlsvReader(self.__FileName)
-            fn = os.path.basename(self.__FileName)
+         if filename is not None:
+            if self.__reader is None:
+               self.__reader = pt.vlsvfile.VlsvReader(self.__FileName)
+            else:
+               raise ValueError("Tried to change an existing reader ("+self.__reader.file_name+") to " + filename)
+            self.Modified()
+            self.__FileName = filename
             self.__metafile = os.path.join(self.__reader.get_cache_folder(),"vlsvvtkcache.pkl")
 
+   def SetReader(self, reader):
+      self.__reader = reader
+      self.Modified()
+      self.__FileName = reader.file_name
+      self.__metafile = os.path.join(self.__reader.get_cache_folder(),"vlsvvtkcache.pkl")
+
    def GetFileName(self):
       return self.__FileName
 

scripts/magnetopause.py

@@ -51,9 +51,16 @@
 
 
 def write_vtk_surface_to_file(vtkSurface, outfilen):
-    writer = vtk.vtkXMLPolyDataWriter() 
-    writer.SetInputConnection(vtkSurface.GetOutputPort())
+
+    writer = vtk.vtkXMLPolyDataWriter()
     writer.SetFileName(outfilen)
+
+    try:
+        writer.SetInputConnection(vtkSurface.GetOutputPort())
+    except:
+        writer.SetInputData(vtkSurface)
+
+
     writer.Write()
     logging.info("wrote ", outfilen)
 
@@ -67,7 +74,7 @@ def write_SDF_to_file(SDF, datafilen, outfilen):
 
 
 
-def magnetopause(datafilen, method="beta_star_with_connectivity", own_tresholds=None, return_surface=True, return_SDF=True, SDF_points=None, Delaunay_alpha=None, beta_star_range=[0.4, 0.5], method_args={}): # TODO: separate streamline suface and vtkDelaunay3d surface in streamline method
+def magnetopause(datafilen, method="beta_star_with_connectivity", own_tresholds=None, return_surface=True, return_SDF=True, SDF_points=None, Delaunay_alpha=None, beta_star_range=[0.0, 0.5], method_args={}): # TODO: separate streamline suface and vtkDelaunay3d surface in streamline method
     """Finds the magnetopause using the specified method. Surface is constructed using vtk's Delaunay3d triangulation which results in a convex hull if no Delaunay_alpha is given.
         Returns vtk.vtkDataSetSurfaceFilter object and/or signed distances (negative -> inside magnetopause) (=SDF) to all cells
         Note that using alpha for Delaunay might make SDF different from expected inside the magnetosphere, especially if surface is constructed with points not everywhere in the magnetosphere (e.g. beta* 0.4-0.5) or if simulation grid size is larger than alpha
@@ -79,7 +86,7 @@ def magnetopause(datafilen, method="beta_star_with_connectivity", own_tresholds=
     :kwarg return_SDF: True/False, return array of distances in m to SDF_points in point input order, negative distance inside the surface
     :kwarg SDF_points: optionally give array of own points to calculate signed distances to. If not given, distances will be to cell centres in the order of f.read_variable("CellID") output
     :kwarg Delaunay_alpha: alpha (float) to give to vtkDelaunay3d, None -> convex hull, alpha=__: surface egdes longer than __ will be excluded
-    :kwarg beta_star_range: [min, max] treshold rage to use with methods "beta_star" and "beta_star_with_connectivity"
+    :kwarg beta_star_range: [min, max] treshold range to use with methods "beta_star" and "beta_star_with_connectivity"
     :kwarg method_args: dict of keyword arguments to be passed down to external functions (for streamlines and shue)
     :returns: vtkDataSetSurfaceFilter object of convex hull or alpha shape if return_surface=True, signed distance field of convex hull or alpha shape of magnetopause if return_SDF=True
     """
@@ -118,20 +125,59 @@ def magnetopause(datafilen, method="beta_star_with_connectivity", own_tresholds=
 
     elif method == "beta_star_with_connectivity":
         # magnetopause from beta_star, with connectivity if possible
-        betastar_region = regions.treshold_mask(f.read_variable("vg_beta_star"), beta_star_range)
-        try:
-            connectivity_region = regions.treshold_mask(f.read_variable("vg_connection"), 0) # closed-closed magnetic field lines
-            magnetosphere_proper =  np.where((connectivity_region==1) | (betastar_region==1), 1, 0)
-            contour_coords = f.get_cell_coordinates(cellids[magnetosphere_proper==1])
-            np.save("pointcloud.npy", contour_coords)
-        except:
-            logging.warning("using field line connectivity for magnetosphere did not work, using only beta*")
-            #condition_dict = {"beta_star": [0.5, 0.6]} # FIC: [0.4, 0.5]) # EGE: [0.9, 1.0]) # max 0.6 in FHA to not take flyaways from outside magnetopause
-            mpause_flags = np.where(betastar_region==1, 1, 0)
-            contour_coords = f.get_cell_coordinates(cellids[mpause_flags!=0])
+        vg_beta_star =f.read_variable("vg_beta_star")
+        vg_conn = f.read_variable("vg_connection")
+        vg_classifier = vg_beta_star*np.minimum(vg_conn,1) # if closed-closed fieldlines, set var to zero
+        if True:
+
+            vtkreader = pt.vlsvfile.VlsvVtkReader()
+            vtkreader.SetReader(f)
+            f.add_cached_variable(vg_classifier, "vg_betastar_classifier")
+            vtkreader.Update()
+
+            # vars =vtkreader.findVariablesFromVlsv(getReducers=False)
+            # add here more/other datareducer outputs for downstream use if needed
+            vars = ["vg_betastar_classifier", "cellid"] 
+
+            for var in [v for v in vars if ("vg_" in v.lower()) or (v.lower() == "cellid")]:
+                vtkreader.addArrayFromVlsv(var)
+
+            vtkreader.Modified() 
+            vtkreader.Update()
+            dataport = vtkreader.GetOutputPort()
+
+            dual = vtk.vtkHyperTreeGridToDualGrid()
+            dual.SetInputConnection(dataport)
+            vtkreader.Update()
+            dual.Update()
+
+            data = dual.GetOutputDataObject(0)
+
+            threshold0 = vtk.vtkThreshold()
+            threshold0.SetInputArrayToProcess(0,0,0, vtk.vtkDataObject.FIELD_ASSOCIATION_POINTS, "vg_betastar_classifier")
+            threshold0.SetInputData(data)
+            threshold0.SetLowerThreshold(beta_star_range[0])
+            threshold0.SetUpperThreshold(beta_star_range[1])
+            threshold0.Update()
+
+            vtkSurface, SDF = regions.vtkSDF(query_points, threshold0.GetOutputDataObject(0))
+
+        else:
+            betastar_region = regions.treshold_mask(vg_classifier, beta_star_range)
+
+            try:
+                connectivity_region = regions.treshold_mask(f.read_variable("vg_connection"), 0) # closed-closed magnetic field lines
+                magnetosphere_proper =  np.where((connectivity_region==1) | (betastar_region==1), 1, 0)
+                contour_coords = f.get_cell_coordinates(cellids[magnetosphere_proper==1])
+                np.save("pointcloud.npy", contour_coords)
+            except:
+                logging.warning("using field line connectivity for magnetosphere did not work, using only beta*")
+                #condition_dict = {"beta_star": [0.5, 0.6]} # FIC: [0.4, 0.5]) # EGE: [0.9, 1.0]) # max 0.6 in FHA to not take flyaways from outside magnetopause
+                mpause_flags = np.where(betastar_region==1, 1, 0)
+                contour_coords = f.get_cell_coordinates(cellids[mpause_flags!=0])
 
-        # make a convex hull surface with vtk's Delaunay
-        vtkSurface, SDF = regions.vtkDelaunay3d_SDF(query_points, contour_coords, Delaunay_alpha)
+            # make a convex hull surface with vtk's Delaunay
+            vtkSurface, SDF = regions.vtkDelaunay3d_SDF(query_points, contour_coords, Delaunay_alpha)
 
     #elif method == "beta_star_with_fieldlines": # either incredibly slow or does not work, don't use without fixing #TODO proprer measure of actual field line backwall point
     #    # magnetopause from beta_star, with field lines connecting to back wall if possible

scripts/regions.py

@@ -67,6 +67,27 @@ def vtkDelaunay3d_SDF(query_points, coordinates, alpha=None):
     return surface, convexhull_sdf
 
 
+def vtkSDF(query_points, dualgrid):
+    ''' Obtain the SDF in relation to a waterproof volumetric grid "dualgrid"
+    '''
+
+    # print(dualgrid)
+    surface = vtk.vtkDataSetSurfaceFilter()
+    surface.SetInputData(dualgrid)
+    surface.Update()
+
+    # print(surface.GetOutput())
+
+
+    implicitPolyDataDistance = vtk.vtkImplicitPolyDataDistance()
+    implicitPolyDataDistance.SetInput(surface.GetOutput())
+    import sys
+    convexhull_sdf = np.zeros(len(query_points))
+    for i,coord in enumerate(query_points):
+        convexhull_sdf[i] = implicitPolyDataDistance.EvaluateFunction(coord)
+        # sys.exit()
+
+    return surface.GetOutputDataObject(0), convexhull_sdf
 
 
 def treshold_mask(data_array, value):
@@ -347,9 +368,9 @@ def errormsg(varstr): logging.warning("{} could not be read, will be ignored".fo
         if magnetopause_kwargs:
             __, magnetopause_SDF = magnetopause.magnetopause(datafile, **magnetopause_kwargs)
         else:
-            __, magnetopause_SDF = magnetopause.magnetopause(datafile, method="beta_star_with_connectivity", Delaunay_alpha=None)  # default magnetopause: beta*+ B connectivity convex hull 
-        write_flags(writer, magnetopause_SDF, 'SDF_magnetopause')
-        write_flags(writer, np.where(np.abs(magnetopause_SDF) < 5e6, 1, 0), "flag_magnetopause")
+            __, magnetopause_SDF = magnetopause.magnetopause(datafile, method="beta_star_with_connectivity", Delaunay_alpha=2e6, )  # default magnetopause: beta*+ B connectivity convex hull 
+        write_flags(writer, magnetopause_SDF, 'vg_SDF_magnetopause')
+        write_flags(writer, np.where(np.abs(magnetopause_SDF) < 5e6, 1, 0), "vg_flag_magnetopause")
 
         # save some magnetopause values for later
         magnetopause_density = np.mean(variables["density"][np.abs(magnetopause_SDF) < 5e6])
@@ -361,7 +382,7 @@ def errormsg(varstr): logging.warning("{} could not be read, will be ignored".fo
     # magnetosphere from magnetopause SDF
     if "magnetosphere" in regions:
         magnetosphere = np.where(magnetopause_SDF<0, 1, 0)
-        write_flags(writer, magnetosphere, 'flag_magnetosphere')
+        write_flags(writer, magnetosphere, 'vg_flag_magnetosphere')
 
 
     ## BOW SHOCK ## #TODO: similar kwargs system as magnetopause?
@@ -371,8 +392,8 @@ def errormsg(varstr): logging.warning("{} could not be read, will be ignored".fo
             bowshock = bowshock_SDF(f, variables, all_points, own_condition_dict=region_conditions["bowshock"])
         else:
             bowshock = bowshock_SDF(f, variables, all_points) # default upstream rho method, might fail with foreshock
-        write_flags(writer, bowshock, 'SDF_bowshock')
-        write_flags(writer, np.where(np.abs(bowshock) < 5e6, 1, 0), "flag_bowshock")
+        write_flags(writer, bowshock, 'vg_SDF_bowshock')
+        write_flags(writer, np.where(np.abs(bowshock) < 5e6, 1, 0), "vg_flag_bowshock")
 
         # magnetosphere+magnetosheath -area
         inside_bowshock = np.where(bowshock<0, 1, 0)
@@ -381,7 +402,7 @@ def errormsg(varstr): logging.warning("{} could not be read, will be ignored".fo
     if "magnetosheath" in regions:
         # magnetosheath from bow shock-magnetosphere difference
         magnetosheath_flags = np.where((inside_bowshock & 1-magnetosphere), 1, 0)
-        write_flags(writer, magnetosheath_flags, 'flag_magnetosheath')
+        write_flags(writer, magnetosheath_flags, 'vg_flag_magnetosheath')
 
         # save magnetosheath density and temperature for further use
         #magnetosheath_density = np.mean(variables["density"][magnetosheath_flags == 1])
@@ -392,7 +413,7 @@ def errormsg(varstr): logging.warning("{} could not be read, will be ignored".fo
     ## UPSTREAM ##
     if "upstream" in regions:
         # upstream from !bowshock
-        write_flags(writer, 1-inside_bowshock, 'flag_upstream')
+        write_flags(writer, 1-inside_bowshock, 'vg_flag_upstream')
         #write_flags(writer, inside_bowshock, 'flag_inside_bowshock')
 
 
@@ -423,7 +444,7 @@ def errormsg(varstr): logging.warning("{} could not be read, will be ignored".fo
                             }
 
         cusp_flags = make_region_flags(variables, cusp_conditions, flag_type=region_flag_type, mask=mask_inMagnetosphere)
-        write_flags(writer, cusp_flags, 'flag_cusps', mask_inMagnetosphere)
+        write_flags(writer, cusp_flags, 'vg_flag_cusps', mask_inMagnetosphere)
 
 
     # magnetotail lobes 
@@ -453,11 +474,11 @@ def errormsg(varstr): logging.warning("{} could not be read, will be ignored".fo
                                 }
             lobe_S_flags = make_region_flags(variables, lobe_S_conditions, flag_type=region_flag_type)
 
-            write_flags(writer, lobe_N_flags, 'flag_lobe_N')
-            write_flags(writer, lobe_S_flags, 'flag_lobe_S')
+            write_flags(writer, lobe_N_flags, 'vg_flag_lobe_N')
+            write_flags(writer, lobe_S_flags, 'vg_flag_lobe_S')
 
         lobes_flags = make_region_flags(variables, lobes_conditions, flag_type=region_flag_type)
-        write_flags(writer, lobes_flags, 'flag_lobes')
+        write_flags(writer, lobes_flags, 'vg_flag_lobes')
 
 
         # lobe density from median densities?
@@ -477,7 +498,7 @@ def errormsg(varstr): logging.warning("{} could not be read, will be ignored".fo
                                              }
 
         central_plasma_sheet_flags = make_region_flags(variables, central_plasma_sheet_conditions,flag_type=region_flag_type, mask=mask_inMagnetosphere)
-        write_flags(writer, central_plasma_sheet_flags, 'flag_central_plasma_sheet', mask_inMagnetosphere)
+        write_flags(writer, central_plasma_sheet_flags, 'vg_flag_central_plasma_sheet', mask_inMagnetosphere)
 
 
     ## Other boundary layers, PSBL sometimes works
@@ -495,10 +516,16 @@ def errormsg(varstr): logging.warning("{} could not be read, will be ignored".fo
 
 def main():
 
-    datafile = "/wrk-vakka/group/spacephysics/vlasiator/3D/EGE/bulk/bulk.0002000.vlsv"
-    outfilen = "EGE_regions_t2000.vlsv"
+    import sys
+
+    try:
+        fileid = int(sys.argv[1])
+    except Exception as e:
+        print("Need fileid")
+    datafile = "/wrk-vakka/group/spacephysics/vlasiator/3D/FID/bulk1/bulk1.{:07d}.vlsv".format(fileid) 
+    outfilen = "/wrk-vakka/group/spacephysics/vlasiator/3D/FID/postprocessing/magnetopause_sdfs/FID_mpause_{:07d}.vlsv".format(fileid) 
 
-    RegionFlags(datafile, outfilen, regions=["all"])
+    RegionFlags(datafile, outfilen, regions=["magnetopause"])
 
 
 if __name__ == "__main__":