#103 - abbreviated config, final multi-object demo #161

Merged
merged 20 commits into from May 24, 2012
Commits
+318 −72
Split
View
@@ -34,7 +34,7 @@ def Script1():
- Noise is poisson using a nominal sky value of 1.e6.
- Galaxies are sersic profiles.
"""
- logger = logging.getLogger("Script1")
+ logger = logging.getLogger("Script1")
# Define some parameters we'll use below.
# Normally these would be read in from some parameter file.
@@ -52,15 +52,15 @@ def Script1():
psf_file_name = os.path.join('output','g08_psf.fits')
psf_beta = 3 #
psf_fwhm = 2.85 # arcsec (=pixels)
- psf_trunc = 2. # FWHM
+ psf_trunc = 2. # FWHM
psf_g1 = -0.019 #
psf_g2 = -0.007 #
psf_centroid_shift = 1.0 # arcsec (=pixels)
gal_file_name = os.path.join('output','g08_gal.fits')
gal_signal_to_noise = 200 # Great08 "LowNoise" run
gal_n = 1 # Use n=1 for "disk" type
- # Great08 mixed pure bulge and pure disk for its LowNoise run.
+ # Great08 mixed pure bulge and pure disk for its LowNoise run.
# We're just doing disks to make things simpler.
gal_resolution = 1.4 # r_obs / r_psf (use r = half_light_radius)
gal_centroid_shift = 1.0 # arcsec (=pixels)
@@ -72,7 +72,7 @@ def Script1():
logger.info('Starting multi-object script 1 using:')
logger.info(' - image with %d x %d postage stamps',nx_stamps,ny_stamps)
logger.info(' - postage stamps of size %d x %d pixels',nx_pixels,ny_pixels)
- logger.info(' - Moffat PSF (beta = %.1f, FWHM = %.2f, trunc = %.2f),',
+ logger.info(' - Moffat PSF (beta = %.1f, FWHM = %.2f, trunc = %.2f),',
psf_beta,psf_fwhm,psf_trunc)
logger.info(' - PSF shear = (%.3f,%.3f)',psf_g1,psf_g2)
logger.info(' - PSF centroid shifts up to = %.2f pixels',psf_centroid_shift)
@@ -106,7 +106,7 @@ def Script1():
logger.info('i,j,x,y')
for ix in range(nx_stamps):
for iy in range(ny_stamps):
- # The -2's in the next line rather than -1 are to provide a border of
+ # The -2's in the next line rather than -1 are to provide a border of
# 1 pixel between postage stamps
b = galsim.BoundsI(ix*nx_pixels , (ix+1)*nx_pixels -2,
iy*ny_pixels , (iy+1)*ny_pixels -2)
@@ -141,7 +141,7 @@ def Script1():
# Now write the image to disk.
psf_image.write(psf_file_name, clobber=True)
logger.info('Wrote PSF file %s',psf_file_name)
-
+
# Define the galaxy profile
# TODO: This is a hack. We should have a getHalfLightRadius() method.
@@ -150,7 +150,7 @@ def Script1():
# First figure out the size we need from the resolution
# great08 resolution was defined as Rgp / Rp where Rp is the FWHM of the PSF
# and Rgp is the FWHM of the convolved galaxy.
- # We make the approximation here that the FWHM adds in quadrature during the
+ # We make the approximation here that the FWHM adds in quadrature during the
# convolution, so we can get the unconvolved size as:
# Rg^2 = Rgp^2 - Rp^2 = Rp^2 * (resolution^2 - 1)
gal_re = psf_re * math.sqrt( gal_resolution**2 - 1)
@@ -163,7 +163,7 @@ def Script1():
# We use a weighted integral of the flux:
# S = sum W(x,y) I(x,y) / sum W(x,y)
# N^2 = Var(S) = sum W(x,y)^2 Var(I(x,y)) / (sum W(x,y))^2
- # Now we assume that Var(I(x,y)) is dominated by the sky noise, so
+ # Now we assume that Var(I(x,y)) is dominated by the sky noise, so
# Var(I(x,y)) = sky_level
# We also assume that we are using a matched filter for W, so W(x,y) = I(x,y).
# Then a few things cancel and we find that
@@ -174,7 +174,7 @@ def Script1():
gal *= gal_signal_to_noise / sn_meas
logger.info('Made galaxy profile')
- # This profile is placed with different orientations and noise realizations
+ # This profile is placed with different orientations and noise realizations
# at each postage stamp in the gal image.
gal_image = galsim.ImageF(nx_pixels * nx_stamps , ny_pixels * ny_stamps)
gal_image.setOrigin(0,0) # For my convenience -- switch to C indexing convention.
@@ -184,10 +184,10 @@ def Script1():
logger.info('i,j,x,y,ellip,theta')
for ix in range(nx_stamps):
for iy in range(ny_stamps):
- # The -2's in the next line rather than -1 are to provide a border of
+ # The -2's in the next line rather than -1 are to provide a border of
# 1 pixel between postage stamps
b = galsim.BoundsI(ix*nx_pixels , (ix+1)*nx_pixels -2,
- iy*ny_pixels , (iy+1)*ny_pixels -2)
+ iy*ny_pixels , (iy+1)*ny_pixels -2)
sub_image = gal_image[b]
# Great08 randomized the locations of the two galaxies in each pair,
@@ -200,7 +200,7 @@ def Script1():
# Apply a random orientation:
#theta = rng() * 2. * math.pi * galsim.radians
- theta = rng() * 2. * math.pi
+ theta = rng() * 2. * math.pi
first_in_pair = False
else:
#theta += math.pi/2 * galsim.radians
@@ -273,7 +273,7 @@ def Script2():
- Applied shear is the same for each file
- Noise is poisson using a nominal sky value of 1.e6
"""
- logger = logging.getLogger("Script2")
+ logger = logging.getLogger("Script2")
# Define some parameters we'll use below.
@@ -302,58 +302,74 @@ def Script2():
rng = galsim.UniformDeviate(random_seed)
# Setup the config object
-
config = galsim.Config()
+ # The configuration should set up several top level things:
+ # config.psf defines the PSF
+ # config.pix defines the pixel response
+ # config.gal defines the galaxy
+ # They are all currently required, although eventually we will probably make it so
+ # they are each optional.
+
+ # Each type of profile is specified by a type. e.g. Moffat:
config.psf.type = 'Moffat'
- config.psf.beta.type = 'InputCatalog'
- config.psf.beta.col = 5
+
+ # The various parameters are typically specified as well
+ config.psf.beta = 3.5
+
+ # These parameters do not need to be constant. There are a number of ways to
+ # specify variables that might change from object to object.
+ # In this case, the parameter specification also has a "type".
+ # For now we only have InputCatalog, which means read the value from a catalog:
config.psf.fwhm.type = 'InputCatalog'
+
+ # InputCatalog requires the extra value of which column to use in the catalog:
config.psf.fwhm.col = 6
- config.psf.ellip.type = 'E1E2'
- config.psf.ellip.e1.type = 'InputCatalog'
- config.psf.ellip.e1.col = 7
- config.psf.ellip.e2.type = 'InputCatalog'
- config.psf.ellip.e2.col = 8
- # TODO rename truncationFWHM to something saner, like trunc
- # And probably make it in units of arcsec rather than FWHM.
- config.psf.truncationFWHM.type = 'InputCatalog'
- config.psf.truncationFWHM.col = 9
-
- config.pix.type = 'SquarePixel'
- config.pix.size = pixel_scale
+ # You can also specify both of these on the same line as a single string:
+ config.psf.truncationFWHM = 'InputCatalog col=9'
+
+ # You can even nest string values using angle brackets:
+ config.psf.ellip = 'E1E2 e1=<InputCatalog col=7> e2=<InputCatalog col=8>'
+
+ # If you don't specify a parameter, and there is a reasonable default, then it
+ # will be used instead. If there is no reasonable default, you will get an error.
+ #config.psf.flux = 1 # Unnecessary
+
+ # If you want to use a variable in your string, you can use Python's % notation:
+ config.pix = 'SquarePixel size=%f'%pixel_scale
+
+ # A profile can be the sum of several components, each with its own type and parameters:
config.gal.type = 'Sum'
# TODO: [galsim.Config()]*2 doesn't work, since shallow copies.
# I guess we need a nicer way to initialize this.
config.gal.items = [galsim.Config(), galsim.Config()]
config.gal.items[0].type = 'Exponential'
- config.gal.items[0].half_light_radius.type = 'InputCatalog'
- config.gal.items[0].half_light_radius.col = 10
- config.gal.items[0].ellip.type = 'E1E2'
- config.gal.items[0].ellip.e1.type = 'InputCatalog'
- config.gal.items[0].ellip.e1.col = 11
- config.gal.items[0].ellip.e2.type = 'InputCatalog'
- config.gal.items[0].ellip.e2.col = 12
+ config.gal.items[0].half_light_radius = 'InputCatalog col=10'
+ config.gal.items[0].ellip = 'E1E2 e1=<InputCatalog col=11> e2=<InputCatalog col=12>'
config.gal.items[0].flux = 0.6
config.gal.items[1].type = 'DeVaucouleurs'
- config.gal.items[1].half_light_radius.type = 'InputCatalog'
- config.gal.items[1].half_light_radius.col = 13
- config.gal.items[1].ellip.type = 'E1E2'
- config.gal.items[1].ellip.e1.type = 'InputCatalog'
- config.gal.items[1].ellip.e1.col = 14
- config.gal.items[1].ellip.e2.type = 'InputCatalog'
- config.gal.items[1].ellip.e2.col = 15
+ config.gal.items[1].half_light_radius = 'InputCatalog col=13'
+ config.gal.items[1].ellip = 'E1E2 e1=<InputCatalog col=14> e2=<InputCatalog col=15>'
config.gal.items[1].flux = 0.4
+
+ # When a composite object (like a Sum) has a flux specified, the "flux" values of the
+ # components are taken to be relative fluxes, and the full object's value sets the
+ # overall normalization. If this is omitted, the overall flux is taken to be the
+ # sum of the component fluxes.
config.gal.flux = gal_flux
- config.gal.shift.type = 'DXDY'
- config.gal.shift.dx.type = 'InputCatalog'
- config.gal.shift.dx.col = 16
- config.gal.shift.dy.type = 'InputCatalog'
- config.gal.shift.dy.col = 17
- config.gal.shear.type = 'G1G2'
- config.gal.shear.g1 = gal_g1
- config.gal.shear.g2 = gal_g2
+
+ # An object may have an ellip and a shear, each of which can be specified in terms
+ # of either E1E2 (distortion) or G1G2 (reduced shear).
+ # The only difference between the two is if there is also a rotation specified.
+ # The order of the various modifications are:
+ # - ellip
+ # - rotation
+ # - shear
+ # - shift
+ config.gal.shear = 'G1G2 g1=%f g2=%f'%(gal_g1,gal_g2)
+ config.gal.shift = 'DXDY dx=<InputCatalog col=16> dy=<InputCatalog col=17>'
+
# Read the catalog
input_cat = galsim.io.ReadInputCat(config,cat_file_name)
@@ -398,10 +414,109 @@ def Script2():
# Now write the image to disk.
galsim.fits.writeMulti(all_images, multi_file_name, clobber=True)
- logger.info('Wrote images to multi-extension fits file %r',multi_file_name)
+ logger.info('Wrote images to multi-extension fits file %r',multi_file_name)
+
+ galsim.fits.writeCube(all_images, cube_file_name, clobber=True)
+ logger.info('Wrote image to fits data cube %r',cube_file_name)
+
+ print
+
+# Script 3: Simulations with real galaxies from a catalog
+def Script3():
+ """
+ Make a fits image cube using real COSMOS galaxies from a catalog describing the training
+ sample.
+
+ - The number of images in the cube matches the number of rows in the catalog.
+ - Each image size is computed automatically by GalSim based on the Nyquist size.
+ - Both galaxies and stars.
+ - PSF is a double Gaussian (eventually Kolmogorov), the same for each galaxy.
+ - Galaxies are randomly rotated to remove the imprint of any lensing shears in the COSMOS
+ data.
+ - The same shear is applied to each galaxy.
+ - Noise is poisson using a nominal sky value of 1.e6
+ the noise in the original COSMOS data.
+ """
+ logger = logging.getLogger("Script3")
+
+ # Define some parameters we'll use below.
+
+ cat_file_name = os.path.join('data','real_galaxy_catalog_example.fits')
+ image_dir = os.path.join('data')
+ multi_file_name = os.path.join('output','multi_real.fits')
+ cube_file_name = os.path.join('output','cube_real.fits')
+ psf_file_name = os.path.join('output','psf_script3.fits')
+
+ random_seed = 1512413
+ sky_level = 1.e6 # ADU
+ pixel_scale = 0.15 # arcsec
+ gal_flux = 1.e5 # arbitrary choice, makes nice (not too) noisy images
+ gal_g1 = -0.027 #
+ gal_g2 = 0.031 #
+ psf_inner_fwhm = 0.6 # FWHM of inner Gaussian of the double Gaussian PSF, arcsec
+ psf_outer_fwhm = 2*psf_inner_fwhm
+ psf_inner_fraction = 0.8 # fraction of total PSF flux in the inner Gaussian
+
+ logger.info('Starting multi-object script 3 using:')
+ logger.info(' - real galaxies from catalog %r',cat_file_name)
+ logger.info(' - double Gaussian PSF')
+ logger.info(' - pixel scale = %.2f',pixel_scale)
+ logger.info(' - Applied gravitational shear = (%.3f,%.3f)',gal_g1,gal_g2)
+ logger.info(' - Poisson noise (sky level = %.1e).', sky_level)
+
+ # Initialize the random number generator we will be using.
+ rng = galsim.UniformDeviate(random_seed)
+
+ # Read in galaxy catalog
+ real_galaxy_catalog = galsim.RealGalaxyCatalog(cat_file_name, image_dir)
+ n_gal = real_galaxy_catalog.n
+ logger.info('Read in %d real galaxies from catalog', n_gal)
+
+ ## Make the ePSF
+ # first make the double Gaussian PSF
+ psf = galsim.atmosphere.DoubleGaussian(psf_inner_fraction, 1.0-psf_inner_fraction, fwhm1 =
+ psf_inner_fwhm, fwhm2 = psf_outer_fwhm)
+ # make the pixel response
+ pix = galsim.Pixel(xw = pixel_scale, yw = pixel_scale)
+ # convolve PSF and pixel response function to get the effective PSF (ePSF)
+ epsf = galsim.Convolve(psf, pix)
+ epsf_image = epsf.draw(dx = pixel_scale)
+ # write to file
+ epsf_image.write(psf_file_name, clobber = True)
+ logger.info('Created ePSF and wrote to file %r',psf_file_name)
+
+ # Build the images
+ all_images = []
+ for i in range(n_gal):
+ logger.info('Start work on image %d',i)
+
+ gal = galsim.RealGalaxy(real_galaxy_catalog, index = i)
+ logger.info(' Read in training sample galaxy and PSF from file')
+
+ sim_image = galsim.simReal(gal, epsf, pixel_scale, g1 = gal_g1, g2 = gal_g2,
+ uniform_deviate = rng, target_flux = gal_flux)
+ logger.info(' Made image of galaxy after deconvolution, rotation, shearing, ')
+ logger.info(' convolving with target PSF, and resampling')
+ xsize, ysize = sim_image.array.shape
+ logger.info(' Final image size: (%d, %d)',xsize, ysize)
+
+ # Add Poisson noise
+ sim_image += sky_level
+ sim_image.addNoise(galsim.CCDNoise(rng))
+ sim_image -= sky_level
+ logger.info(' Added Poisson noise')
+
+ # Store that into the list of all images
+ all_images += [sim_image]
+
+ logger.info('Done making images of galaxies')
+
+ # Now write the image to disk.
+ galsim.fits.writeMulti(all_images, multi_file_name, clobber=True)
+ logger.info('Wrote images to multi-extension fits file %r',multi_file_name)
galsim.fits.writeCube(all_images, cube_file_name, clobber=True)
- logger.info('Wrote image to fits data cube %r',cube_file_name)
+ logger.info('Wrote image to fits data cube %r',cube_file_name)
print
@@ -413,7 +528,7 @@ def main(argv):
except Exception as err:
print __doc__
raise err
-
+
# Output files are put in the directory output. Make sure it exists.
if not os.path.isdir('output'):
os.mkdir('output')
@@ -428,12 +543,6 @@ def main(argv):
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: Great08-like image
if scriptNum == 0 or scriptNum == 1:
@@ -443,6 +552,8 @@ def main(argv):
if scriptNum == 0 or scriptNum == 2:
Script2()
+ if scriptNum == 0 or scriptNum == 3:
+ Script3()
if __name__ == "__main__":
main(sys.argv)
Oops, something went wrong.