Add astropy.units support to config layer.

docs/conf.py

@@ -68,6 +68,7 @@
     'sphinx.ext.autosectionlabel',
     'sphinx.ext.napoleon',
     'sphinx.ext.coverage',
+    'sphinx.ext.intersphinx',
     'breathe',
     'gh-link',
 ]
@@ -77,6 +78,11 @@
 breathe_default_project = "GalSim"
 breathe_default_members = ('members', 'undoc-members')
 
+# Link to astropy
+intersphinx_mapping = {
+    'astropy': ('http://docs.astropy.org/en/stable/', None),
+}
+
 # Add any paths that contain templates here, relative to this directory.
 #templates_path = ['_templates']
 

docs/config_image.rst

@@ -344,9 +344,9 @@ other types, including custom Bandpass types.
 * 'FileBandpass' is the default type here, and you may omit the type name when using it.
 
     * ``file_name`` = *str_value* (required)  The file to read in.
-    * ``wave_type`` = *str_value* (required)  The unit of the wavelengths in the file ('nm' or 'Ang' or variations on these -- cf. `Bandpass`)
-    * ``blue_limit`` = *float_value* (optional)  Hard cut off on the blue side.
-    * ``red_limit`` = *float value* (optional)  Hard cut off on the red side.
+    * ``wave_type`` = *str_value* or *unit_value* (required)  The unit of the wavelengths in the file ('nm' or 'Ang' or variations on these -- cf. `Bandpass`)
+    * ``blue_limit`` = *float_value* or *quantity_value* (optional)  Hard cut off on the blue side.
+    * ``red_limit`` = *float value* or *quantity_value* (optional)  Hard cut off on the red side.
     * ``zeropoint`` = *float_value* (optional)  The zero-point to use.
     * ``thin`` = *float_value* (optional)  If given, call `Bandpass.thin` on the Bandpass after reading in from the file, using this for the ``rel_err``.
 

docs/config_objects.rst

@@ -33,26 +33,26 @@ PSF Types
 * 'Airy'  A simple Airy disk. (Typically one would convolve this by some model of the atmospheric component of the PSF.  cf. 'Convolution' below.)
 
     * ``lam_over_diam`` = *float_value* (either ``lam_over_diam`` or both ``lam`` and ``diam`` required)  Lambda / telescope_diameter converted to units of arcsec (or whatever units you want your profile to use).
-    * ``lam`` = *float_value* (either ``lam_over_diam`` or both ``lam`` and ``diam`` required).  This should be the wavelength in nanometers.
-    * ``diam`` = *float_value* (either ``lam_over_diam`` or both ``lam`` and ``diam`` required).  This should be the telescope diameter in meters.
+    * ``lam`` = *float_value* or *quantity_value* (either ``lam_over_diam`` or both ``lam`` and ``diam`` required).  This should be the wavelength in nanometers.
+    * ``diam`` = *float_value* or *quantity_value* (either ``lam_over_diam`` or both ``lam`` and ``diam`` required).  This should be the telescope diameter in meters.
     * ``obscuration`` = *float_value* (default = 0)  The linear size of an obstructing secondary mirror as a fraction of the full mirror size.
     * ``scale_unit`` = *str_value* (default = 'arcsec') Units to be used for internal calculations when calculating lam/diam.
 
 * 'Kolmogorov'  A Kolmogorov turbulent spectrum: :math:`T(k) \sim \exp(-D(k)/2)`, where :math:`D(k) = 6.8839 (\lambda k/2\pi r0)^{5/3}`.
 
     * ``lam_over_r0`` = *float_value* (exactly one of ``lam_over_r0``, ``fwhm`` or ``half_light_radius`` or both ``lam`` and ``r0`` is required) Lambda / r0 converted to units of arcsec (or whatever units you want your profile to use).
-    * ``lam`` = *float_value* (exactly one of ``lam_over_r0``, ``fwhm`` or ``half_light_radius`` or both ``lam`` and ``r0`` is required) The wavelength in nanometers.
-    * ``r0`` = *float_value* (exactly one of ``lam_over_r0``, ``fwhm`` or ``half_light_radius`` or both ``lam`` and ``r0`` is required) The Fried parameter in meters.
-    * ``r0_500`` = *float_value* (optional, in lieu of ``r0``).  The Fried parameter in meters at a wavelength of 500 nm.  The correct ``r0`` value will be calculated using the standard relation r0 = r0_500 * (lam/500)``1.2.
+    * ``lam`` = *float_value* or *quantity_value* (exactly one of ``lam_over_r0``, ``fwhm`` or ``half_light_radius`` or both ``lam`` and ``r0`` is required) The wavelength in nanometers.
+    * ``r0`` = *float_value* or *quantity_value* (exactly one of ``lam_over_r0``, ``fwhm`` or ``half_light_radius`` or both ``lam`` and ``r0`` is required) The Fried parameter in meters.
+    * ``r0_500`` = *float_value* or *quantity_value* (optional, in lieu of ``r0``).  The Fried parameter in meters at a wavelength of 500 nm.  The correct ``r0`` value will be calculated using the standard relation r0 = r0_500 * (lam/500)``1.2.
     * ``fwhm`` = *float_value* (exactly one of ``lam_over_r0``, ``fwhm`` or ``half_light_radius`` or both ``lam`` and ``r0`` is required)
     * ``half_light_radius`` = *float_value* (exactly one of ``lam_over_r0``, ``fwhm`` or ``half_light_radius`` or both ``lam`` and ``r0`` is required)
     * ``scale_unit`` = *str_value* (default = 'arcsec') Units to be used for internal calculations when calculating lam/r0.
 
 * 'OpticalPSF'  A PSF from aberrated telescope optics.
 
     * ``lam_over_diam`` = *float_value* (either ``lam_over_diam`` or both ``lam`` and ``diam`` required)
-    * ``lam`` = *float_value* (either ``lam_over_diam`` or both ``lam`` and ``diam`` required).  This should be the wavelength in nanometers.
-    * ``diam`` = *float_value* (either ``lam_over_diam`` or both ``lam`` and ``diam`` required).  This should be the telescope diameter in meters.
+    * ``lam`` = *float_value* or *quantity_value* (either ``lam_over_diam`` or both ``lam`` and ``diam`` required).  This should be the wavelength in nanometers.
+    * ``diam`` = *float_value* or *quantity_value* (either ``lam_over_diam`` or both ``lam`` and ``diam`` required).  This should be the telescope diameter in meters.
     * ``defocus`` = *float_value* (default = 0) The defocus value, using the Noll convention for the normalization. (Noll index 4)
     * ``astig1`` = *float_value* (default = 0) The astigmatism in the y direction, using the Noll convention for the normalization. (Noll index 5)
     * ``astig2`` = *float_value* (default = 0) The astigmatism in the x direction, using the Noll convention for the normalization. (Noll index 6)
@@ -85,9 +85,9 @@ PSF Types
     * ``zenith_coord`` = *CelestialCoord*  (optional)  The (ra,dec) coordinate of the zenith.
     * ``HA`` = *Angle_value* (optional)  Hour angle of the observation.
     * ``latitude`` = *Angle_value* (optional)  Latitude of the observatory.
-    * ``pressure`` = *float_value* (default = 69.328)  Air pressure in kPa.
-    * ``temperature`` = *float_value* (default = 293.15)  Temperature in K.
-    * ``H2O_pressure`` = *float_value* (default = 1.067)  Water vapor pressure in kPa.
+    * ``pressure`` = *float_value* or *quantity_value* (default = 69.328)  Air pressure in kPa.
+    * ``temperature`` = *float_value* or *quantity_value* (default = 293.15)  Temperature in K.
+    * ``H2O_pressure`` = *float_value* or *quantity_value* (default = 1.067)  Water vapor pressure in kPa.
 
 Galaxy Types
 ------------
@@ -410,11 +410,11 @@ other types, including custom SED types.
 * 'FileSED' is the default type here, and you may omit the type name when using it.
 
     * ``file_name`` = *str_value* (required)  The file to read in.
-    * ``wave_type`` = *str_value* (required)  The unit of the wavelengths in the file ('nm' or 'Ang' or variations on these -- cf. `SED`)
+    * ``wave_type`` = *str_value* or *unit_value* (required)  The unit of the wavelengths in the file ('nm' or 'Ang' or variations on these -- cf. `SED`)
     * ``flux_type`` = *str_value* (required)  The type of spectral density or dimensionless normalization used in the file ('flambda', 'fnu', 'fphotons' or '1' -- cf. `SED`)
     * ``redshift`` = *float_value* (optional)  If given, shift the spectrum to the given redshift.  You can also specify the redshift as an object-level parameter if preferred.
-    * ``norm_flux_density`` = *float_value* (optional)  Set a normalization value of the flux density at a specific wavelength.  If given, ``norm_wavelength`` is required.
-    * ``norm_wavelength`` = *float_value* (optional)  The wavelength to use for the normalization flux density.
+    * ``norm_flux_density`` = *float_value* or *quantity_value* (optional)  Set a normalization value of the flux density at a specific wavelength.  If given, ``norm_wavelength`` is required.
+    * ``norm_wavelength`` = *float_value* or *quantity_value* (optional)  The wavelength to use for the normalization flux density.
     * ``norm_flux`` = *float_value* (optional)  Set a normalization value of the flux over a specific bandpass.  If given, ``norm_bandpass`` is required.
     * ``norm_bandpass`` = *Bandpass* (optional)  The bandpass to use for the normalization flux.
 

docs/config_stamp.rst

@@ -182,16 +182,16 @@ The photon operator types defined by GalSim are:
     * ``zenith_coord`` = *sky_value* (optional; see above) the celestial coordinates of the zenith.
     * ``HA`` = *angle_value* (optional; see above) the local hour angle.
     * ``latitude`` = *angle_value* (optional; see above) the latitude of the telescope.
-    * ``pressure`` = *float_value* (default = 69.328) the pressure in kPa.
-    * ``temperature`` = *float_value* (default = 293.15) the temperature in Kelvin.
-    * ``H2O_pressure`` = *float_value* (default = 1.067) the water vapor pressure in kPa.
+    * ``pressure`` = *float_value* or *quantity_value* (default = 69.328) the pressure in kPa.
+    * ``temperature`` = *float_value* or *quantity_value* (default = 293.15) the temperature in Kelvin.
+    * ``H2O_pressure`` = *float_value* or *quantity_value* (default = 1.067) the water vapor pressure in kPa.
 
 * 'FocusDepth' adjusts the positions of the photons at the surface of the sensor to account for
   the nominal focus being either above or below the sensor surface.  The depth value is typically
   negative, since the best focus is generally somewhere in the bulk of the sensor (although for
   short wavelengths it is often very close to the surface).
 
-    * ``depth`` = *float_value* (required)  The distance above the surface where the photons are
+    * ``depth`` = *float_value* (required)  The distance (in pixels) above the surface where the photons are
       nominally in focus.  A negative value means the focus in below the surface of the sensor.
 
 * 'Refraction' adjusts the incidence angles to account for refraction at the surface of the

docs/config_values.rst

@@ -13,6 +13,8 @@ have been designating using the following italic names:
 - *pos_value* corresponds to a GalSim `PositionD` instance.  See `pos_value`
 - *sky_value* corresponds to a GalSim `CelestialCoord` instance.  See `sky_value`
 - *table_value* corresponds to a GalSim `LookupTable` instance.  See `table_value`
+- *quantity_value* corresponds to an Astropy :class:`~astropy.units.Quantity` instance.  See `quantity_value`
+- *unit_value* corresponds to an Astropy :class:`~astropy.units.Unit` instance.  See `unit_value`
 
 Each of the Python types can be given as a constant value using the normal Python conventions
 for how to specify such a value.  The GalSim *angle_value* and *pos_value* also have
@@ -525,7 +527,7 @@ Options are:
 
 * A dict with:
 
-    * ``type`` = *str* (required)  There is currently only one valid option:
+    * ``type`` = *str* (required)  Valid options are:
 
         * 'RADec' Specify ra and dec separately.
 
@@ -571,6 +573,49 @@ Options are:
 
 * A string that starts with '$' or '@'.  See `Shorthand Notation`.
 
+quantity_value
+--------------
+
+These represent `astropy.units.Quantity` values, which are a combination of a float and a unit (specifically, an `astropy.units.Unit`).
+
+Options are:
+
+* An `astropy.units.Quantity` value (e.g. '8.7*u.m', where 'u' is the `astropy.units` module)
+* A string interpretable by `astropy.units.Quantity` (e.g. '8.7 m')
+* A dict with:
+
+    * ``type`` = *str* (required)  Valid options are:
+
+        * 'Quantity' Specify the value and unit separately.
+
+            * ``value`` = *float_value* (required)
+            * ``unit`` = *unit_value* (required)
+
+        * 'Eval'  Evaluate a string.  See `Eval type`.
+
+* A string that starts with '$' or '@'.  See `Shorthand Notation`.
+
+unit_value
+----------
+
+These represent `astropy.units.Unit` values.
+
+Options are:
+
+* An `astropy.units.Unit` value (e.g. 'u.m', where 'u' is the `astropy.units` module)
+* A string interpretable by `astropy.units.Unit` (e.g. 'm')
+* A dict with:
+
+    * ``type`` = *str* (required)  Valid options are:
+
+        * 'Unit' Specify the unit.
+
+            * ``unit`` = *str_value* (required)
+
+        * 'Eval'  Evaluate a string.  See `Eval type`.
+
+* A string that starts with '$' or '@'.  See `Shorthand Notation`.
+
 
 Eval type
 ---------

examples/demo11.py

@@ -81,6 +81,7 @@ def main(argv):
         observation.  However, we whiten the noise of the final image so the final image has
         stationary Gaussian noise, rather than correlated noise.
     """
+    import astropy.units as u
     logging.basicConfig(format="%(message)s", level=logging.INFO, stream=sys.stdout)
     logger = logging.getLogger("demo11")
 
@@ -95,7 +96,10 @@ def main(argv):
                                       # (This corresponds to 8 galaxies / arcmin^2)
     grid_spacing = 90.0               # The spacing between the samples for the power spectrum
                                       # realization (arcsec)
-    tel_diam = 4                      # Let's figure out the flux for a 4 m class telescope
+    tel_diam = 400*u.cm               # Let's figure out the flux for a 4 m class telescope.  We'll
+                                      # also show how astropy units can be used here (see especially
+                                      # the eval_variables section of the yaml version of this
+                                      # demo).
     exp_time = 300                    # exposing for 300 seconds.
     center_ra = 19.3*galsim.hours     # The RA, Dec of the center of the image on the sky
     center_dec = -33.1*galsim.degrees
@@ -106,7 +110,7 @@ def main(argv):
     # is 2.4 but there is obscuration (a linear factor of 0.33).  Here, we assume that the telescope
     # we're simulating effectively has no obscuration factor.  We're also ignoring the pi/4 factor
     # since it appears in the numerator and denominator, so we use area = diam^2.
-    hst_eff_area = 2.4**2 * (1.-0.33**2)
+    hst_eff_area = (2.4*u.m)**2 * (1.-0.33**2)
     flux_scaling = (tel_diam**2/hst_eff_area) * exp_time
 
     # random_seed is used for both the power spectrum realization and the random properties

examples/demo11.yaml

@@ -44,6 +44,7 @@
 # - image.noise : symmetrize
 # - wcs type : Tan(dudx, dudy, dvdx, dvdy, units, origin, ra, dec)
 # - top level field eval_variables
+# - quantity letter type for eval variable.
 #
 # - Power spectrum shears and magnifications for non-gridded positions.
 # - Reading a compressed FITS image (using BZip2 compression).
@@ -74,8 +75,10 @@ eval_variables :
     # An unusual prefix: a = Angle.
     atheta : &theta 0.17 degrees
 
-    # tel_diam = the telescope diameter in meters.
-    ftel_diam : &tel_diam 4
+    # tel_diam = the telescope diameter.  Prefixing this with a 'q' means it is an
+    # astropy.units.Quantity.  We'll specify centimeters here, note the interaction with
+    # meters down below when we use it to calculate the scale_flux.
+    qtel_diam : &tel_diam 400 cm
 
     # exp_time = the exposure time in seconds.
     fexp_time : &exp_time 300
@@ -167,7 +170,7 @@ gal :
     #     noise_pad_size = 40 * sqrt(2) * 0.2 arcsec/pixel = 11.3 arcsec
     noise_pad_size : 11.3
 
-    # We will select a galaxy at random from the catalog. One could easily do this by setting 
+    # We will select a galaxy at random from the catalog. One could easily do this by setting
     # index : { type : Random }
     # but we will instead allow the catalog to choose a random galaxy for us.  It will remove any
     # selection effects involved in postage stamp creation using weights that are stored in the
@@ -192,7 +195,7 @@ gal :
         # second exposures.  The COSMOS galaxies have their flux set to be appropriate for HST
         # (a 2.4 m telescope with a linear obscuration factor of 0.33) with 1 second exposures.
         # So the flux should be scaled by (4**2/(2.4*(1-0.33**2))) * 300
-        "$(tel_diam**2 / (2.4**2*(1.-0.33**2))) * exp_time"
+        "$(tel_diam**2 / ((2.4*u.m)**2*(1.-0.33**2))) * exp_time"
 
 
 # Define some other information about the images

examples/json/demo11.json

@@ -19,7 +19,7 @@
 "eval_variables" : {
     "fpixel_scale" : 0.2,
     "atheta" : "0.17 degrees",
-    "ftel_diam" : 4,
+    "qtel_diam" : "400 cm",
     "fexp_time" : 300,
     "fimage_size" : 2048,
     "inobjects" : 288,
@@ -44,7 +44,7 @@
     "shear" : { "type" : "PowerSpectrumShear" },
     "magnification" : { "type" : "PowerSpectrumMagnification" },
     "rotation" : { "type" : "Random" },
-    "scale_flux" : "$(tel_diam**2 / (2.4**2*(1.-0.33**2))) * exp_time"
+    "scale_flux" : "$(tel_diam**2 / ((2.4*u.m)**2*(1.-0.33**2))) * exp_time"
 },
 
 "image" : {

galsim/airy.py

@@ -19,6 +19,7 @@
 __all__ = [ 'Airy' ]
 
 import math
+import astropy.units as u
 
 from . import _galsim
 from .gsobject import GSObject
@@ -63,12 +64,12 @@ class Airy(GSObject):
     Parameters:
         lam_over_diam:  The parameter that governs the scale size of the profile.
                         See above for details about calculating it.
-        lam:            Lambda (wavelength) in units of nanometers.  Must be supplied with
-                        ``diam``, and in this case, image scales (``scale``) should be specified
-                        in units of ``scale_unit``.
-        diam:           Telescope diameter in units of meters.  Must be supplied with
-                        ``lam``, and in this case, image scales (``scale``) should be specified
-                        in units of ``scale_unit``.
+        lam:            Lambda (wavelength) either as an astropy Quantity, or as a float in units
+                        of nanometers.  Must be supplied with ``diam``, and in this case, image
+                        scales (``scale``) should be specified in units of ``scale_unit``.
+        diam:           Telescope diameter either as an astropy Quantity, or as a float in units of
+                        meters.  Must be supplied with ``lam``, and in this case, image scales
+                        (``scale``) should be specified in units of ``scale_unit``.
         obscuration:    The linear dimension of a central obscuration as a fraction of the
                         pupil dimension.  [default: 0]
         flux:           The flux (in photons/cm^2/s) of the profile. [default: 1]
@@ -80,13 +81,16 @@ class Airy(GSObject):
         gsparams:       An optional `GSParams` argument. [default: None]
     """
     _req_params = { }
-    _opt_params = { "flux" : float , "obscuration" : float, "diam" : float,
-                    "scale_unit" : str }
+    _opt_params = { "flux" : float ,
+                    "obscuration" : float,
+                    "diam" : (float, u.Quantity),
+                    "scale_unit" : str
+                  }
     # Note that this is not quite right; it's true that either lam_over_diam or lam should be
     # supplied, but if lam is supplied then diam is required.  Errors in which parameters are used
     # may be caught either by config or by the python code itself, depending on the particular
     # error.
-    _single_params = [{ "lam_over_diam" : float , "lam" : float } ]
+    _single_params = [{ "lam_over_diam" : float , "lam" : (float, u.Quantity) } ]
 
     # For an unobscured Airy, we have the following factor which can be derived using the
     # integral result given in the Wikipedia page (http://en.wikipedia.org/wiki/Airy_disk),
@@ -108,6 +112,11 @@ def __init__(self, lam_over_diam=None, lam=None, diam=None, obscuration=0., flux
         self._flux = float(flux)
         self._gsparams = GSParams.check(gsparams)
 
+        if isinstance(lam, u.Quantity):
+            lam = lam.to_value(u.nm)
+        if isinstance(diam, u.Quantity):
+            diam = diam.to_value(u.m)
+
         # Parse arguments: either lam_over_diam in arbitrary units, or lam in nm and diam in m.
         # If the latter, then get lam_over_diam in units of scale_unit, as specified in
         # docstring.