Translating right ascension and declination onto image

Question

I want to read in the right ascension (in hour angles), declination (in degrees) and size (in arcmin) of a catalogue of galaxies and draw all of them in a large image of specified pixel size.

I tried converting the ra, dec and size into pixels to create a Bounds object for each galaxy, but get an error that "BoundsI must be initialized with integer values." I understand that pixels have to be integers...

But is there a way to center the large image at a specified ra and dec, then input the ra and dec of each galaxy as parameters to draw it in?

Thank you in advance!

Answer 1

GalSim uses the CelestialCoord class to handle coordinates in the sky and any of a number of WCS classes to handle the conversion from pixels to celestial coordinates.

The two demos in the tutorial series that use a CelestialWCS (the base class for WCS classes that use celestial coordinates for their world coordinate system) are demo11 and demo13. So you might want to take a look at them. However, neither one does something very close to what you're doing.

So here's a script that more or less does what you described.

import galsim
import numpy

# Make some random input data so we can run this.
# You would use values from your input catalog.
ngal = 20
numpy.random.seed(123)
ra = 15 + 0.02*numpy.random.random( (ngal) )    # hours
dec = -34 + 0.3*numpy.random.random( (ngal) )   # degrees
size = 0.1 * numpy.random.random( (ngal) )      # arcmin
e1 = 0.5 * numpy.random.random( (ngal) ) - 0.25
e2 = 0.5 * numpy.random.random( (ngal) ) - 0.25

# arcsec is usually the more natural units for sizes, so let's
# convert to that here to make things simpler later.
# There are options throughout GalSim to do things in different
# units, such as arcmin, but arcsec is the default, so it will
# be simpler if we don't have to worry about that.
size *= 60  # size now in arcsec

# Some plausible location at which to center the image.
# Note that we are now attaching the right units to these
# so GalSim knows what angle they correspond to.
cen_ra = numpy.mean(ra) * galsim.hours
cen_dec = numpy.mean(dec) * galsim.degrees

# GalSim uses CelestialCoord to handle celestial coordinates.
# It knows how to do all the correct spherical geometry calculations.
cen_coord = galsim.CelestialCoord(cen_ra, cen_dec)
print 'cen_coord = ',cen_coord.ra.hms(), cen_coord.dec.dms()

# Define some reasonable pixel size.
pixel_scale = 0.4  # arcsec / pixel

# Make the full image of some size.
# Powers of two are typical, but not required.
image_size = 2048
image = galsim.Image(image_size, image_size)

# Define the WCS we'll use to connect pixels to celestial coords.
# For real data, this would usually be read from the FITS header.
# Here, we'll need to make our own.  The simplest one that properly
# handles celestial coordinates is TanWCS.  It first goes from
# pixels to a local tangent plane using a linear affine transformation.
# Then it projects that tangent plane into the spherical sky coordinates.
# In our case, we can just let the affine transformation be a uniform
# square pixel grid with its origin at the center of the image.
affine_wcs = galsim.PixelScale(pixel_scale).affine().withOrigin(image.center())
wcs = galsim.TanWCS(affine_wcs, world_origin=cen_coord)
image.wcs = wcs  # Tell the image to use this WCS

for i in range(ngal):
    # Get the celestial coord of the galaxy
    coord = galsim.CelestialCoord(ra[i]*galsim.hours, dec[i]*galsim.degrees)
    print 'gal coord = ',coord.ra.hms(), coord.dec.dms()

    # Where is it in the image?
    image_pos = wcs.toImage(coord)
    print 'position in image = ',image_pos

    # Make some model of the galaxy.
    flux = size[i]**2 * 1000  # Make bigger things brighter...
    gal = galsim.Exponential(half_light_radius=size[i], flux=flux)
    gal = gal.shear(e1=e1[i],e2=e2[i])

    # Pull out a cutout around where we want the galaxy to be.
    # The bounds needs to be in integers.
    # The fractional part of the position will go into offset when we draw.
    ix = int(image_pos.x)
    iy = int(image_pos.y)
    bounds = galsim.BoundsI(ix-64, ix+64, iy-64, iy+64)

    # This might be (partially) off the full image, so get the overlap region.
    bounds = bounds & image.bounds
    if not bounds.isDefined():
        print '    This galaxy is completely off the image.'
        continue

    # This is the portion of the full image where we will draw.  If you try to
    # draw onto the full image, it will use a lot of memory, but if you go too
    # small, you might see artifacts at the edges.  You might need to 
    # experiment a bit with what is a good size cutout.
    sub_image = image[bounds]

    # Draw the galaxy.  
    # GalSim by default will center the object at the "true center" of the
    # image.  We actually want it centered at image_pos, so provide the
    # difference as the offset parameter.
    # Also, the default is to overwrite the image.  But we want to add to
    # the existing image in case galaxies overlap.  Hence add_to_image=True
    gal.drawImage(image=sub_image, offset=image_pos - sub_image.trueCenter(),
                  add_to_image=True)

# Probably want to add a little noise...
image.addNoise(galsim.GaussianNoise(sigma=0.5))

# Write to a file.
image.write('output.fits')

Answer 2

GalSim deals with image bounds and locations using image coordinates. The way to connect true positions on the sky (RA, dec) into image coordinates is using the World Coordinate System (WCS) functionality in GalSim. I gather from your description that there is a simple mapping from RA/dec into pixel coordinates (i.e., there are no distortions).

So basically, you would set up a simple WCS defining the (RA, dec) center of the big image and its pixel scale. Then for a given galaxy (RA, dec), you can use the "toImage" method of the WCS to figure out where on the big image the galaxy should live. Any subimage bounds can be constructed using that information.

For a simple example with a trivial world coordinate system, you can check out demo10 in the GalSim repository.