Merge pull request #10 from Smithsonian/issue_09

Issue 09 and analysis

src/graspfile/analysis/grid.py

@@ -0,0 +1,79 @@
+import numpy as np
+import numpy.ma as ma
+
+
+def find_peak(field, comp=0, max_radius=None, min_radius=None):
+    """Find the peak magnitude of a component in the field.
+
+    Parameters:
+        field ``GraspField``: The field to work on.
+        comp int: The field component to look at.
+        max_radius float: Ignore portions of the grid outside this radius from the center of the field.
+        min_radius float: Ignore portions of the grid inside this radius from the center fo the field.
+
+    Returns:
+        x_peak float:, y_peak float: The x and y values of the peak value."""
+    x_vals, y_vals = field.positions_1d
+
+    f = abs(field.field[:, :, comp])
+    if max_radius is not None:
+        rad = field.radius_grid()
+        rad_max_mask = ma.masked_greater(rad, max_radius)
+        f = ma.array(f, mask=rad_max_mask.mask)
+
+    if min_radius is not None:
+        rad = field.radius_grid()
+        rad_min_mask = ma.masked_less(rad, min_radius)
+        f = ma.array(f, mask=rad_min_mask.mask)
+
+    ny, nx = np.unravel_index(np.argmax(abs(f)), f.shape)
+    x_peak = x_vals[nx]
+    y_peak = y_vals[ny]
+
+    return x_peak, y_peak
+
+
+# find the center of illumination of the field
+def find_center(field, comp=0, trunc_level=0.0, max_radius=None, min_radius=None):
+    """Find the center of illumination by finding the "center of mass" of the field.
+
+    Parameters:
+        field ``GraspField``: The field to work on.
+        comp int: The field component to look at.
+        trunc_level float: Ignore the contributions from portions of the grid below this field level.
+        max_radius float: Ignore portions of the grid outside this radius from the center of the field.
+        min_radius float: Ignore portions of the grid inside this radius from the center fo the field.
+
+    Returns:
+        x_cent float, y_cent float: The x and y values of the center of the field."""
+    xv, yv = field.positions
+
+    f = abs(field.field[:, :, comp])
+    if trunc_level != 0.0:
+        f = ma.masked_less_equal(f, trunc_level)
+        xv = ma.array(xv, mask=f.mask)
+        yv = ma.array(yv, mask=f.mask)
+
+    if max_radius is not None:
+        rad = field.radius_grid()
+        rad_max_mask = ma.masked_greater(rad, max_radius)
+        f = ma.array(f, mask=rad_max_mask.mask)
+        xv = ma.array(xv, mask=rad_max_mask.mask)
+        yv = ma.array(yv, mask=rad_max_mask.mask)
+
+    if min_radius is not None:
+        rad = field.radius_grid()
+        rad_min_mask = ma.masked_less(rad, min_radius)
+        f = ma.array(f, mask=rad_min_mask.mask)
+        xv = ma.array(xv, mask=rad_min_mask.mask)
+        yv = ma.array(yv, mask=rad_min_mask.mask)
+
+    x_illum = xv * f
+    y_illum = yv * f
+
+    norm = np.sum(f)
+
+    x_cent = np.sum(x_illum) / norm
+    y_cent = np.sum(y_illum) / norm
+
+    return x_cent, y_cent

src/graspfile/grid.py

@@ -1,6 +1,12 @@
 """This is the module for manipulating grid files containing one or more field cuts from TICRA Tools, GRASP and CHAMP
 """
 
+try:
+    import configparser
+except ImportError:
+    import ConfigParser as configparser
+
+import graspfile.numpy_utilities as numpy_utilities
 import numpy
 
 
@@ -123,8 +129,29 @@ def index_radial_dist(self, i, j):
     @property
     def positions(self):
         """Return meshed grids of the x and y positions of each point in the field"""
-        return numpy.meshgrid(numpy.linspace(self.grid_min_x, self.grid_max_x, self.grid_n_x),
-                              numpy.linspace(self.grid_min_y, self.grid_max_y, self.grid_n_y))
+        x_positions, y_positions = self.positions_1d
+        return numpy.meshgrid(x_positions, y_positions)
+
+    @property
+    def positions_1d(self):
+        """Return numpy arrays of the x and y positions used in the field"""
+        return (numpy.linspace(self.grid_min_x, self.grid_max_x, self.grid_n_x),
+                numpy.linspace(self.grid_min_y, self.grid_max_y, self.grid_n_y))
+
+    def get_value(self, xv, yv):
+        """Return the value of the field at the nearest point to xv, yv
+
+        Args:
+            xv: float containing the x coordinate of the point to get.
+            yv: float containing the y coordinate of the point to get.
+        Returns:
+            ndarray: containing self.field_components values of the field at xv, yv"""
+        x_vals, y_vals = self.positions_1d
+
+        nx = numpy_utilities.find_nearest_idx(x_vals, xv)
+        ny = numpy_utilities.find_nearest_idx(y_vals, yv)
+
+        return self.field[nx, ny, :]
 
     def radius_grid(self, center=None):
         """Return an array holding the radii of each point from the beam centre.
@@ -233,6 +260,54 @@ def __init__(self):
         self.beam_centers = []
         """list: list of beam centers for individual fields in the file."""
 
+    def parse_header(self, header):
+        """Parse the header of a grid file for useful information
+
+        Args:
+            header: list of lines in the header text.
+        """
+        # Use configparser to interpret the header info.
+        # TO-DO:
+        #   This is very dodgy, as it ignores the possibility of different frequency sets for different
+        #   sources in the file, and erase the first source's information
+        #   We should build a real parser for this that can handle multiple copies of keys
+        config = configparser.ConfigParser(strict=False, allow_no_value=True)
+        config.read_string(u"[Default]\n" + "\n".join(header))
+        config = dict(config["Default"])
+        # Parse the header to get the frequency information
+        if "frequency" in config.keys():
+            # This works for TICRA GRASP version before TICRA Tools
+            res = config["frequency"]
+            first, arg, rest = res.partition(":")
+            if first.strip() == "start_frequency":
+                # print rest
+                # We have a frequency range
+                start, stop, num_freq = rest.rsplit(",")
+                self.freqs = numpy.linspace(float(start.split()[0]), float(stop.split()[0]), int(num_freq))
+            else:
+                # We probably have a list of frequencies
+                # print res
+                freq_strs = res.rsplit("'")
+                freqs = []
+                for f in freq_strs:
+                    freqs.append(float(f.split()[0]))
+                self.freqs = numpy.array(freqs)
+        else:
+            search_key = "frequencies"
+            term = [key for key, val in config.items() if search_key in key][0]
+            value = config[term]
+
+            # This works for TICRA Tools versions > 19.0
+            #
+            # If the frequency list is long, it may spread over more than one line
+            self.freq_unit = term.strip().split()[1].strip("[]")
+
+            freq_str_list = value.split()
+            freqs = []
+            for f in freq_str_list:
+                freqs.append(float(f))
+            self.freqs = numpy.array(freqs)
+
     def read(self, fi):
         """Reads GRASP output grid files from file object and fills a number of variables
         and numpy arrays with the data"""
@@ -247,37 +322,7 @@ def read(self, fi):
             else:
                 self.header.append(line)
 
-        # Parse the header to get the frequency information
-        for (l, line) in enumerate(self.header):
-            term, arg, res = line.partition(":")
-            if term.strip() == "FREQUENCY":
-                # This works for TICRA GRASP version before TICRA Tools
-                first, arg, rest = res.partition(":")
-                if first.strip() == "start_frequency":
-                    # print rest
-                    # We have a frequency range
-                    start, stop, num_freq = rest.rsplit(",")
-                    self.freqs = numpy.linspace(float(start.split()[0]), float(stop.split()[0]), int(num_freq))
-                else:
-                    # We probably have a list of frequencies
-                    # print res
-                    freq_strs = res.rsplit("'")
-                    freqs = []
-                    for f in freq_strs:
-                        freqs.append(float(f.split()[0]))
-                    self.freqs = numpy.array(freqs)
-                break
-            elif term.strip().split()[0] == "FREQUENCIES":
-                # This works for TICRA Tools versions > 19.0
-                #
-                # If the frequency list is long, it may spread over more than one line
-                self.freq_unit = term.strip().split()[1].strip("[]")
-
-                freq_str_list = "  ".join(self.header[l+1:]).split()
-                freqs = []
-                for f in freq_str_list:
-                    freqs.append(float(f))
-                self.freqs = numpy.array(freqs)
+        self.parse_header(self.header)
 
         # We've now got through the header text and are ready to read the general
         # field type parameters

src/graspfile/numpy_utilities.py

@@ -0,0 +1,12 @@
+import numpy as np
+
+
+def find_nearest_idx(array, value):
+    """Return the index of nearest value in an array to the given value"""
+    idx = (np.abs(array-value)).argmin()
+    return idx
+
+
+def find_nearest(array, value):
+    """Return the nearest value in an array to the given value"""
+    return array[find_nearest_idx(array, value)]

src/graspfile/torfile.py

@@ -73,13 +73,10 @@ class GraspTorMember:
     """A container for the member parameter of an GraspTorObject """
 
     def __init__(self, tor_member=None):
-        #: str: The name of the member of a GraspTorObject
         self.name = None
 
-        #: str: The type of the member
         self._type = None
 
-        #: various GraspTor values: The value of the Member parameter
         self._value = None
         if tor_member:
             self.fill(tor_member)