keepout_path_graphics

Implementation details of util/keepout_path_graphics.py.

keepout_path_graphics: Plot targets, planets, and keepout regions, with optional DRM overlay

For usage, use the -h option. Some options may be described there but not here.

Produces a time-stepped map, in sky coordinates, of stellar targets, planets, keepout regions, and, if a DRM is provided (see below), detections, characterizations, and slews. Output is to a movie (.mp4), or individual frames (.png). Also, some cumulative keepout summary statistics, and plots, can optionally be generated.

DRM or observing-tour display: If pickles summarizing a DRM are supplied (--drm), then the observing tour is loaded and the observations are shown. It is assumed that --drm and --spc are used together.

Typical usage:

keepout_path_graphics.py -s 0 -l 0.2 -d 0.5 -m $HOME/keepout.mp4 ./sampleScript_coron.json

where:

• -s is the start-time of the movie if zero, the start-time from the script is used. if a float in [0,1], the movie is started that proportion of the way through the mission. if a number <= 10000 is given, this is the offset from mission start in days if a number > 10000 is given, this is the MJD start time.
• -l is the length in years if zero is given (-l 0), the duration from the script is used.
• -d is the delta-t between frames in days
• -m is the name of the movie

optionally:

• -e -- use equatorial coordinates as opposed to ra/dec, HIGHLY recommended
• -C -- make a coronagraph-only movie, even if a starshade is present. (only affects movie and final frame, not cumulative keepout)
• -f DIR -- the directory name for .png frame-by-frame output
• -c DIR -- the directory name for cumulative keepout output
• --drm FILE -- the file containing the DRM as a pickle
• --spc FILE -- the file containing the SPC as a pickle

Option for experts only:

• -x SCRIPT -- the given SCRIPT filename is loaded on top of the argument script if the given SCRIPT name begins with !, it is treated as a json literal rather than a filename

GraphicsStyle

Bases: object

Singleton class that holds graphics styles.

Source code in util/keepout_path_graphics.py

class GraphicsStyle(object):
    r'''Singleton class that holds graphics styles.'''
    # mapping of conditions to colors
    # formerly: char_ok   = (0.0, 0.5, 1.0) # sky blue
    char_HZ   = (0.0, 0.7, 0.0) # green
    char_ok   = (0.5, 0.2, 1.0) # purple
    char_fail = (0.8, 0.1, 0.0) # brick-red
    det_HZ    = (0.0, 0.7, 0.0) # green
    det_ok    = (0.5, 0.2, 1.0) # purple
    det_fail  = (0.8, 0.1, 0.0) # brick-red
    # line/shape styles
    char_mec  = 'black' # characterization, marker edge color
    det_mec = None # detection, marker edge color
    char_lw = 0.5 # characterization, linewidth
    det_lw = 0.0 # detection, linewidth

    def __init__(self): pass

    def count2pt(self, count):
        r'''Display "points" (i.e., circle radius in point units) per detection.'''
        return np.sqrt(0.3 + 1.7 * count)

    def count2pt2(self, count):
        r'''Display "points" (i.e., markersize point^2 units) per detection.
        The relationship between count2pt vs. count2pt2 is empirical.'''
        return 1.7 + 1.0 * count

count2pt(count)

Display "points" (i.e., circle radius in point units) per detection.

Source code in util/keepout_path_graphics.py

def count2pt(self, count):
    r'''Display "points" (i.e., circle radius in point units) per detection.'''
    return np.sqrt(0.3 + 1.7 * count)

count2pt2(count)

Display "points" (i.e., markersize point^2 units) per detection. The relationship between count2pt vs. count2pt2 is empirical.

Source code in util/keepout_path_graphics.py

def count2pt2(self, count):
    r'''Display "points" (i.e., markersize point^2 units) per detection.
    The relationship between count2pt vs. count2pt2 is empirical.'''
    return 1.7 + 1.0 * count

ObserveInfo

Bases: object

Observations made by a mission, in order, as loaded from an external DRM pickle.

Source code in util/keepout_path_graphics.py

class ObserveInfo(object):
    r"""Observations made by a mission, in order, as loaded from an external DRM pickle."""

    def load_from_spc_file(self, loc, spc):
        r'''Load the DRM, and star/planet info from a "spc" file given as an argument.
        This drm + spc file transfer is compatible the Exosims ipyparallel output.'''
        # these will fail noisily if there is no file present
        print('Loading DRM from', loc)
        self.drm = pickle.load(open(loc, 'rb'), **PICKLE_ARGS)
        # spc file contains a dict with many fields - save them all
        self.spc = pickle.load(open(spc, 'rb'), **PICKLE_ARGS)
        # and then extract the ones we need
        self.snm = self.spc['Name']
        #self.sco = self.spc['coords'] # we have the script, so coords are not needed
        #self.spc_coord = self.sco.heliocentrictrueecliptic
        # planet Equivalent Insolation Distance (EID), for HZ test
        self.EID_planet = (self.spc['a'] / np.sqrt(self.spc['L'][self.spc['plan2star']])).value
        # planet radius, in EarthRad
        self.Rp_planet = self.spc['Rp'].value

    def __init__(self, loc, spc):

        # load DRM and Star-Planet info
        if spc:
            self.load_from_spc_file(loc, spc)
        # number of observations
        N = len(self.drm)
        # stars we visited, in DRM order
        self.drm_s_ind = [r['star_ind'] for r in self.drm]
        # visit times
        self.times = np.array([extract_time(r['arrival_time']) for r in self.drm])
        # summaries
        # Note, there is room for improvement with recording successful observations,
        # and HZ observations.  Specifically, we could record number of dets/chars/HZ,
        # or we could record only unique dets/chars/HZ.  It really depends on what we
        # want to plot.
        self.was_det  = np.zeros((N,), dtype=bool) # detection or char?
        self.success  = np.zeros((N,), dtype=bool) # was it successful for any planet?
        self.hab_zone = np.zeros((N,), dtype=bool) # were there successful HZ observations?
        for i, obs in enumerate(self.drm):
            # deal with various key possibilities
            if 'det_status0' in obs:
                Dkey = 'det_status0'
            elif 'det_status1' in obs:
                Dkey = 'det_status1'
            else:
                Dkey = 'det_status'
            # was the observation a detection?
            was_det = (Dkey in obs)
            self.was_det[i] = was_det
            # query detections/characterizations for status in the same way
            if was_det:
                obs_status = obs[Dkey]
            else:
                obs_status = get_char_status(obs)
            # 03-Nov-2021, turmon // Nov-2023, turmon, believe the below is obsolete
            # inserted due to some DRM's containing records that are neither det nor char
            ##    print('Masked error: set DRM[%d][%s] = 0 because key is absent (star_ind = %d)' %
            ##              (i, st_key, obs['star_ind']))
            # was the observation successful?
            self.success[i] = np.any(np.array(obs_status) > 0)
            # was any successful observation done on a HZ planet?
            # 1/ planet-indexes of successful observations, if any
            success_plan_inds = np.array(obs['plan_inds'],dtype=int)[np.where(obs_status > 0)[0]]
            # 2/ indicator for each successful planet, Is the planet in the HZ?
            #    turmon 11/2021: changed the upper boundary of L_planet from 1.55 to 1.67
            hab_zone = np.logical_and(self.EID_planet[success_plan_inds] >= 0.40,
                                      self.EID_planet[success_plan_inds] <= 1.67)
            # 3/ was any successful planet observation in the HZ?
            #    Note: perhaps we actually want a count?
            self.hab_zone[i] = np.any(hab_zone)

    def tour_summary(self, args):
        '''Compute and optionally dump tour summary statistics.'''
        Nstar = len(self.snm)
        tour = [obs for i, obs in enumerate(self.drm) if not self.was_det[i]]
        # list of all stars visited
        self.visited = np.array([d['star_ind'] for d in tour], dtype=int)
        # total cumulative number-of-visits, by star
        self.vcount = np.bincount(self.visited, minlength=Nstar)
        # slew adjacency matrix
        self.visit2 = np.zeros((Nstar, Nstar), dtype=int)
        # slew transition list (on what chars did we go from d1->d2)
        visit_when = np.empty((Nstar, Nstar), dtype=object)
        for d1 in range(Nstar):
            for d2 in range(Nstar):
                visit_when[d1,d2] = []
        # if only one thing in tour, there was no slew
        if len(tour) > 1:
            for when, d1, d2 in zip(list(range(len(tour)-1)), tour[0:-1], tour[1:]):
                self.visit2[d1['star_ind'], d2['star_ind']] += 1
                visit_when[ d1['star_ind'], d2['star_ind']].append(when)

        # dump some info - should use a csv-writer
        if args.out_cume:
            # visit counts
            fn = os.path.join(args.out_cume, 'path-visits.csv')
            with open(fn, 'w') as f:
                arr = ['name', 'visit_mean', 'lon', 'lat', 'dist']
                f.write(','.join(arr) + '\n')
                for i, ct in enumerate(self.vcount):
                    arr = [
                        '%s' % self.snm[i],
                        '%.1f' % ct,
                        '%.2f' % self.xpos[i],  # [deg]
                        '%.2f' % self.ypos[i],  # [deg]
                        '%.2f' % self.distance[i], # [pc]
                        ]
                    f.write(','.join(arr) + '\n')
                print("Visits written to `%s'" % fn)
            # slew counts
            fn = os.path.join(args.out_cume, 'path-slews.csv')
            with open(fn, 'w') as f:
                pos1, pos2 = np.where(self.visit2 > 0)
                arr = ['source', 'dest', 'slews', 'label']
                f.write(','.join(arr) + '\n')
                for p1, p2 in zip(pos1, pos2):
                    # e.g., "Slew Numbers: 17; 39"
                    label_string = ('"Slew Number%s: %s"' % (
                        '' if len(visit_when[p1,p2]) == 1 else "s",
                        '; '.join(str(when+1) for when in visit_when[p1,p2])))
                    arr = [
                        '%d' % p1,
                        '%d' % p2,
                        '%d' % self.visit2[p1,p2],
                        label_string
                        ]
                    f.write(','.join(arr) + '\n')
                print("Slews written to `%s'" % fn)

    def summary(self):
        s = f'Loaded DRM:\n  {len(self.drm)} observations\n  {len(set(self.drm_s_ind))} stars observed\n  {len(self.snm)} stars total'
        return s

    def lookup(self, day):
        r"""Return index of observation closest to day."""
        return np.argmin(np.abs(self.times - day))

    def obs_history(self, day):
        r"""Return observation history up to the given day: counts of detections and charaterizations.
        Each observation fits into one of the cases below.
        Coded as below, where each is a Counter indexed by star number. 
          dets0 = observation with no detections.
          dets1 = observation with >=1 detection.
          chars0 = characterization attempt, failed.
          chars1 = characterization attempt, success.
        Could enhance to add another counter, for habitable-zone detections vs. non-HZ detections.
        """
        index = self.lookup(day)
        # TODO: prefer a single dict of Counter(), instead of multiple Counter's
        detsHZ  = Counter(); dets1  = Counter(); dets0  = Counter()
        charsHZ = Counter(); chars1 = Counter(); chars0 = Counter()
        for i in range(index):
            # star index for this observation
            star_ind = self.drm[i]['star_ind']
            # FIXME: Some DRMs have chars and dets in the same observation
            if self.was_det[i]:
                # detection
                # FIXME: HZ detections will wipe out non-HZ detections
                # could fix this with a count of HZ and non-HZ detections
                if self.hab_zone[i]:
                    detsHZ[star_ind] += 1
                elif self.success[i]:
                    dets1[star_ind] += 1
                else:
                    dets0[star_ind] += 1
            else:
                # characterization
                if self.hab_zone[i]:
                    charsHZ[star_ind] += 1
                elif self.success[i]:
                    chars1[star_ind] += 1
                else:
                    chars0[star_ind] += 1
        return detsHZ, dets1, dets0, charsHZ, chars1, chars0

    def obs_summary(self, day, normal, z0):
        r"""Return a summary of observation history up to given day.

        If normal is False, it is the last day, and slightly different summary rules hold.
        The summary is a list of tuples, intended to be: star_number, color, size, extra.
        In this case, `extra' is a dictionary of point properties, which will be given 
        to the matplotlib circle-plotter.
        To control plot order, we need to have the base z-order (z0) of the returned
        summaries.  When we need objects to come out on top, we use z0+(dz), as follows:
          dets: z0; dets-on-top: z0+1; chars: z0+2; chars-on-top: z0:3
        Because of closely-packed binaries, it is important to put chars on top.
        """
        # observations to emphasize:
        #  set of star indexes in the present slew
        stars_important = set(self.stars_on_slew(day))
        # obtain graphics styles
        GS = GraphicsStyle()

        # characterization marker style - always on top of dets
        c_style = dict(ec=GS.char_mec, lw=GS.char_lw, zorder=z0+2)
        # make rv, the list of observations-to-date
        rv = []
        detsHZ, dets1, dets0, charsHZ, chars1, chars0 = self.obs_history(day)
        for star in (detsHZ + dets1 + dets0 + charsHZ + chars1 + chars0).keys():
            # if not "important", detected stars will be dimmer
            # if last frame (not normal), no detected stars will be dimmed
            # in either case, set no edge: the edge overlaps the body of the circle,
            # and then, when made transparent, donut shapes result.
            if (star in stars_important) or (not normal):
                d_style = dict(ec=GS.det_mec, lw=GS.det_lw, zorder=z0)
            else:
                d_style = dict(ec=GS.det_mec, lw=GS.det_lw, zorder=z0, alpha=0.3)
            # order of *if* statements encodes display preference order - show only one glyph per star
            if chars1[star] > 0 and charsHZ[star] > 0:
                # HZ and regular: circle within a circle. Omit fails.
                # count2pt(1+...) gives extra room for the bounding circle
                rv.append( (star, GS.char_HZ, GS.count2pt(1+chars1[star]+charsHZ[star]), c_style) )
                rv.append( (star, GS.char_ok, GS.count2pt(chars1[star]), dict(c_style, zorder=z0+3)) )
            elif charsHZ[star] > 0:
                rv.append( (star, GS.char_HZ, GS.count2pt(charsHZ[star]), c_style) )
            elif chars1[star] > 0:
                rv.append( (star, GS.char_ok, GS.count2pt(chars1[star]),  c_style) )
            elif chars0[star] > 0:
                rv.append( (star, GS.char_fail, GS.count2pt(chars0[star]), c_style) )
            elif dets1[star] > 0 and detsHZ[star] > 0:
                # HZ and regular: circle within a circle. Omit fails.
                rv.append( (star, GS.det_HZ, GS.count2pt(1+dets1[star]+detsHZ[star]), d_style) )
                rv.append( (star, GS.det_ok, GS.count2pt(dets1[star]), dict(d_style, zorder=z0+1)) )
            elif detsHZ[star] > 0:
                rv.append( (star, GS.det_HZ, GS.count2pt(detsHZ[star]), d_style) )
            elif dets1[star] > 0:
                rv.append( (star, GS.det_ok, GS.count2pt(dets1[star]), d_style) )
            elif dets0[star] > 0:
                # all detections failed
                rv.append( (star, GS.det_fail, GS.count2pt(dets0[star]), d_style) )
            else:
                assert False, 'Condition should be unreachable'
        return rv

    def char_bracket(self, day):
        r"""Return pair of characterization observations bracketing a given day."""
        center = self.lookup(day)
        hi = np.min(np.where(self.was_det[center:] == False)[0])
        # char_status is not at top-level of all the DRM's
        # lo = max([i for i in range(center) if 'char_status' in self.drm[i]])
        lo = max([i for i in range(center) if not self.was_det[i]])
        if not hi or not lo:
            return None, None
        else:
            return lo, hi

    def char_before(self, day):
        r'''Return index of the previous characterization BEFORE day'''
        center = self.lookup(day)
        try:
            endpoint = np.where(self.was_det[0:center] == False)[0][-1]
        except IndexError:
            endpoint = 0
        return endpoint

    def char_after(self, day):
        r'''Return index of the next characterization AFTER day'''
        center = self.lookup(day)
        try:
            endpoint = center + np.where(self.was_det[center:] == False)[0][0]
        except IndexError:
            endpoint = len(self.was_det) - 1
        return endpoint

    def char_to_date(self, day):
        r"""Return all characterization observations up to a given day."""
        # catch the next characterization AFTER day
        endpoint = self.char_after(day)
        # return DRM entries from 0:endpoint-1 that were NOT detections
        return np.take(self.drm, np.where(self.was_det[:endpoint+1] == False)[0])

    def stars_on_slew(self, day):
        r"""Return list of all stars in the characterization-slew-interval containing day."""
        i0 = self.char_before(day)
        i1 = self.char_after(day)
        return self.drm_s_ind[i0:i1+1]

char_after(day)

Return index of the next characterization AFTER day

Source code in util/keepout_path_graphics.py

def char_after(self, day):
    r'''Return index of the next characterization AFTER day'''
    center = self.lookup(day)
    try:
        endpoint = center + np.where(self.was_det[center:] == False)[0][0]
    except IndexError:
        endpoint = len(self.was_det) - 1
    return endpoint

char_before(day)

Return index of the previous characterization BEFORE day

Source code in util/keepout_path_graphics.py

def char_before(self, day):
    r'''Return index of the previous characterization BEFORE day'''
    center = self.lookup(day)
    try:
        endpoint = np.where(self.was_det[0:center] == False)[0][-1]
    except IndexError:
        endpoint = 0
    return endpoint

char_bracket(day)

Return pair of characterization observations bracketing a given day.

Source code in util/keepout_path_graphics.py

def char_bracket(self, day):
    r"""Return pair of characterization observations bracketing a given day."""
    center = self.lookup(day)
    hi = np.min(np.where(self.was_det[center:] == False)[0])
    # char_status is not at top-level of all the DRM's
    # lo = max([i for i in range(center) if 'char_status' in self.drm[i]])
    lo = max([i for i in range(center) if not self.was_det[i]])
    if not hi or not lo:
        return None, None
    else:
        return lo, hi

char_to_date(day)

Return all characterization observations up to a given day.

Source code in util/keepout_path_graphics.py

def char_to_date(self, day):
    r"""Return all characterization observations up to a given day."""
    # catch the next characterization AFTER day
    endpoint = self.char_after(day)
    # return DRM entries from 0:endpoint-1 that were NOT detections
    return np.take(self.drm, np.where(self.was_det[:endpoint+1] == False)[0])

load_from_spc_file(loc, spc)

Load the DRM, and star/planet info from a "spc" file given as an argument. This drm + spc file transfer is compatible the Exosims ipyparallel output.

Source code in util/keepout_path_graphics.py

def load_from_spc_file(self, loc, spc):
    r'''Load the DRM, and star/planet info from a "spc" file given as an argument.
    This drm + spc file transfer is compatible the Exosims ipyparallel output.'''
    # these will fail noisily if there is no file present
    print('Loading DRM from', loc)
    self.drm = pickle.load(open(loc, 'rb'), **PICKLE_ARGS)
    # spc file contains a dict with many fields - save them all
    self.spc = pickle.load(open(spc, 'rb'), **PICKLE_ARGS)
    # and then extract the ones we need
    self.snm = self.spc['Name']
    #self.sco = self.spc['coords'] # we have the script, so coords are not needed
    #self.spc_coord = self.sco.heliocentrictrueecliptic
    # planet Equivalent Insolation Distance (EID), for HZ test
    self.EID_planet = (self.spc['a'] / np.sqrt(self.spc['L'][self.spc['plan2star']])).value
    # planet radius, in EarthRad
    self.Rp_planet = self.spc['Rp'].value

lookup(day)

Return index of observation closest to day.

Source code in util/keepout_path_graphics.py

def lookup(self, day):
    r"""Return index of observation closest to day."""
    return np.argmin(np.abs(self.times - day))

obs_history(day)

Return observation history up to the given day: counts of detections and charaterizations. Each observation fits into one of the cases below. Coded as below, where each is a Counter indexed by star number. dets0 = observation with no detections. dets1 = observation with >=1 detection. chars0 = characterization attempt, failed. chars1 = characterization attempt, success. Could enhance to add another counter, for habitable-zone detections vs. non-HZ detections.

Source code in util/keepout_path_graphics.py

def obs_history(self, day):
    r"""Return observation history up to the given day: counts of detections and charaterizations.
    Each observation fits into one of the cases below.
    Coded as below, where each is a Counter indexed by star number. 
      dets0 = observation with no detections.
      dets1 = observation with >=1 detection.
      chars0 = characterization attempt, failed.
      chars1 = characterization attempt, success.
    Could enhance to add another counter, for habitable-zone detections vs. non-HZ detections.
    """
    index = self.lookup(day)
    # TODO: prefer a single dict of Counter(), instead of multiple Counter's
    detsHZ  = Counter(); dets1  = Counter(); dets0  = Counter()
    charsHZ = Counter(); chars1 = Counter(); chars0 = Counter()
    for i in range(index):
        # star index for this observation
        star_ind = self.drm[i]['star_ind']
        # FIXME: Some DRMs have chars and dets in the same observation
        if self.was_det[i]:
            # detection
            # FIXME: HZ detections will wipe out non-HZ detections
            # could fix this with a count of HZ and non-HZ detections
            if self.hab_zone[i]:
                detsHZ[star_ind] += 1
            elif self.success[i]:
                dets1[star_ind] += 1
            else:
                dets0[star_ind] += 1
        else:
            # characterization
            if self.hab_zone[i]:
                charsHZ[star_ind] += 1
            elif self.success[i]:
                chars1[star_ind] += 1
            else:
                chars0[star_ind] += 1
    return detsHZ, dets1, dets0, charsHZ, chars1, chars0

obs_summary(day, normal, z0)

Return a summary of observation history up to given day.

If normal is False, it is the last day, and slightly different summary rules hold. The summary is a list of tuples, intended to be: star_number, color, size, extra. In this case, `extra' is a dictionary of point properties, which will be given to the matplotlib circle-plotter. To control plot order, we need to have the base z-order (z0) of the returned summaries. When we need objects to come out on top, we use z0+(dz), as follows: dets: z0; dets-on-top: z0+1; chars: z0+2; chars-on-top: z0:3 Because of closely-packed binaries, it is important to put chars on top.

Source code in util/keepout_path_graphics.py

def obs_summary(self, day, normal, z0):
    r"""Return a summary of observation history up to given day.

    If normal is False, it is the last day, and slightly different summary rules hold.
    The summary is a list of tuples, intended to be: star_number, color, size, extra.
    In this case, `extra' is a dictionary of point properties, which will be given 
    to the matplotlib circle-plotter.
    To control plot order, we need to have the base z-order (z0) of the returned
    summaries.  When we need objects to come out on top, we use z0+(dz), as follows:
      dets: z0; dets-on-top: z0+1; chars: z0+2; chars-on-top: z0:3
    Because of closely-packed binaries, it is important to put chars on top.
    """
    # observations to emphasize:
    #  set of star indexes in the present slew
    stars_important = set(self.stars_on_slew(day))
    # obtain graphics styles
    GS = GraphicsStyle()

    # characterization marker style - always on top of dets
    c_style = dict(ec=GS.char_mec, lw=GS.char_lw, zorder=z0+2)
    # make rv, the list of observations-to-date
    rv = []
    detsHZ, dets1, dets0, charsHZ, chars1, chars0 = self.obs_history(day)
    for star in (detsHZ + dets1 + dets0 + charsHZ + chars1 + chars0).keys():
        # if not "important", detected stars will be dimmer
        # if last frame (not normal), no detected stars will be dimmed
        # in either case, set no edge: the edge overlaps the body of the circle,
        # and then, when made transparent, donut shapes result.
        if (star in stars_important) or (not normal):
            d_style = dict(ec=GS.det_mec, lw=GS.det_lw, zorder=z0)
        else:
            d_style = dict(ec=GS.det_mec, lw=GS.det_lw, zorder=z0, alpha=0.3)
        # order of *if* statements encodes display preference order - show only one glyph per star
        if chars1[star] > 0 and charsHZ[star] > 0:
            # HZ and regular: circle within a circle. Omit fails.
            # count2pt(1+...) gives extra room for the bounding circle
            rv.append( (star, GS.char_HZ, GS.count2pt(1+chars1[star]+charsHZ[star]), c_style) )
            rv.append( (star, GS.char_ok, GS.count2pt(chars1[star]), dict(c_style, zorder=z0+3)) )
        elif charsHZ[star] > 0:
            rv.append( (star, GS.char_HZ, GS.count2pt(charsHZ[star]), c_style) )
        elif chars1[star] > 0:
            rv.append( (star, GS.char_ok, GS.count2pt(chars1[star]),  c_style) )
        elif chars0[star] > 0:
            rv.append( (star, GS.char_fail, GS.count2pt(chars0[star]), c_style) )
        elif dets1[star] > 0 and detsHZ[star] > 0:
            # HZ and regular: circle within a circle. Omit fails.
            rv.append( (star, GS.det_HZ, GS.count2pt(1+dets1[star]+detsHZ[star]), d_style) )
            rv.append( (star, GS.det_ok, GS.count2pt(dets1[star]), dict(d_style, zorder=z0+1)) )
        elif detsHZ[star] > 0:
            rv.append( (star, GS.det_HZ, GS.count2pt(detsHZ[star]), d_style) )
        elif dets1[star] > 0:
            rv.append( (star, GS.det_ok, GS.count2pt(dets1[star]), d_style) )
        elif dets0[star] > 0:
            # all detections failed
            rv.append( (star, GS.det_fail, GS.count2pt(dets0[star]), d_style) )
        else:
            assert False, 'Condition should be unreachable'
    return rv

stars_on_slew(day)

Return list of all stars in the characterization-slew-interval containing day.

Source code in util/keepout_path_graphics.py

def stars_on_slew(self, day):
    r"""Return list of all stars in the characterization-slew-interval containing day."""
    i0 = self.char_before(day)
    i1 = self.char_after(day)
    return self.drm_s_ind[i0:i1+1]

tour_summary(args)

Compute and optionally dump tour summary statistics.

Source code in util/keepout_path_graphics.py

def tour_summary(self, args):
    '''Compute and optionally dump tour summary statistics.'''
    Nstar = len(self.snm)
    tour = [obs for i, obs in enumerate(self.drm) if not self.was_det[i]]
    # list of all stars visited
    self.visited = np.array([d['star_ind'] for d in tour], dtype=int)
    # total cumulative number-of-visits, by star
    self.vcount = np.bincount(self.visited, minlength=Nstar)
    # slew adjacency matrix
    self.visit2 = np.zeros((Nstar, Nstar), dtype=int)
    # slew transition list (on what chars did we go from d1->d2)
    visit_when = np.empty((Nstar, Nstar), dtype=object)
    for d1 in range(Nstar):
        for d2 in range(Nstar):
            visit_when[d1,d2] = []
    # if only one thing in tour, there was no slew
    if len(tour) > 1:
        for when, d1, d2 in zip(list(range(len(tour)-1)), tour[0:-1], tour[1:]):
            self.visit2[d1['star_ind'], d2['star_ind']] += 1
            visit_when[ d1['star_ind'], d2['star_ind']].append(when)

    # dump some info - should use a csv-writer
    if args.out_cume:
        # visit counts
        fn = os.path.join(args.out_cume, 'path-visits.csv')
        with open(fn, 'w') as f:
            arr = ['name', 'visit_mean', 'lon', 'lat', 'dist']
            f.write(','.join(arr) + '\n')
            for i, ct in enumerate(self.vcount):
                arr = [
                    '%s' % self.snm[i],
                    '%.1f' % ct,
                    '%.2f' % self.xpos[i],  # [deg]
                    '%.2f' % self.ypos[i],  # [deg]
                    '%.2f' % self.distance[i], # [pc]
                    ]
                f.write(','.join(arr) + '\n')
            print("Visits written to `%s'" % fn)
        # slew counts
        fn = os.path.join(args.out_cume, 'path-slews.csv')
        with open(fn, 'w') as f:
            pos1, pos2 = np.where(self.visit2 > 0)
            arr = ['source', 'dest', 'slews', 'label']
            f.write(','.join(arr) + '\n')
            for p1, p2 in zip(pos1, pos2):
                # e.g., "Slew Numbers: 17; 39"
                label_string = ('"Slew Number%s: %s"' % (
                    '' if len(visit_when[p1,p2]) == 1 else "s",
                    '; '.join(str(when+1) for when in visit_when[p1,p2])))
                arr = [
                    '%d' % p1,
                    '%d' % p2,
                    '%d' % self.visit2[p1,p2],
                    label_string
                    ]
                f.write(','.join(arr) + '\n')
            print("Slews written to `%s'" % fn)

circle_marker(**given_props)

Return a Line2D with circular markers for an axis legend, given_props over-rides. Some key properties to pass in: markeredgecolor (mec), markeredgewidth (mew), markerfacecolor (mfc), markersize.

Source code in util/keepout_path_graphics.py

def circle_marker(**given_props):
    r'''Return a Line2D with circular markers for an axis legend, given_props over-rides.
    Some key properties to pass in: markeredgecolor (mec), markeredgewidth (mew), 
    markerfacecolor (mfc), markersize.'''
    # x, y are not used when in the legend
    # lw=0 => turn off lines between markers
    # markers set to circles
    # mew=0 => turn off marker edges by default
    props = dict(xdata=[0], ydata=[0], lw=0, marker='o', markeredgewidth=0)
    props.update(given_props)
    return plt.Line2D(**props)

ensure_dir(directory)

Ensure enclosing dir exists, especially in multiprocessing context.

Source code in util/keepout_path_graphics.py

def ensure_dir(directory):
    r'''Ensure enclosing dir exists, especially in multiprocessing context.'''
    # the loop guards against another process creating the dir after we checked,
    # but before the os.makedirs() call.
    tries, done = 0, False
    while not os.path.exists(directory) and not done:
        try:
            os.makedirs(directory, 0o775)
            done = True # ensure we exit the while
        except OSError:
            # hope it was just concurrent creation
            tries += 1
            time.sleep(random.uniform(0.1, 0.2))
            if tries == 2:
                raise # give up

extract_time(tm)

Unbox a time, if it is boxed, for compatibility with multiple DRMs.

Source code in util/keepout_path_graphics.py

def extract_time(tm):
    r'''Unbox a time, if it is boxed, for compatibility with multiple DRMs.'''
    try:
        return tm.value
    except AttributeError:
        return tm

gcp(lon1, lat1, lon2, lat2)

GCP = great circle points from (lon1,lat1) -> (lon2,lat2).

Point-count is auto-determined for smoothness without waste.

Source code in util/keepout_path_graphics.py

def gcp(lon1, lat1, lon2, lat2):
    r"""GCP = great circle points from (lon1,lat1) -> (lon2,lat2).

    Point-count is auto-determined for smoothness without waste."""
    Npt = 200
    p1 = np.array(lonlat2xyz(lon1, lat1))
    p2 = np.array(lonlat2xyz(lon2, lat2))
    delta = p2 - p1
    # length of the line: between 0 and 2, generally < 1
    r = np.linalg.norm(delta)
    d = np.linspace(0, 1.0, int(Npt*r + 10))
    # sample points along the line connecting p1 and p2
    pts = p1 + np.outer(d, delta)
    # convert these points to lon,lat
    pts_ll = [xyz2lonlat(p[0], p[1], p[2]) for p in pts]
    return pts_ll

get_arrow(lines, length)

"Pick one of the line segments within lines, and make an arrow point along it.

Source code in util/keepout_path_graphics.py

def get_arrow(lines, length):
    r""""Pick one of the line segments within lines, and make an arrow point along it."""
    n0 = len(lines) - 1
    arrow_0 = np.array(lines[n0][0])
    arrow_1 = np.array(lines[n0][1])
    delta = arrow_1 - arrow_0
    # make it unit norm - safeguard in case line is degenerate
    delta = delta / np.maximum(1e-6, np.linalg.norm(delta))
    return arrow_0, length * delta

get_char_status(obs)

Utility function, gets char status from a drm observation.

Source code in util/keepout_path_graphics.py

def get_char_status(obs):
    r'''Utility function, gets char status from a drm observation.
    '''
    if 'char_status' in obs:
        # starshade char - has a unitary char_status key
        return obs['char_status']
    else:
        # coronagraph-only - get per-planet char status by combining channels
        # Each char status is a vector of -1, 0, +1 for each planet.
        # Planets are combined separately: 2 planets -> return a length-2 array.
        # There are 3x3 values to specify in combining 2 status values across
        # channels (wavelengths).  We combine according to these rules,
        # where s, s' are statuses {-1 (partial), 0 (fail), +1 (success)} --
        #   c(s,s) = s; c(s,s') = c(s',s); c(1, s) = 1; c(0, s) = s.
        # This turns out to be the same as the following:
        c = lambda s1, s2: np.sign(s1 + s2 + np.maximum(s1, s2))
        rv = 0.0 # start with the identity element
        for islice in obs['char_info']:
            rv = c(rv, islice['char_status'])
        return rv

get_koangles(OS)

Get koangles array from information in the OpticalSystem. The (un-pythonic) code here is copied from Prototypes/SurveySimulation so that we generate exactly the koangles as expected by the caching mechanism and the Observatory.keepout() method.

Source code in util/keepout_path_graphics.py

def get_koangles(OS):
    r'''Get koangles array from information in the OpticalSystem.
    The (un-pythonic) code here is copied from Prototypes/SurveySimulation so that 
    we generate exactly the koangles as expected by the caching mechanism and
    the Observatory.keepout() method.'''
    # choose observing modes selected for detection (default marked with a flag)
    allModes = OS.observingModes

    nSystems  = len(allModes)
    systNames = np.unique([allModes[x]['syst']['name'] for x in np.arange(nSystems)]).tolist()
    koStr     = ["koAngles_Sun", "koAngles_Moon", "koAngles_Earth", "koAngles_Small"]
    koangles  = np.zeros([len(systNames),4,2])
    tmpNames  = list(systNames)
    cnt = 0
    for x in np.arange(nSystems):
        name = allModes[x]['syst']['name']
        if name in tmpNames:
            koangles[cnt] = np.asarray([allModes[x]['syst'][k] for k in koStr])
            cnt += 1
            tmpNames.remove(name)
    return koangles

get_slew_begin(obs)

Return the time of the beginning of the slew, if any, for obs

Source code in util/keepout_path_graphics.py

def get_slew_begin(obs):
    r'''Return the time of the beginning of the slew, if any, for obs'''
    slew_time = strip_units(obs.get('slew_time', 0.0))
    arrival_time = strip_units(obs['arrival_time'])
    return arrival_time - slew_time

run_max(x)

Return the maximal run of True along each row of a 2d matrix x.

It is assumed that all rows start and end False, with one or more transitions to True along the way. You can False-pad x on left and right to ensure this.

Source code in util/keepout_path_graphics.py

def run_max(x):
    r"""Return the maximal run of True along each row of a 2d matrix x.

    It is assumed that all rows start and end False, with one or more
    transitions to True along the way.  You can False-pad x on left
    and right to ensure this."""
    # 0->1 and 1->0 transitions, and their row,col locations in x
    row,up   = np.where(np.diff(x * 1) > 0)
    row,down = np.where(np.diff(x * 1) < 0)
    # run-lengths for each run (0, 1, or >1 runs per row of x)
    run_lengths = down - up
    # tabulate maximal run per row of x
    x_runs = np.zeros(x.shape[0])
    for r in np.unique(row):
        x_runs[r] = np.max(run_lengths[row == r])
    return x_runs

spherical_cap(lon1, lat1, theta)

Return lon-lat point-list of a spherical cap of theta degrees around (lon,lat).

Envision the (lon,lat) point with a line connecting it to the origin of a unit sphere, and then a plane orthogonal to that line. To generate the cap, we construct a basis spanning the plane, and then draw a circle of radius tan(theta) in that plane. Then we extract the lon/lat values of that circle, and turn them into a list of connected segments. This is needed because the lon/lat circle can overflow the edge of the lon/lat plot, and appear as two half-circles.

Source code in util/keepout_path_graphics.py

def spherical_cap(lon1, lat1, theta):
    r'''Return lon-lat point-list of a spherical cap of theta degrees around (lon,lat).

    Envision the (lon,lat) point with a line connecting it to the origin of a unit sphere,
    and then a plane orthogonal to that line.  To generate the cap, we construct a basis 
    spanning the plane, and then draw a circle of radius tan(theta) in that plane.  Then
    we extract the lon/lat values of that circle, and turn them into a list of connected
    segments.  This is needed because the lon/lat circle can overflow the edge of the
    lon/lat plot, and appear as two half-circles.'''
    p_xyz = np.array(lonlat2xyz(lon1, lat1)).reshape((3,1))
    # orthogonal rotation matrix with p as one basis
    # the first column of R is p_xyz, and the other two (the ones we want) span the
    # orthongonal complement of p_xyz.  We use the latter two as a basis for the circle
    R = np.linalg.svd(p_xyz)[0]
    # radius of circle
    delta = np.tan(np.radians(theta))
    # sample [0,2pi], larger circles => more samples, very high-latitude => even more 
    twopi = np.linspace(0.0, 2*np.pi, int(theta*5+12) + (200 if theta > 50 else 0))
    # circle in the y-z plane of radius delta, with x == 0
    circle = np.stack((0*twopi, delta*np.cos(twopi), delta*np.sin(twopi)))
    # q = p + R * circle
    # rotate the circle to live in the plane perpendicular to p_xyz, and
    # translate the rotated circle to be centered about p_xyz
    q_xyz = p_xyz + np.dot(R, circle)
    # find the lon/lat corresponding to q
    lo, la = xyz2lonlat(q_xyz[0,:], q_xyz[1,:], q_xyz[2,:])
    # cut the circle if, e.g., lon jumps from 359.9 to 0.1,
    # making the circle into a list of line segments
    lola, closed = wrap_path(lo, la)
    return lola, closed

strip_units(x)

Strip astropy units from x.

Source code in util/keepout_path_graphics.py

def strip_units(x):
    r'''Strip astropy units from x.'''
    # TODO: allow coercing units to a supplied value
    if hasattr(x, 'value'):
        return x.value
    else:
        return x

wrap_path(x, y)

Given a path as lists of x's and y's, separate it into contiguous path segments.

The input paths are lists of (lon,lat), and if the jump between adjacent values is too large, this indicates a jump such as a wrap-around from 359.9 -> 0.1 in longitude. Think of a circle in (lon,lat) near the right edge of the plot. We separate the path at these points. The return value looks like: lines = [[(x1, y1), (x2, y2), ...], [(xP, yP), (xP+1, yP+1), ...] , ... ] We also return a boolean, "closed", indicating if the resulting line segments each enclose a closed boundary.

Source code in util/keepout_path_graphics.py

def wrap_path(x, y):
    r'''Given a path as lists of x's and y's, separate it into contiguous path segments.

    The input paths are lists of (lon,lat), and if the jump between adjacent values is
    too large, this indicates a jump such as a wrap-around from 359.9 -> 0.1 in longitude.
    Think of a circle in (lon,lat) near the right edge of the plot.
    We separate the path at these points.  The return value looks like:
       lines = [[(x1, y1), (x2, y2), ...], [(xP, yP), (xP+1, yP+1), ...] , ... ]
    We also return a boolean, "closed", indicating if the resulting line segments
    each enclose a closed boundary.'''
    lines = []
    line1 = []
    for i in range(len(x)):
        line1.append((x[i], y[i]))
        # index of next point, wrapping around
        i_next = i+1 if i+1 < len(x) else 0
        # start a new line segment if the next point will cause a big jump
        # note, we can jump in longitude, but not in latitude
        if abs(x[i] - x[i_next]) > 90:
            # finish the line appropriately depending on how the wrapping goes
            midpoint = (y[i] + y[i_next]) * 0.5
            if x[i] > 340 and x[i_next] < 20:
                # wrapping upward from 360 -> 0: cut line at 360, start at 0
                line1.append((360.0, midpoint))
                lines.append(line1)
                line1 = [(0.0, midpoint)]
            elif x[i] < 20 and x[i_next] > 340:
                # wrapping down from 0 to 360: cut at 0, start at 360
                line1.append((0.0, midpoint))
                lines.append(line1)
                line1 = [(360.0, midpoint)]
            else:
                # harder cases: near poles, longitude can jump but not wrap,
                # so in this case make no assumption: just delete the jump
                lines.append(line1)
                line1 = []
    # leftovers
    if len(line1) > 0:
        lines.append(line1)
    # perform some ad hoc shape joins
    if len(lines) == 3:
        # this special case joins circles that have been split
        # beginning <jump> middle <jump> end
        #   ==>
        # end + beginning <jump> middle
        lines = [lines[2] + lines[0], lines[1]]
        # typically this is a circle split across lon=360, so it is two closed shapes now
        closed = True
    elif len(lines) == 2:
        lines = [lines[1] + lines[0]]
        # typically this is a circle enclosing the pole, and it is not one closed shape
        # (note, it would be possible to close this by including the polar strip)
        closed = False
    else:
        # no jumps => boundary is closed
        closed = True
    return lines, closed

wrap_paths(p1, p2)

Construct a path from p1 to p2, allowing for wrap-around of longitude.

Note that latitude does not wrap, because +90 is distinct from -90.

Source code in util/keepout_path_graphics.py

def wrap_paths(p1, p2):
    r'''Construct a path from p1 to p2, allowing for wrap-around of longitude.

    Note that latitude does not wrap, because +90 is distinct from -90.'''
    # destructure the pairs
    x1, y1 = p1
    x2, y2 = p2

    # Construct a great-circle path from p1 to p2: a series of segments
    #   for reference, the one-segment path from p1 -> p2 is:
    #     lines = [[(x1, y1), (x2, y2)]]
    pts = gcp(x1, y1, x2, y2)
    lines = []
    # allow for longitude wrapping
    for i in range(len(pts)-1):
        p1 = pts[i]
        p2 = pts[i+1]
        # bisect line segments that wrap around
        if abs(p1[0] - p2[0]) > 90:
            # a midway latitude, not distance-weighted: ok for short segment
            midlat = (p1[1] + p2[1]) * 0.5
            if p1[0] > p2[0]:
                # 360 -> 0 wrap
                lines.append([(p1[0], p1[1] ), (360.0, midlat)])
                lines.append([(  0.0, midlat), (p2[0], p2[1] )])
            else:
                # 0 -> 360 wrap
                lines.append([(p1[0],  p1[1]), (  0.0, midlat)])
                lines.append([(360.0, midlat), (p2[0], p2[1] )])
        else:
            lines.append([(p1[0], p1[1]), (p2[0], p2[1])])
    return lines