Merge pull request #89 from lsst/u/yoachim/add-solar-elong

U/yoachim/add solar elong

python/lsst/sims/featureScheduler/basis_functions/basis_functions.py

@@ -74,7 +74,7 @@ def check_feasibility(self, conditions):
         """
         return True
 
-    def _calc_value(self, conditions, **kwarge):
+    def _calc_value(self, conditions, **kwargs):
         self.value = 0
         # Update the last time we had an mjd
         self.mjd_last = conditions.mjd + 0
@@ -557,16 +557,18 @@ class Near_sun_twilight_basis_function(Base_basis_function):
         The maximum airmass to try and observe (unitless)
     """
 
-    def __init__(self, nside=None, max_airmass=2.5):
+    def __init__(self, nside=None, max_airmass=2.5, penalty=np.nan):
         super(Near_sun_twilight_basis_function, self).__init__(nside=nside)
         self.max_airmass = int_rounded(max_airmass)
-        self.result = np.zeros(hp.nside2npix(self.nside))
+        self.result = np.empty(hp.nside2npix(self.nside))
+        self.result.fill(penalty)
 
     def _calc_value(self, conditions, indx=None):
         result = self.result.copy()
-        good_pix = np.where((int_rounded(conditions.airmass) < self.max_airmass) &
-                            (int_rounded(conditions.az_to_sun) < int_rounded(np.pi/2.)))
-        result[good_pix] = conditions.airmass[good_pix] / self.max_airmass.value
+        good_pix = np.where((conditions.airmass >= 1.) & 
+                            (int_rounded(conditions.airmass) < self.max_airmass) &
+                            (int_rounded(np.abs(conditions.az_to_sun)) < int_rounded(np.pi/2.)))
+        result[good_pix] = conditions.airmass[good_pix] / self.max_airmass.initial
         return result
 
 

python/lsst/sims/featureScheduler/basis_functions/mask_basis_funcs.py

@@ -8,7 +8,68 @@
 
 __all__ = ['Zenith_mask_basis_function', 'Zenith_shadow_mask_basis_function',
            'Moon_avoidance_basis_function', 'Map_cloud_basis_function',
-           'Planet_mask_basis_function', 'Mask_azimuth_basis_function']
+           'Planet_mask_basis_function', 'Mask_azimuth_basis_function',
+           'Solar_elongation_mask_basis_function', 'Area_check_mask_basis_function']
+
+
+
+class Area_check_mask_basis_function(Base_basis_function):
+    """Take a list of other mask basis functions, and do an additional check for area available
+    """
+    def __init__(self, bf_list, nside=32, min_area=1000.):
+        super(Area_check_mask_basis_function, self).__init__(nside=nside)
+        self.bf_list = bf_list
+        self.result = np.zeros(hp.nside2npix(self.nside), dtype=float)
+        self.min_area = min_area
+
+    def check_feasibility(self, conditions):
+        result = True
+        for bf in self.bf_list:
+            if not bf.check_feasibility(conditions):
+                return False
+
+        area_map = self.result.copy()
+        for bf in self.bf_list:
+            area_map *= bf(conditions)
+
+        good_pix = np.where(area_map == 0)[0]
+        if hp.nside2pixarea(self.nside, degrees=True)*good_pix.size < self.min_area:
+            result = False
+        return result
+
+    def _calc_value(self, conditions, **kwargs):
+        result = self.result.copy()
+        for bf in self.bf_list:
+            result *= bf(conditions)
+        return result
+
+
+
+class Solar_elongation_mask_basis_function(Base_basis_function):
+    """Mask things at various solar elongations
+
+    Parameters
+    ----------
+    min_elong : float (0)
+        The minimum solar elongation to consider (degrees).
+    max_elong : float (60.)
+        The maximum solar elongation to consider (degrees).
+    """
+
+    def __init__(self, min_elong=0., max_elong=60., nside=None, penalty=np.nan):
+        super(Solar_elongation_mask_basis_function, self).__init__(nside=nside)
+        self.min_elong = np.radians(min_elong)
+        self.max_elong = np.radians(max_elong)
+        self.penalty = penalty
+        self.result = np.empty(hp.nside2npix(self.nside), dtype=float)
+        self.result.fill(self.penalty)
+
+    def _calc_value(self, conditions, indx=None):
+        result = self.result.copy()
+        in_range = np.where((int_rounded(conditions.solar_elongation) >= int_rounded(self.min_elong)) &
+                            (int_rounded(conditions.solar_elongation) <= int_rounded(self.max_elong)))[0]
+        result[in_range] = 1
+        return result
 
 
 class Zenith_mask_basis_function(Base_basis_function):
@@ -21,8 +82,8 @@ class Zenith_mask_basis_function(Base_basis_function):
     max_alt : float (82.)
         The maximum allowed altitude (degrees)
     """
-    def __init__(self, min_alt=20., max_alt=82.):
-        super(Zenith_mask_basis_function, self).__init__()
+    def __init__(self, min_alt=20., max_alt=82., nside=None):
+        super(Zenith_mask_basis_function, self).__init__(nside=nside)
         self.update_on_newobs = False
         self.min_alt = np.radians(min_alt)
         self.max_alt = np.radians(max_alt)

python/lsst/sims/featureScheduler/features/conditions.py

@@ -1,5 +1,5 @@
 import numpy as np
-from lsst.sims.utils import _approx_RaDec2AltAz, Site, _hpid2RaDec, m5_flat_sed, _approx_altaz2pa
+from lsst.sims.utils import _approx_RaDec2AltAz, Site, _hpid2RaDec, m5_flat_sed, _approx_altaz2pa, _angularSeparation
 import healpy as hp
 from lsst.sims.featureScheduler.utils import set_default_nside, match_hp_resolution, season_calc, smallest_signed_angle
 
@@ -151,6 +151,8 @@ def __init__(self, nside=None, site='LSST', exptime=30., mjd_start=59853.5, seas
             Healpix map of the azimuthal distance to the sun for each healpixel (radians)
         az_to_anitsun : np.array
             Healpix map of the azimuthal distance to the anit-sun for each healpixel (radians)
+        solar_elongation : np.array
+            Healpix map of the solar elongation (angular distance to the sun) for each healpixel (radians)
 
         Attributes (set by the scheduler)
         -------------------------------
@@ -338,6 +340,7 @@ def mjd(self, value):
         self._az_to_sun = None
         self._az_to_antisun = None
         self._season = None
+        self._solar_elongation = None
 
     @property
     def skybrightness(self):
@@ -378,6 +381,15 @@ def calc_M5Depth(self):
                                                           self.exptime,
                                                           self._airmass[good])
 
+    def calc_solar_elongation(self):
+        self._solar_elongation = _angularSeparation(self.ra, self.dec, self.sunRA, self.sunDec)
+
+    @property
+    def solar_elongation(self):
+        if self._solar_elongation is None:
+            self.calc_solar_elongation()
+        return self._solar_elongation
+
     def calc_az_to_sun(self):
         self._az_to_sun = smallest_signed_angle(self.ra, self.sunRA)
 

python/lsst/sims/featureScheduler/modelObservatory/model_observatory.py

@@ -443,7 +443,7 @@ def observe(self, observation):
         # obsevation['rotSkyPos'] will be ignored.
         slewtime, visittime = self.observatory.observe(observation, self.mjd, rotTelPos=observation['rotTelPos'])
 
-        # inf slewtime means the observation failed (probably outsire alt limits)
+        # inf slewtime means the observation failed (probably outside alt limits)
         if ~np.all(np.isfinite(slewtime)):
             return None, False
 

python/lsst/sims/featureScheduler/sim_runner.py

@@ -52,6 +52,8 @@ def sim_runner(observatory, scheduler, filter_scheduler=None, mjd_start=None, su
     nskip = 0
     new_night = False
 
+    mjd_last_flush = -1
+
     while mjd < end_mjd:
         if not scheduler._check_queue_mjd_only(observatory.mjd):
             scheduler.update_conditions(observatory.return_conditions())
@@ -69,7 +71,12 @@ def sim_runner(observatory, scheduler, filter_scheduler=None, mjd_start=None, su
             observations.append(completed_obs)
             filter_scheduler.add_observation(completed_obs[0])
         else:
+            # An observation failed to execute, usually it was outside the altitude limits.
+            if observatory.mjd == mjd_last_flush:
+                raise RuntimeError("Scheduler has failed to provide a valid observation multiple times.")
+            # if this is a first offence, might just be that targets set. Flush queue and get some new targets.
             scheduler.flush_queue()
+            mjd_last_flush = observatory.mjd + 0
         if new_night:
             # find out what filters we want mounted
             conditions = observatory.return_conditions()
@@ -88,7 +95,7 @@ def sim_runner(observatory, scheduler, filter_scheduler=None, mjd_start=None, su
             if len(observations) == n_visit_limit:
                 break
         # XXX--handy place to interupt and debug
-        # if len(observations) > 3:
+        #if len(observations) > 25:
         #    import pdb ; pdb.set_trace()
     runtime = time.time() - t0
     print('Skipped %i observations' % nskip)

python/lsst/sims/featureScheduler/surveys/base_survey.py

@@ -174,11 +174,13 @@ class BaseMarkovDF_survey(BaseSurvey):
         Random number seed, used for randomly orienting sky tessellation.
     camera : str ('LSST')
         Should be 'LSST' or 'comcam'
+    area_required : float (None)
+        The valid area that should be present in the reward function (square degrees).
     """
     def __init__(self, basis_functions, basis_weights, extra_features=None,
                  smoothing_kernel=None,
                  ignore_obs=None, survey_name='', nside=None, seed=42,
-                 dither=True, detailers=None, camera='LSST'):
+                 dither=True, detailers=None, camera='LSST', area_required=None):
 
         super(BaseMarkovDF_survey, self).__init__(basis_functions=basis_functions,
                                                   extra_features=extra_features,
@@ -197,7 +199,7 @@ def __init__(self, basis_functions, basis_weights, extra_features=None,
         elif self.camera == 'comcam':
             self.fields_init = comcamTessellate()
         else:
-            ValueError('camera %s unknown, should be "LSST" or "comcam"' %camera)
+            ValueError('camera %s unknown, should be "LSST" or "comcam"' % camera)
         self.fields = self.fields_init.copy()
         self.hp2fields = np.array([])
         self._hp2fieldsetup(self.fields['RA'], self.fields['dec'])
@@ -207,13 +209,34 @@ def __init__(self, basis_functions, basis_weights, extra_features=None,
         else:
             self.smoothing_kernel = None
 
+        if area_required is None:
+            self.area_required = area_required
+        else:
+            self.area_required = area_required * (np.pi/180.)**2  # To steradians
+
         # Start tracking the night
         self.night = -1
 
         # Set the seed
         np.random.seed(seed)
         self.dither = dither
 
+    def _check_feasibility(self, conditions):
+        """
+        Check if the survey is feasable in the current conditions
+        """
+        for bf in self.basis_functions:
+            result = bf.check_feasibility(conditions)
+            if not result:
+                return result
+        if self.area_required is not None:
+            reward = self.calc_reward_function(conditions)
+            good_pix = np.where(np.isfinite(reward) == True)[0]
+            area = hp.nside2pixarea(self.nside) * np.size(good_pix)
+            if area < self.area_required:
+                return False
+        return result
+
     def _hp2fieldsetup(self, ra, dec, leafsize=100):
         """Map each healpixel to nearest field. This will only work if healpix
         resolution is higher than field resolution.
@@ -307,6 +330,11 @@ def calc_reward_function(self, conditions):
         if self.smoothing_kernel is not None:
             self.smooth_reward()
 
+        if self.area_required is not None:
+            good_area = np.where(np.abs(self.reward) >= 0)[0].size * hp.nside2pixarea(self.nside)
+            if good_area < self.area_required:
+                self.reward = -np.inf
+
         return self.reward
 
     def generate_observations_rough(self, conditions):

python/lsst/sims/featureScheduler/surveys/surveys.py

@@ -19,7 +19,7 @@ class Greedy_survey(BaseMarkovDF_survey):
     def __init__(self, basis_functions, basis_weights, filtername='r',
                  block_size=1, smoothing_kernel=None, nside=None,
                  dither=True, seed=42, ignore_obs=None, survey_name='',
-                 nexp=2, exptime=30., detailers=None, camera='LSST'):
+                 nexp=2, exptime=30., detailers=None, camera='LSST', area_required=None):
 
         extra_features = {}
 
@@ -30,7 +30,8 @@ def __init__(self, basis_functions, basis_weights, filtername='r',
                                             ignore_obs=ignore_obs,
                                             nside=nside,
                                             survey_name=survey_name, dither=dither,
-                                            detailers=detailers, camera=camera)
+                                            detailers=detailers, camera=camera,
+                                            area_required=area_required)
         self.filtername = filtername
         self.block_size = block_size
         self.nexp = nexp
@@ -132,7 +133,7 @@ def __init__(self, basis_functions, basis_weights,
                  dither=True, seed=42, ignore_obs=None,
                  survey_note='blob', detailers=None, camera='LSST',
                  twilight_scale=True, in_twilight=False, check_scheduled=True, min_area=None,
-                 grow_blob=True):
+                 grow_blob=True, area_required=None):
 
         if nside is None:
             nside = set_default_nside()
@@ -142,7 +143,8 @@ def __init__(self, basis_functions, basis_weights,
                                           filtername=None,
                                           block_size=0, smoothing_kernel=smoothing_kernel,
                                           dither=dither, seed=seed, ignore_obs=ignore_obs,
-                                          nside=nside, detailers=detailers, camera=camera)
+                                          nside=nside, detailers=detailers, camera=camera,
+                                          area_required=area_required)
         self.flush_time = flush_time/60./24.  # convert to days
         self.nexp = nexp
         self.nexp_dict = nexp_dict
@@ -305,6 +307,15 @@ def calc_reward_function(self, conditions):
             self.reward[az_out] = np.nan
         else:
             self.reward = -np.inf
+
+        if self.area_required is not None:
+            good_area = np.where(np.abs(self.reward) >= 0)[0].size * hp.nside2pixarea(self.nside)
+            if good_area < self.area_required:
+                self.reward = -np.inf
+
+        #if ('twi' in self.survey_note) & (np.any(np.isfinite(self.reward))):
+        #    import pdb ; pdb.set_trace()
+
         self.reward_checked = True
         return self.reward
 

python/lsst/sims/featureScheduler/utils/footprints.py

@@ -22,7 +22,7 @@
            'WFD_healpixels', 'WFD_no_gp_healpixels', 'WFD_bigsky_healpixels', 'WFD_no_dust_healpixels',
            'SCP_healpixels', 'NES_healpixels',
            'galactic_plane_healpixels', #'low_lat_plane_healpixels', 'bulge_healpixels',
-           'magellanic_clouds_healpixels',
+           'magellanic_clouds_healpixels', 'Constant_footprint',
            'generate_goal_map', 'standard_goals',
            'calc_norm_factor', 'filter_count_ratios', 'Step_line', 'Footprints', 'Footprint',
            'Step_slopes']
@@ -176,6 +176,21 @@ def __call__(self, mjd, array=False):
         return self.arr2struc(self.current_footprints)
 
 
+class Constant_footprint(Footprint):
+    def __init__(self, nside=32,
+                 filters={'u': 0, 'g': 1, 'r': 2, 'i': 3, 'z': 4, 'y': 5}):
+        self.nside = nside
+        self.filters = filters
+        self.npix = hp.nside2npix(nside)
+        self.footprints = np.zeros((len(filters), self.npix), dtype=float)
+        self.out_dtype = list(zip(filters, [float]*len(filters)))
+        self.to_return = self.arr2struc(self.footprints)
+
+    def __call__(self, mjd, array=False):
+        return self.to_return
+
+
+
 class Footprints(Footprint):
     """An object to combine multiple Footprint objects.
     """