Permalink
Browse files

merged Mike's changes in with mine

  • Loading branch information...
2 parents 560b595 + 8e0f110 commit 873f1ec53f55cd4ea875e9e30e7fda2cb63fd613 @barnabytprowe barnabytprowe committed Apr 10, 2012
Showing with 128 additions and 124 deletions.
  1. +121 −118 examples/Demo.py
  2. +7 −6 galsim/base.py
View
@@ -3,56 +3,6 @@
Some example scripts to see some basic usage of the GalSim library.
"""
-# List is issues that came up while making this. These should be spawned into
-# issues on GitHub.
-#
-# - Need a way to easily specify the size of the output image.
-# This should work to write to a particular portion of a larger image as
-# well, since that will be the normal thing we'll be doing.
-#
-# - Probably want to use the name Pixel for what is now Boxcar.
-# This means that SBPixel should probably also be renamed.
-#
-# - Probably don't want noise.addGaussian to have mean as a variable, since
-# it doesn't really make sense. Also, sigma should be a required variable,
-# not one with a default value of 1.
-#
-# - Should we move addGaussian to main galsim namespace?
-# If so, probably change the name to addGaussianNoise.
-# Maybe even make it a method of image, rather than a free function:
-# image.addGaussianNoise(rng, sigma = 0.3)
-# This seems clearer than galsim.noise.addGaussian(image,rng,simga)
-# Likewise:
-# image.addPoissonNoise(rng, gain = 1.7)
-#
-# - For Moffat, Sersic (others?) constructors:
-# re is not a very clear name for the half-light radius.
-# Should switch this to something more verbose, like half_light_radius.
-#
-# - Could we overload the + operator for GSObject?
-# So we could define a double gaussian as (atmos1 + atmos2)?
-# Then each gaussian could be separately sheared for example.
-# Or you might have an Airy + Gaussian. Or a triple Gaussian.
-# Better than defining specific new classes for every combination someone
-# might want.
-#
-# - If not the above, then we need to make GSAdd more user-friendly.
-# Currenly, its only constructor takes an SBAdd object, and we don't
-# really want the users to have to deal with that. (Or anything that
-# starts with SB I think.)
-#
-# - GSConvolve should drop the GS. GSAdd too.
-#
-# - Should get python versions of hsm, so we can do the hsm checks
-# directly here in the python code.
-#
-# - The applyShear function has the "wrong" conventions for the shear.
-# The arguments are taken to be (e1,e2) which are a distortion, not
-# a shear. Need to fix that. We can also have applyDistortion if we
-# want to have something that uses this convention, but it shouldn't
-# be what applyShear does.
-
-
import sys
import os
import subprocess
@@ -95,10 +45,14 @@ def __init__(self, file_name):
# These are distortions e1,e2
# Find the corresponding shear:
esq = self.e1*self.e1 + self.e2*self.e2
- e = math.sqrt(esq)
- g = math.tanh(0.5 * math.atanh(e))
- self.g1 = self.e1 * (g/e)
- self.g2 = self.e2 * (g/e)
+ if esq > 0.:
+ e = math.sqrt(esq)
+ g = math.tanh(0.5 * math.atanh(e))
+ self.g1 = self.e1 * (g/e)
+ self.g2 = self.e2 * (g/e)
+ else:
+ self.g1 = 0.
+ self.g2 = 0.
class HSM_Regauss:
"""
@@ -144,19 +98,25 @@ def Script1():
# In non-script code, use getLogger(__name__) at module scope instead.
logger = logging.getLogger("Script1")
- logger.info('Script 1:')
- logger.info('Starting script to convolve circular Gaussian galaxy (flux=1000, sigma=2),'\
- ' circular Gaussian PSF (flux=1, sigma=1),'\
- ' and pixel response (dx=0.2), then add Gaussian noise.')
+ gal_flux = 1.e5 # ADU
+ gal_sigma = 2. # pixels
+ psf_sigma = 1. # pixels
+ pixel_scale = 0.2 # arcsec / pixel
+ noise = 0.03 # ADU / pixel
+
+ logger.info('Starting script 1 using:')
+ logger.info(' - circular Gaussian galaxy (flux = %.1e, sigma = %.1f),',gal_flux,gal_sigma)
+ logger.info(' - circular Gaussian PSF (sigma = %.1f),',psf_sigma)
+ logger.info(' - pixel scale = %.2f,',pixel_scale)
+ logger.info(' - Gaussian noise (sigma = %.2f).',noise)
# Define the galaxy profile
- gal = galsim.Gaussian(flux=1000, sigma=2.)
+ gal = galsim.Gaussian(flux=gal_flux, sigma=gal_sigma)
# Define the PSF profile
- psf = galsim.Gaussian(flux=1., sigma=1.) # psf flux should always = 1
+ psf = galsim.Gaussian(flux=1., sigma=psf_sigma) # psf flux should always = 1
# Define the pixel size
- pixel_scale = 0.2 # arcsec / pixel
# Boxcar function to represent this pixellation
# The pixels could be rectangles, but normally xw = yw = pixel_scale
pix = galsim.Boxcar(xw=pixel_scale, yw=pixel_scale)
@@ -174,7 +134,7 @@ def Script1():
# Defaut seed is set from the current time.
rng = galsim.UniformDeviate()
# Use this to add Gaussian noise with specified sigma
- galsim.noise.addGaussian(image, rng, sigma=0.01)
+ galsim.noise.addGaussian(image, rng, sigma=noise)
# Write the image to a file
if not os.path.isdir('output'):
@@ -185,12 +145,12 @@ def Script1():
moments = HSM_Moments(file_name)
- logger.info('HSM reports that the image has measured moments Mxx, Myy, Mxy:'\
- ' %f, %f, %f', moments.mxx, moments.myy, moments.mxy)
- logger.info('e1,e2 = %f,%f', moments.e1, moments.e2)
- logger.info('g1,g2 = %f,%f', moments.g1, moments.g2)
- logger.info('Expected value for moments in limit that pixel response and noise are negligible:'\
- ' %f, %f, 0', (1.0**2+2.0**2)/(pixel_scale**2), (1.0**2+2.0**2)/(pixel_scale**2))
+ logger.info('HSM reports that the image has measured moments:')
+ logger.info(' Mxx = %.3f, Myy = %.3f, Mxy = %.3f', moments.mxx, moments.myy, moments.mxy)
+ logger.info('Expected values in the limit that pixel response and noise are negligible:')
+ mxx_exp = (gal_sigma**2+psf_sigma**2)/(pixel_scale**2) # == expected myy
+ logger.info(' Mxx = %.3f, Myy = %.3f, Mxy = %.3f', mxx_exp, mxx_exp,0)
+ print
# Script 2: Sheared, exponential galaxy, Moffat PSF, Poisson noise
def Script2():
@@ -202,23 +162,34 @@ def Script2():
"""
# In non-script code, use getLogger(__name__) at module scope instead.
logger = logging.getLogger("Script2")
- logger.info('Starting script to convolve sheared (0.1, 0.2) exponential galaxy,'\
- ' circular Moffat PSF,'\
- ' and pixel response, then add Poisson noise.')
+
+ gal_flux = 1.e5 # ADU
+ gal_r0 = 2.7 # pixels
+ g1 = 0.1 #
+ g2 = 0.2 #
+ psf_beta = 5 #
+ psf_re = 1.0 # pixels
+ pixel_scale = 0.2 # arcsec / pixel
+ sky_level = 1.e3 # ADU / pixel
+ gain = 1.0 # ADU / e-
+
+ logger.info('Starting script 2 using:')
+ logger.info(' - sheared (%.2f,%.2f) exponential galaxy (flux = %.1e, r0 = %.2f),',
+ g1, g2, gal_flux, gal_r0)
+ logger.info(' - circular Moffat PSF (beta = %.1f, re = %.2f),', psf_beta,psf_re)
+ logger.info(' - pixel scale = %.2f,', pixel_scale)
+ logger.info(' - Poisson noise (sky level = %.1e, gain = %.1f).', sky_level, gain)
# Define the galaxy profile.
- gal = galsim.Exponential(flux=1.e5, r0=2.7)
+ gal = galsim.Exponential(flux=gal_flux, r0=gal_r0)
# Shear the galaxy by some value.
- g1 = 0.1
- g2 = 0.2
gal.applyShear(g1, g2)
# Define the PSF profile.
- psf = galsim.Moffat(beta=5, flux=1., re=1.0)
+ psf = galsim.Moffat(beta=psf_beta, flux=1., re=psf_re)
# Define the pixel size
- pixel_scale = 0.2 # arcsec / pixel
pix = galsim.Boxcar(xw=pixel_scale, yw=pixel_scale)
# Final profile is the convolution of these.
@@ -230,7 +201,6 @@ def Script2():
image_epsf = final_epsf.draw(dx=pixel_scale)
# Add a constant sky level to the image.
- sky_level = 1.e4
# Create an image with the same bounds as image, with a constant
# sky level.
sky_image = galsim.ImageF(bounds=image.getBounds(), initValue=sky_level)
@@ -239,7 +209,7 @@ def Script2():
# This time use a particular seed, so it the image is deterministic.
rng = galsim.UniformDeviate(1534225)
# Use this to add Poisson noise.
- galsim.noise.addPoisson(image, rng, gain=1.)
+ galsim.noise.addPoisson(image, rng, gain=gain)
# Subtract off the sky.
image -= sky_image
@@ -257,12 +227,12 @@ def Script2():
moments = HSM_Moments(file_name)
moments_corr = HSM_Regauss(file_name, file_name_epsf, image.array.shape)
logger.info('HSM reports that the image has measured moments:')
- logger.info(' %f, %f, %f', moments.mxx, moments.myy, moments.mxy)
- logger.info('e1,e2 = %f,%f', moments.e1, moments.e2)
- logger.info('g1,g2 = %f,%f', moments.g1, moments.g2)
+ logger.info(' Mxx = %.3f, Myy = %.3f, Mxy = %.3f', moments.mxx, moments.myy, moments.mxy)
logger.info('When carrying out Regaussianization PSF correction, HSM reports')
- logger.info('e1,e2 = %f,%f', moments_corr.e1, moments_corr.e2)
- logger.info('g1,g2 = %f,%f', moments_corr.g1, moments_corr.g2)
+ logger.info(' g1,g2 = %.3f,%.3f', moments_corr.g1, moments_corr.g2)
+ logger.info('Expected values in the limit that noise and non-Gaussianity are negligible:')
+ logger.info(' g1,g2 = %.3f,%.3f', g1,g2)
+ print
# Script 3: Sheared, Sersic galaxy, Gaussian + OpticalPSF (atmosphere + optics) PSF, Poisson noise
@@ -280,75 +250,101 @@ def Script3():
"""
# In non-script code, use getLogger(__name__) at module scope instead.
logger = logging.getLogger("Script3")
- logger.info('Starting script to convolve sheared (0.1, -0.2) Sersic n=3 galaxy,'\
- ' optical PSF with defocus, coma, astigmatism'\
- ' non-circular double Gaussian atmospheric PSF'\
- ' pixel response + distortion,'\
- ' then add Poisson pixel noise and Gaussian read noise.')
-
+ gal_flux = 1.e5 # ADU
+ gal_n = 3.5 #
+ gal_re = 2.7 # pixels
+ g1 = -0.23 #
+ g2 = 0.15 #
+ atmos_sigma1=2.1 # pixels
+ atmos_sigma2=0.9 # pixels
+ atmos_f1=0.3 #
+ atmos_g1 = -0.13 #
+ atmos_g2 = -0.09 #
+ opt_defocus=5.0 # wavelengths
+ opt_a1=-2.9 # wavelengths
+ opt_a2=1.2 # wavelengths
+ opt_c1=6.4 # wavelengths
+ opt_c2=-3.3 # wavelengths
+ opt_padFactor=6 # multiples of Airy padding required to avoid folding for aberrated PSFs
+ lam = 800 # nm NB: don't use lambda - that's a reserved word.
+ tel_diam = 4. # meters
+ pixel_scale = 0.23 # arcsec / pixel
+ wcs_g1 = -0.02 #
+ wcs_g2 = 0.01 #
+ sky_level = 1.e3 # ADU / pixel
+ gain = 1.7 # ADU / e-
+ read_noise = 0.3 # ADU / pixel
+
+ logger.info('Starting script 3 using:')
+ logger.info(' - sheared (%.2f,%.2f) Sersic galaxy (flux = %.1e, n = %.1f, re = %.2f),',
+ g1, g2, gal_flux, gal_n, gal_re)
+ logger.info(' - sheared (%.2f,%.2f) double-Gaussian atmospheric PSF', atmos_g1,atmos_g2)
+ logger.info(' (sigma1 = %.2f, sigma2 = %.2f, frac1 = %.2f)',
+ atmos_sigma1, atmos_sigma2, atmos_f1)
+ logger.info(' - optical PSF with defocus = %.1f, astigmatism = (%.1f,%.1f),',
+ opt_defocus, opt_a1, opt_a2)
+ logger.info(' coma = (%.1f,%.1f), lambda = %.0f nm, D = %.1f m',
+ opt_c1, opt_c2, lam, tel_diam)
+ logger.info(' - pixel scale = %.2f,',pixel_scale)
+ logger.info(' - WCS distortion = (%.2f,%.2f),',wcs_g1,wcs_g2)
+ logger.info(' - Poisson noise (sky level = %.1e, gain = %.1f).',sky_level, gain)
+ logger.info(' - Gaussian read noise (sigma = %.2f).',read_noise)
+
+
# Define the galaxy profile.
- gal = galsim.Sersic(3.5, flux=1.e5, re=4.)
+ gal = galsim.Sersic(gal_n, flux=gal_flux, re=gal_re)
# Shear the galaxy by some value.
- g1 = 0.1
- g2 = -0.2
gal.applyShear(g1, g2)
# Define the atmospheric part of the PSF.
atmos = galsim.atmosphere.DoubleGaussian(
- flux1=0.3, sigma1=2.1, flux2=0.7, sigma2=0.9)
- atmos_g1 = -0.13
- atmos_g2 = -0.09
+ flux1=atmos_f1, sigma1=atmos_sigma1, flux2=1-atmos_f1, sigma2=atmos_sigma2)
atmos.applyShear(atmos_g1 , atmos_g2)
- # Define the pixel scale here, since we need it for the optical PSF
- pixel_scale = 0.23 # arcsec / pixel
-
# Define the optical part of the PSF.
# The first argument of OpticalPSF below is lambda/D,
# which needs to be in pixel units, so do the calculation:
- lam = 800 * 1.e-9 # 800 nm. NB: don't use lambda - that's a reserved word.
- D = 4. # 4 meter telescope, say.
- lam_over_D = lam / D # radians
+ lam_over_D = lam * 1.e-9 / tel_diam # radians
lam_over_D *= 206265 # arcsec
lam_over_D *= pixel_scale # pixels
- logger.info('lambda over D = %f', lam_over_D)
+ logger.info('Calculated lambda over D = %f pixels', lam_over_D)
# The rest of the values here should be given in units of the
# wavelength of the incident light. padFactor is used to here to reduce 'folding' for these
# quite strong aberration values
optics = galsim.OpticalPSF(lam_over_D,
- defocus=5., coma1=6.4, coma2=-3.3, astig1=-2.9, astig2=1.2,
- padFactor=6)
+ defocus=opt_defocus, coma1=opt_c1, coma2=opt_c2, astig1=opt_a1,
+ astig2=opt_a2, padFactor=opt_padFactor)
# Start with square pixels
pix = galsim.Boxcar(xw=pixel_scale, yw=pixel_scale)
- # Then shear them slightly
- wcs_g1 = -0.02
- wcs_g2 = 0.01
- pix.applyShear(wcs_g1, wcs_g2)
+ # Then shear them slightly by the negative of the wcs shear.
+ # This way the later distortion of the full image will bring them back to square.
+ pix.applyShear(-wcs_g1, -wcs_g2)
# Final profile is the convolution of these.
final = galsim.GSConvolve([gal, atmos, optics, pix])
final_epsf = galsim.GSConvolve([atmos, optics, pix])
+ # Now apply the wcs shear to the final image.
+ final.applyShear(wcs_g1, wcs_g2)
+ final_epsf.applyShear(wcs_g1, wcs_g2)
+
# Draw the image with a particular pixel scale.
image = final.draw(dx=pixel_scale)
image_epsf = final_epsf.draw(dx=pixel_scale)
# Draw the optical PSF component at its Nyquist sample rate
image_opticalpsf = optics.draw(dx=lam_over_D/2.)
# Add a constant sky level to the image.
- sky_level = 1.e4
sky_image = galsim.ImageF(bounds=image.getBounds(), initValue=sky_level)
image += sky_image
# Add Poisson noise to the image.
- gain = 1.7
rng = galsim.UniformDeviate(1314662)
galsim.noise.addPoisson(image, rng, gain=gain)
# Also add (Gaussian) read noise.
- read_noise = 0.3
galsim.noise.addGaussian(image, rng, sigma=read_noise)
# Subtract off the sky.
@@ -370,13 +366,14 @@ def Script3():
moments = HSM_Moments(file_name)
moments_corr = HSM_Regauss(file_name, file_name_epsf, image.array.shape)
- logger.info('HSM reports that the image has measured moments:'\
- ' %f, %f, %f', moments.mxx, moments.myy, moments.mxy)
- logger.info('e1,e2 = %f,%f', moments.e1, moments.e2)
- logger.info('g1,g2 = %f,%f', moments.g1, moments.g2)
+
+ logger.info('HSM reports that the image has measured moments:')
+ logger.info(' Mxx = %.3f, Myy = %.3f, Mxy = %.3f', moments.mxx, moments.myy, moments.mxy)
logger.info('When carrying out Regaussianization PSF correction, HSM reports')
- logger.info('e1,e2 = %f,%f', moments_corr.e1, moments_corr.e2)
- logger.info('g1,g2 = %f,%f', moments_corr.g1, moments_corr.g2)
+ logger.info(' g1,g2 = %f,%f', moments_corr.g1, moments_corr.g2)
+ logger.info('Expected values in the limit that noise and non-Gaussianity are negligible:')
+ logger.info(' g1,g2 = %f,%f', g1+wcs_g1,g2+wcs_g2)
+ print
def main(argv):
try:
@@ -392,10 +389,16 @@ def main(argv):
# For fancy logging setups (e.g. when running on a big cluster) we could
# use logging.fileConfig to use a config file to control logging.
logging.basicConfig(
- format="%(name)s[%(levelname)s]: %(message)s", # could also add date/time, pid, etc...
+ format="%(message)s",
level=logging.DEBUG,
stream=sys.stdout
)
+ # We do some fancier logging for Script3, just to demonstrate that we can:
+ # - we log to both stdout and to a log file
+ # - the log file has a lot more (mostly redundant) information
+ logFile = logging.FileHandler(os.path.join("output", "script3.log"))
+ logFile.setFormatter(logging.Formatter("%(name)s[%(levelname)s] %(asctime)s: %(message)s"))
+ logging.getLogger("Script3").addHandler(logFile)
# Script 1: Gaussian galaxy, Gaussian PSF, Gaussian noise.
if scriptNum == 0 or scriptNum == 1:
Oops, something went wrong.

0 comments on commit 873f1ec

Please sign in to comment.