galsim - 将赤经和赤纬转化为图像

Question

我想读取星系目录的赤经（以小时角为单位）、赤纬（以度为单位）和大小（以弧分为单位），并将它们全部绘制在指定像素大小的大图像中。

我尝试将 ra、dec 和 size 转换为像素来为每个星系创建一个 Bounds 对象，但得到一个错误“BoundsI must be initialized with integer values”。我知道像素必须是整数......

但是有没有办法将大图像以指定的 ra 和 dec 为中心，然后输入每个星系的 ra 和 dec 作为参数来绘制它？

先感谢您！

Answer 1

GalSim 使用图像坐标处理图像边界和位置。将天空上的真实位置（RA，dec）连接到图像坐标的方法是使用 GalSim 中的世界坐标系（WCS）功能。我从您的描述中得知，有一个从 RA/dec 到像素坐标的简单映射（即没有失真）。

因此，基本上，您将设置一个简单的 WCS，定义大图像的 (RA, dec) 中心及其像素比例。然后对于给定的星系（RA，dec），您可以使用 WCS 的“toImage”方法来确定星系应该在大图像上的位置。可以使用该信息构建任何子图像边界。

对于一个简单的世界坐标系的简单示例，您可以查看GalSim存储库中的 demo10。

Answer 2

GalSim 使用 CelestialCoord 类来处理天空中的坐标，并使用许多 WCS 类来处理从像素到天体坐标的转换。

本教程系列中使用 CelestialWCS（使用天体坐标作为其世界坐标系的 WCS 类的基类）的两个演示是 demo11 和 demo13。所以你可能想看看它们。但是，没有一个人做的事情与您正在做的事情非常接近。

所以这是一个或多或少符合您所描述的脚本。

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')