World Coordinate Systems

The World Coordinate System (or WCS) is the traditional term for the mapping from pixel coordinates to the coordinate system on the sky. (I know, the world’s down here, and the sky’s up there, so you’d think it would be reversed, but that’s the way it goes. Astronomy is full of terms that don’t quite make sense when you look at them too closely.)

There are two kinds of world coordinates that we use here:

  • Celestial coordinates are defined in terms of right ascension (RA) and declination (Dec). They are a spherical coordinate system on the sky, akin to longitude and latitude on Earth. cf. http://en.wikipedia.org/wiki/Celestial_coordinate_system

  • Euclidean coordinates are defined relative to a tangent plane projection of the sky. If you imagine the sky coordinates on an actual sphere with a particular radius, then the tangent plane is tangent to that sphere. We use the labels (u,v) for the coordinates in this system, where +v points north and +u points west. (Yes, west, not east. As you look up into the sky, if north is up, then west is to the right.)

The classes in this file provide a mapping from image coordinates (in pixels) to one of these two kinds of world coordinates. We use the labels (x,y) for the image coordinates.

WCS Base Classes

class galsim.BaseWCS[source]

The base class for all other kinds of WCS transformations.

All the functions the user will typically need are defined here. Most subclasses just define helper functions to implement each particular WCS definition. So this base class defines the common interface for all WCS classes.

There are several types of WCS classes that we implement. The basic class hierarchy is:

`BaseWCS`
    --- `EuclideanWCS`
            --- `UniformWCS`
                    --- `LocalWCS`
    --- `CelestialWCS`

These base classes are not constructible. They do not have __init__ defined.

  1. LocalWCS classes are those which really just define a pixel size and shape. They implicitly have the origin in image coordinates correspond to the origin in world coordinates. They are primarily designed to handle local transformations at the location of a single galaxy, where it should usually be a good approximation to consider the pixel shape to be constant over the size of the galaxy.

    Currently we define the following LocalWCS classes:

    - `PixelScale`
- `ShearWCS`
- `JacobianWCS`

  2. UniformWCS classes have a constant pixel size and shape, but they have an arbitrary origin in both image coordinates and world coordinates. A LocalWCS class can be turned into a non-local UniformWCS class when an image has its bounds changed, e.g. by the commands Image.setCenter, Image.setOrigin or Image.shift.

    Currently we define the following non-local, UniformWCS classes:

    - `OffsetWCS`
- `OffsetShearWCS`
- `AffineTransform`

  3. EuclideanWCS classes use a regular Euclidean coordinate system for the world coordinates, using PositionD for the world positions. We use the notation (u,v) for the world coordinates and (x,y) for the image coordinates.

    Currently we define the following non-uniform, EuclideanWCS class:

    - `UVFunction`

  4. CelestialWCS classes are defined with their world coordinates on the celestial sphere in terms of right ascension (RA) and declination (Dec). The pixel size and shape are always variable. We use CelestialCoord for the world coordinates, which helps facilitate the spherical trigonometry that is sometimes required.

    Currently we define the following CelestialWCS classes: (All but the first are defined in the file fitswcs.py.)

    • RaDecFunction

    • AstropyWCS – requires astropy.wcs python module to be installed

    • PyAstWCS – requires starlink.Ast python module to be installed

    • WcsToolsWCS – requires wcstools command line functions to be installed

    • GSFitsWCS – native code, but has less functionality than the above

There are also a few factory functions in fitswcs.py intended to act like class initializers:

  • FitsWCS tries to read a fits file using one of the above classes and returns an instance of whichever one it found was successful. It should always be successful, since its final attempt uses AffineTransform, which has reasonable defaults when the WCS key words are not in the file, but of course this will only be a very rough approximation of the true WCS.

  • TanWCS constructs a simple tangent plane projection WCS directly from the projection parameters instead of from a fits header.

  • FittedSIPWCS constructs a TAN-SIP WCS by fitting to a list of reference celestial and image coordinates.

Some things you can do with a WCS instance:

  • Convert positions between image coordinates and world coordinates (sometimes referred to as sky coordinates):

    >>> world_pos = wcs.toWorld(image_pos)
>>> image_pos = wcs.toImage(world_pos)

    Note: the transformation from world to image coordinates is not guaranteed to be implemented. If it is not implemented for a particular WCS class, a NotImplementedError will be raised.

    The image_pos parameter should be a PositionD. However, world_pos will be a CelestialCoord if the transformation is in terms of celestial coordinates (if wcs.isCelestial() == True). Otherwise, it will be a PositionD as well.

  • Convert a GSObject that is defined in world coordinates to the equivalent profile defined in terms of image coordinates (or vice versa):

    >>> image_profile = wcs.toImage(world_profile)
>>> world_profile = wcs.toWorld(image_profile)

    For non-uniform WCS types (for which wcs.isUniform() == False), these need either an image_pos or world_pos parameter to say where this conversion should happen:

    >>> image_profile = wcs.toImage(world_profile, image_pos=image_pos)

  • Construct a local linear approximation of a WCS at a given location:

    >>> local_wcs = wcs.local(image_pos = image_pos)
>>> local_wcs = wcs.local(world_pos = world_pos)

    If wcs.toWorld(image_pos) is not implemented for a particular WCS class, then a NotImplementedError will be raised if you pass in a world_pos argument.

    The returned local_wcs is usually a JacobianWCS instance, but see the doc string for local for more details.

  • Construct a full affine approximation of a WCS at a given location:

    >>> affine_wcs = wcs.affine(image_pos = image_pos)
>>> affine_wcs = wcs.affine(world_pos = world_pos)

    This preserves the transformation near the location of image_pos, but it is linear, so the transformed values may not agree as you get farther from the given point.

    The returned affine_wcs is always an AffineTransform instance.

  • Get some properties of the pixel size and shape:

    >>> area = local_wcs.pixelArea()
>>> min_linear_scale = local_wcs.minLinearScale()
>>> max_linear_scale = local_wcs.maxLinearScale()
>>> jac = local_wcs.jacobian()
>>> # Use jac.dudx, jac.dudy, jac.dvdx, jac.dvdy

    Non-uniform WCS types also have these functions, but for them, you must supply either image_pos or world_pos. So the following are equivalent:

    >>> area = wcs.pixelArea(image_pos)
>>> area = wcs.local(image_pos).pixelArea()

  • Query some overall attributes of the WCS transformation:

    >>> wcs.isLocal()       # is this a local WCS?
>>> wcs.isUniform()     # does this WCS have a uniform pixel size/shape?
>>> wcs.isCelestial()   # are the world coordinates on the celestial sphere?
>>> wcs.isPixelScale()  # is this either a PixelScale or an OffsetWCS?
affine(image_pos=None, world_pos=None, color=None)[source]

Return the local AffineTransform of the WCS at a given point.

This returns a linearized version of the current WCS at a given point. It returns an AffineTransform that is locally approximately the same as the WCS in the vicinity of the given point.

It is similar to jacobian(), except that this preserves the offset information between the image coordinates and world coordinates rather than setting both origins to (0,0). Instead, the image origin is taken to be image_pos.

For non-celestial coordinate systems, the world origin is taken to be wcs.toWorld(image_pos). In fact, wcs.affine(image_pos) is really just shorthand for:

>>> wcs.jacobian(image_pos).withOrigin(image_pos, wcs.toWorld(image_pos))

For celestial coordinate systems, there is no well-defined choice for the origin of the Euclidean world coordinate system. So we just take (u,v) = (0,0) at the given position. So, wcs.affine(image_pos) is equivalent to:

>>> wcs.jacobian(image_pos).withOrigin(image_pos)

You can use the returned AffineTransform to access the relevant values of the 2x2 Jacobian matrix and the origins directly:

>>> affine = wcs.affine(image_pos)
>>> x,y = np.meshgrid(np.arange(0,32,1), np.arange(0,32,1))
>>> u = affine.dudx * (x-affine.x0) + jac.dudy * (y-affine.y0) + affine.u0
>>> v = affine.dvdx * (x-affine.x0) + jac.dvdy * (y-affine.y0) + affine.v0
>>> # ... use u,v values to work directly in sky coordinates.

As usual, you may provide either image_pos or world_pos as you prefer to specify the location at which to approximate the WCS.

Parameters:

  • image_pos – The image coordinate position (for non-uniform WCS types)

  • world_pos – The world coordinate position (for non-uniform WCS types)

  • color – For color-dependent WCS’s, the color term for which to evaluate the local affine transform. [default: None]

Returns:

an AffineTransform instance

fixColor(color)[source]

Fix the color to a particular value.

This changes a color-dependent WCS into the corresponding color-independent WCS for the given color.

Parameters:

color – The value of the color term to use.

Returns:

the new color-independent WCS

isCelestial()[source]

Return whether the world coordinates are CelestialCoord (i.e. ra,dec).

wcs.isCelestial() is shorthand for isinstance(wcs, galsim.CelestialWCS).

isLocal()[source]

Return whether the WCS transformation is a local, linear approximation.

wcs.isLocal() is shorthand for isinstance(wcs, galsim.LocalWCS).

isPixelScale()[source]

Return whether the WCS transformation is a simple PixelScale or OffsetWCS.

These are the simplest two WCS transformations. PixelScale is local and OffsetWCS is non-local. If an Image has one of these WCS transformations as its WCS, then im.scale works to read and write the pixel scale. If not, im.scale will raise a TypeError exception.

wcs.isPixelScale() is shorthand for isinstance(wcs, (galsim.PixelScale, galsim.OffsetWCS)).

isUniform()[source]

Return whether the pixels in this WCS have uniform size and shape.

wcs.isUniform() is shorthand for isinstance(wcs, galsim.UniformWCS).

jacobian(image_pos=None, world_pos=None, color=None)[source]

Return the local JacobianWCS of the WCS at a given point.

This is basically the same as local(), but the return value is guaranteed to be a JacobianWCS, which can be useful in some situations, since you can access the values of the 2x2 Jacobian matrix directly:

>>> jac = wcs.jacobian(image_pos)
>>> x,y = np.meshgrid(np.arange(0,32,1), np.arange(0,32,1))
>>> u = jac.dudx * x + jac.dudy * y
>>> v = jac.dvdx * x + jac.dvdy * y
>>> # ... use u,v values to work directly in world coordinates.

If you do not need the extra functionality, then you should use local() instead, since it may be more efficient.

Parameters:

  • image_pos – The image coordinate position (for non-uniform WCS types)

  • world_pos – The world coordinate position (for non-uniform WCS types)

  • color – For color-dependent WCS’s, the color term for which to evaluate the local jacobian. [default: None]

Returns:

a JacobianWCS instance.

local(image_pos=None, world_pos=None, color=None)[source]

Return the local linear approximation of the WCS at a given point.

Parameters:

  • image_pos – The image coordinate position (for non-uniform WCS types)

  • world_pos – The world coordinate position (for non-uniform WCS types)

  • color – For color-dependent WCS’s, the color term for which to evaluate the local WCS. [default: None]

Returns:

a LocalWCS instance.

makeSkyImage(image, sky_level, color=None)[source]

Make an image of the sky, correctly accounting for the pixel area, which might be variable over the image.

Note: This uses finite differences of the wcs mapping to calculate the area of each

pixel in world coordinates. It is usually pretty accurate everywhere except within a few arcsec of the north or south poles.

Parameters:

  • image – The image onto which the sky values will be put.

  • sky_level – The sky level in ADU/arcsec^2 (or whatever your world coordinate system units are, if not arcsec).

  • color – For color-dependent WCS’s, the color term to use for making the sky image. [default: None]

maxLinearScale(image_pos=None, world_pos=None, color=None)[source]

Return the maximum linear scale of the transformation in any direction.

This is basically the semi-major axis of the Jacobian. Sometimes you need a linear scale size for some calculation. This function returns the largest scale in any direction. The function minLinearScale() returns the smallest.

For non-uniform WCS transforms, you must provide either image_pos or world_pos to say where the pixel is located.

Parameters:

  • image_pos – The image coordinate position (for non-uniform WCS types)

  • world_pos – The world coordinate position (for non-uniform WCS types)

  • color – For color-dependent WCS’s, the color term for which to evaluate the scale. [default: None]

Returns:

the maximum pixel area in any direction in arcsec.

minLinearScale(image_pos=None, world_pos=None, color=None)[source]

Return the minimum linear scale of the transformation in any direction.

This is basically the semi-minor axis of the Jacobian. Sometimes you need a linear scale size for some calculation. This function returns the smallest scale in any direction. The function maxLinearScale() returns the largest.

For non-uniform WCS transforms, you must provide either image_pos or world_pos to say where the pixel is located.

Parameters:

  • image_pos – The image coordinate position (for non-uniform WCS types)

  • world_pos – The world coordinate position (for non-uniform WCS types)

  • color – For color-dependent WCS’s, the color term for which to evaluate the scale. [default: None]

Returns:

the minimum pixel area in any direction in arcsec.

pixelArea(image_pos=None, world_pos=None, color=None)[source]

Return the area of a pixel in arcsec**2 (or in whatever units you are using for world coordinates if it is a EuclideanWCS).

For non-uniform WCS transforms, you must provide either image_pos or world_pos to say where the pixel is located.

Parameters:

  • image_pos – The image coordinate position (for non-uniform WCS types)

  • world_pos – The world coordinate position (for non-uniform WCS types)

  • color – For color-dependent WCS’s, the color term for which to evaluate the pixel area. [default: None]

Returns:

the pixel area in arcsec**2.

posToImage(world_pos, color=None)[source]

Convert a position from world coordinates to image coordinates.

This is equivalent to wcs.toImage(world_pos).

Parameters:

  • world_pos – The world coordinate position

  • color – For color-dependent WCS’s, the color term to use. [default: None]

posToWorld(image_pos, color=None, **kwargs)[source]

Convert a position from image coordinates to world coordinates.

This is equivalent to wcs.toWorld(image_pos).

Parameters:

  • image_pos – The position in image coordinates

  • color – For color-dependent WCS’s, the color term to use. [default: None]

  • project_center – (Only valid for CelestialWCS) A CelestialCoord to use for projecting the result onto a tangent plane world system rather than returning a CelestialCoord. [default: None]

  • projection – If project_center != None, the kind of projection to use. See CelestialCoord.project for the valid options. [default: ‘gnomonic’]

Returns:

world_pos

profileToImage(world_profile, image_pos=None, world_pos=None, color=None, flux_ratio=1.0, offset=(0, 0))[source]

Convert a profile from world coordinates to image coordinates.

This is equivalent to wcs.toImage(world_profile, ...).

Parameters:

  • world_profile – The profile in world coordinates to transform.

  • image_pos – The image coordinate position (for non-uniform WCS types)

  • world_pos – The world coordinate position (for non-uniform WCS types)

  • color – For color-dependent WCS’s, the color term to use. [default: None]

  • flux_ratio – An optional flux scaling to be applied at the same time. [default: 1]

  • offset – An optional offset to be applied at the same time. [default: 0,0]

profileToWorld(image_profile, image_pos=None, world_pos=None, color=None, flux_ratio=1.0, offset=(0, 0))[source]

Convert a profile from image coordinates to world coordinates.

This is equivalent to wcs.toWorld(image_profile, ...).

Parameters:

  • image_profile – The profile in image coordinates to transform.

  • image_pos – The image coordinate position (for non-uniform WCS types)

  • world_pos – The world coordinate position (for non-uniform WCS types)

  • color – For color-dependent WCS’s, the color term to use. [default: None]

  • flux_ratio – An optional flux scaling to be applied at the same time. [default: 1]

  • offset – An optional offset to be applied at the same time. [default: 0,0]

shearToImage(world_shear, image_pos=None, world_pos=None, color=None)[source]

Convert a shear measured in world coordinates to image coordinates

This is equivalent to wcs.toImage(shear, ...).

This reverses the process described in shearToWorld.

Parameters:

  • world_shear – The shear in world coordinates to convert.

  • image_pos – The image coordinate position (for non-uniform WCS types)

  • world_pos – The world coordinate position (for non-uniform WCS types)

  • color – For color-dependent WCS’s, the color term to use. [default: None]

Returns:

The shear in image coordinates.

Return type:

image_shear

shearToWorld(image_shear, image_pos=None, world_pos=None, color=None)[source]

Convert a shear measured in image coordinates to world coordinates

This is equivalent to wcs.toWorld(shear, ...).

Specifically, the input shear is taken to be defined such that +e1 means elongated along the x-axis, -e1 is along the y-axis, +e2 is along the y=x line, and -e2 is along y=-x.

If the WCS is a EuclideanWCS, then the returned shear is converted to the equivalent shear in u-v coordinates. I.e. +e1 is along the u-axis, etc.

If the WCS is a CelestialWCS, then the returned shear is converted to sky coordinates where North is up and West is to the right (as appropriate for looking up into the sky from Earth). I.e. +e1 is E-W, -e1 is N-S, +e2 is NW-SE, and -e2 is NE-SW. This is the convention used by many, but not all, codes and catalogs that use or report shears on the spherical sky.

Parameters:

  • image_shear – The shear in image coordinates to convert.

  • image_pos – The image coordinate position (for non-uniform WCS types)

  • world_pos – The world coordinate position (for non-uniform WCS types)

  • color – For color-dependent WCS’s, the color term to use. [default: None]

Returns:

The shear in world coordinates.

Return type:

world_shear

shiftOrigin(origin, world_origin=None, color=None)[source]

Shift the origin of the current WCS function, returning the new WCS.

This function creates a new WCS instance (always a non-local WCS) that shifts the origin by the given amount. In other words, it treats the image position origin the same way the current WCS treats (x,y) = (0,0).

If the current WCS is a local WCS, this essentially declares where on the image you want the origin of the world coordinate system to be. i.e. where is (u,v) = (0,0). So, for example, to set a WCS that has a constant pixel size with the world coordinates centered at the center of an image, you could write:

>>> wcs = galsim.PixelScale(scale).shiftOrigin(im.center)

This is equivalent to the following:

>>> wcs = galsim.OffsetWCS(scale, origin=im.center)

For non-local WCS types, the origin defines the location in the image coordinate system should mean the same thing as (x,y) = (0,0) does for the current WCS. The following example should work regardless of what kind of WCS this is:

>>> world_pos1 = wcs.toWorld(PositionD(0,0))
>>> wcs2 = wcs.shiftOrigin(new_origin)
>>> world_pos2 = wcs2.toWorld(new_origin)
>>> # world_pos1 should be equal to world_pos2

Furthermore, if the current WCS is a EuclideanWCS (wcs.isCelestial() == False) you may also provide a world_origin argument which defines what (u,v) position you want to correspond to the new origin. Continuing the previous example:

>>> wcs3 = wcs.shiftOrigin(new_origin, new_world_origin)
>>> world_pos3 = wcs3.toWorld(new_origin)
>>> # world_pos3 should be equal to new_world_origin
Parameters:

  • origin – The image coordinate position to use for what is currently treated as (0,0).

  • world_origin – The world coordinate position to use at the origin. Only valid if wcs.isCelestial() == False. [default: None]

  • color – For color-dependent WCS’s, the color term to use in the connection between the current origin and world_origin. [default: None]

Returns:

the new shifted WCS

toImage(*args, **kwargs)[source]

Convert from world coordinates to image coordinates

There are essentially three overloaded versions of this function here.

  1. The first converts a position from world coordinates to image coordinates. If the WCS is a EuclideanWCS, the argument may be either a PositionD or PositionI argument. If it is a CelestialWCS, then the argument must be a CelestialCoord. It returns the corresponding position in image coordinates as a PositionD:

    >>> image_pos = wcs.toImage(world_pos)

    Equivalent to posToImage.

  2. The second is nearly the same, but takes either u and v values or ra and dec values (depending on the kind of wcs being used) directly and returns x and y values. For this version, the inputs may be numpy arrays, in which case the returned values are also numpy arrays:

    >>> x, y = wcs.toImage(u, v)                 # For EuclideanWCS types
>>> x, y = wcs.toImage(ra, dec, units=units) # For CelestialWCS types

    Equivalent to uvToxy or radecToxy.

  3. The third converts a surface brightness profile (a GSObject) from world coordinates to image coordinates, returning the profile in image coordinates as a new GSObject. For non-uniform WCS transforms, you must provide either image_pos or world_pos to say where the profile is located so the right transformation can be performed. And optionally, you may provide a flux scaling to be performed at the same time:

    >>> image_profile = wcs.toImage(world_profile, image_pos=None, world_pos=None,
                                flux_ratio=1, offset=(0,0))

    Equivalent to profileToImage.

  4. The fourth converts a shear measurement (a Shear) from image coordinates to world coordinates. As above, for non-uniform WCS transforms, you must provide either image_pos or world_pos. If the wcs is a CelestialWCS, then the returned shear follows the convention used by TreeCorr (among others) for shear on a sphere, namely in the local coordinates where north is up and west is to the right.

    >>> world_shear = wcs.toWorld(image_shear, image_pos=None, world_pos=None)

    Equivalent to shearToImage.

toWorld(*args, **kwargs)[source]

Convert from image coordinates to world coordinates.

There are essentially four overloaded versions of this function here.

  1. The first converts a Position from image coordinates to world coordinates. It returns the corresponding position in world coordinates as a PositionD if the WCS is a EuclideanWCS, or a CelestialCoord if it is a CelestialWCS:

    >>> world_pos = wcs.toWorld(image_pos)

    Equivalent to posToWorld.

  2. The second is nearly the same, but takes x and y values directly and returns either u, v or ra, dec, depending on the kind of wcs being used. For this version, x and y may be numpy arrays, in which case the returned values are also numpy arrays:

    >>> u, v = wcs.toWorld(x, y)                 # For EuclideanWCS types
>>> ra, dec = wcs.toWorld(x, y, units=units) # For CelestialWCS types

    Equivalent to xyTouv or xyToradec.

  3. The third converts a surface brightness profile (a GSObject) from image coordinates to world coordinates, returning the profile in world coordinates as a new GSObject. For non-uniform WCS transforms, you must provide either image_pos or world_pos to say where the profile is located, so the right transformation can be performed. And optionally, you may provide a flux scaling to be performed at the same time:

    >>> world_profile = wcs.toWorld(image_profile, image_pos=None, world_pos=None,
                                flux_ratio=1, offset=(0,0))

    Equivalent to profileToWorld.

  4. The fourth converts a shear measurement (a Shear) from image coordinates to world coordinates. As above, for non-uniform WCS transforms, you must provide either image_pos or world_pos. If the wcs is a CelestialWCS, then the returned shear follows the convention used by TreeCorr (among others) for shear on a sphere, namely in the local coordinates where north is up and west is to the right.

    >>> world_shear = wcs.toWorld(image_shear, image_pos=None, world_pos=None)

    Equivalent to shearToWorld.

writeToFitsHeader(header, bounds)[source]

Write this WCS function to a FITS header.

This is normally called automatically from within the galsim.fits.write() function.

The code will attempt to write standard FITS WCS keys so that the WCS will be readable by other software (e.g. ds9). It may not be able to do so accurately, in which case a linearized version will be used instead. (Specifically, it will use the local affine transform with respect to the image center.)

However, this is not necessary for the WCS to survive a round trip through the FITS header, as it will also write GalSim-specific key words that should allow it to reconstruct the WCS correctly.

Parameters:

  • header – A FitsHeader (or dict-like) object to write the data to.

  • bounds – The bounds of the image.

class galsim.wcs.CelestialWCS[source]

Bases: BaseWCS

A CelestialWCS is a BaseWCS whose world coordinates are on the celestial sphere. We use the CelestialCoord class for the world coordinates.

radecToxy(ra, dec, units, color=None)[source]

Convert ra,dec from world coordinates to image coordinates.

This is equivalent to wcs.toImage(ra,dec, units=units).

It is also equivalent to wcs.posToImage(galsim.CelestialCoord(ra * units, dec * units)) when ra and dec are scalars; however, this routine allows ra and dec to be numpy arrays, in which case, the calculation will be vectorized, which is often much faster than using the pos interface.

Parameters:

  • ra – The ra value(s) in world coordinates

  • dec – The dec value(s) in world coordinates

  • units – The units to use for the input ra, dec values.

  • color – For color-dependent WCS’s, the color term to use. [default: None]

property x0

The x coordinate of self.origin.

xyToradec(x, y, units=None, color=None)[source]

Convert x,y from image coordinates to world coordinates.

This is equivalent to wcs.toWorld(x,y, units=units).

It is also equivalent to wcs.posToWorld(galsim.PositionD(x,y)).rad when x and y are scalars if units is ‘radians’; however, this routine allows x and y to be numpy arrays, in which case, the calculation will be vectorized, which is often much faster than using the pos interface.

Parameters:

  • x – The x value(s) in image coordinates

  • y – The y value(s) in image coordinates

  • units – (Only valid for CelestialWCS, in which case it is required) The units to use for the returned ra, dec values.

  • color – For color-dependent WCS’s, the color term to use. [default: None]

Returns:

ra, dec

property y0

The y coordinate of self.origin.

class galsim.wcs.EuclideanWCS[source]

Bases: BaseWCS

A EuclideanWCS is a BaseWCS whose world coordinates are on a Euclidean plane. We usually use the notation (u,v) to refer to positions in world coordinates, and they use the class PositionD.

property u0

The x component of self.world_origin (aka u).

uvToxy(u, v, color=None)[source]

Convert u,v from world coordinates to image coordinates.

This is equivalent to wcs.toWorld(u,v).

It is also equivalent to wcs.posToImage(galsim.PositionD(u,v)) when u and v are scalars; however, this routine allows u and v to be numpy arrays, in which case, the calculation will be vectorized, which is often much faster than using the pos interface.

Parameters:

  • u – The u value(s) in world coordinates

  • v – The v value(s) in world coordinates

  • color – For color-dependent WCS’s, the color term to use. [default: None]

property v0

The y component of self.world_origin (aka v).

property x0

The x component of self.origin.

xyTouv(x, y, color=None)[source]

Convert x,y from image coordinates to world coordinates.

This is equivalent to wcs.toWorld(x,y).

It is also equivalent to wcs.posToWorld(galsim.PositionD(x,y)) when x and y are scalars; however, this routine allows x and y to be numpy arrays, in which case, the calculation will be vectorized, which is often much faster than using the pos interface.

Parameters:

  • x – The x value(s) in image coordinates

  • y – The y value(s) in image coordinates

  • color – For color-dependent WCS’s, the color term to use. [default: None]

Returns:

ra, dec

property y0

The y component of self.origin.

class galsim.wcs.UniformWCS[source]

Bases: EuclideanWCS

A UniformWCS is a EuclideanWCS which has a uniform pixel size and shape.

inverse()[source]

Return the inverse transformation, i.e. the transformation that swaps the roles of the “image” and “world” coordinates.

class galsim.wcs.LocalWCS[source]

Bases: UniformWCS

A LocalWCS is a UniformWCS in which (0,0) in image coordinates is at the same place as (0,0) in world coordinates

property origin

The image coordinate position to use as the origin.

withOrigin(origin, world_origin=None, color=None)[source]

Recenter the current WCS function at a new origin location, returning the new WCS.

This function creates a new WCS instance (a non-local WCS) with the same local behavior as the current WCS, but with the given origin. In other words, you are declaring where on the image you want the new origin of the world coordinate system to be. i.e. where is (u,v) = (0,0).

So, for example, to set a WCS that has a constant pixel size with the world coordinates centered at the center of an image, you could write:

>>> wcs = galsim.PixelScale(scale).withOrigin(im.center)

This is equivalent to the following:

>>> wcs = galsim.OffsetWCS(scale, origin=im.center)

You may also provide a world_origin argument which defines what (u,v) position you want to correspond to the new origin.

>>> wcs2 = galsim.PixelScale(scale).withOrigin(new_origin, new_world_origin)
>>> world_pos2 = wcs2.toWorld(new_origin)
>>> assert world_pos2 == new_world_origin

Note

This is equivalent to shiftOrigin, but for for local WCS’s, the shift is also the new location of the origin, so withOrigin is a convenient alternate name for this action. Indeed the shiftOrigin function used to be named withOrigin, but that name was confusing for non-local WCS’s, as the action in that case is really shifting the origin, not setting the new value.

Parameters:

  • origin – The image coordinate position to use as the origin.

  • world_origin – The world coordinate position to use as the origin. [default: None]

  • color – For color-dependent WCS’s, the color term to use in the connection between the current origin and world_origin. [default: None]

Returns:

the new recentered WCS

property world_origin

The world coordinate position to use as the origin.

Euclidean WCS’s

class galsim.PixelScale(scale)[source]

Bases: LocalWCS

This is the simplest possible WCS transformation. It only involves a unit conversion from pixels to arcsec (or whatever units you want to take for your world coordinate system).

The conversion functions are:

u = x * scale v = y * scale

A PixelScale is initialized with the command:

>>> wcs = galsim.PixelScale(scale)
Parameters:

scale – The pixel scale, typically in units of arcsec/pixel.

property scale

The pixel scale

class galsim.OffsetWCS(scale, origin=None, world_origin=None)[source]

Bases: UniformWCS

This WCS is similar to PixelScale, except the origin is not necessarily (0,0) in both the image and world coordinates.

The conversion functions are:

u = (x-x0) * scale + u0 v = (y-y0) * scale + v0

An OffsetWCS is initialized with the command:

>>> wcs = galsim.OffsetWCS(scale, origin=None, world_origin=None)
Parameters:

  • scale – The pixel scale, typically in units of arcsec/pixel.

  • origin – Optional origin position for the image coordinate system. If provided, it should be a PositionD or PositionI. [default: PositionD(0., 0.)]

  • world_origin – Optional origin position for the world coordinate system. If provided, it should be a PositionD. [default: galsim.PositionD(0., 0.)]

property origin

The image coordinate position to use as the origin.

property scale

The pixel scale.

property world_origin

The world coordinate position to use as the origin.

class galsim.ShearWCS(scale, shear)[source]

Bases: LocalWCS

This WCS is a uniformly sheared coordinate system.

The shear is given as the shape that a round object has when observed in image coordinates.

The conversion functions in terms of (g1,g2) are therefore:

x = (u + g1 u + g2 v) / scale / sqrt(1-g1**2-g2**2) y = (v - g1 v + g2 u) / scale / sqrt(1-g1**2-g2**2)

or, writing this in the usual way of (u,v) as a function of (x,y):

u = (x - g1 x - g2 y) * scale / sqrt(1-g1**2-g2**2) v = (y + g1 y - g2 x) * scale / sqrt(1-g1**2-g2**2)

A ShearWCS is initialized with the command:

>>> wcs = galsim.ShearWCS(scale, shear)
Parameters:

  • scale – The pixel scale, typically in units of arcsec/pixel.

  • shear – The shear, which should be a Shear instance.

The Shear transformation conserves object area, so if the input scale == 1 then the transformation represented by the ShearWCS will conserve object area also.

property scale

The pixel scale.

property shear

The applied Shear.

class galsim.OffsetShearWCS(scale, shear, origin=None, world_origin=None)[source]

Bases: UniformWCS

This WCS is a uniformly sheared coordinate system with image and world origins that are not necessarily coincident.

The conversion functions are:

x = ( (1+g1) (u-u0) + g2 (v-v0) ) / scale / sqrt(1-g1**2-g2**2) + x0 y = ( (1-g1) (v-v0) + g2 (u-u0) ) / scale / sqrt(1-g1**2-g2**2) + y0

u = ( (1-g1) (x-x0) - g2 (y-y0) ) * scale / sqrt(1-g1**2-g2**2) + u0 v = ( (1+g1) (y-y0) - g2 (x-x0) ) * scale / sqrt(1-g1**2-g2**2) + v0

An OffsetShearWCS is initialized with the command:

>>> wcs = galsim.OffsetShearWCS(scale, shear, origin=None, world_origin=None)
Parameters:

  • scale – The pixel scale, typically in units of arcsec/pixel.

  • shear – The shear, which should be a Shear instance.

  • origin – Optional origin position for the image coordinate system. If provided, it should be a PositionD or PositionI. [default: PositionD(0., 0.)]

  • world_origin – Optional origin position for the world coordinate system. If provided, it should be a PositionD. [default: PositionD(0., 0.)]

property origin

The image coordinate position to use as the origin.

property scale

The pixel scale.

property shear

The applied Shear.

property world_origin

The world coordinate position to use as the origin.

class galsim.JacobianWCS(dudx, dudy, dvdx, dvdy)[source]

Bases: LocalWCS

This WCS is the most general local linear WCS implementing a 2x2 Jacobian matrix.

The conversion functions are:

u = dudx x + dudy y v = dvdx x + dvdy y

A JacobianWCS has attributes dudx, dudy, dvdx, dvdy that you can access directly if that is convenient. You can also access these as a NumPy array directly with:

>>> J = jac_wcs.getMatrix()

Also, JacobianWCS has another method that other WCS classes do not have. The call:

>>> scale, shear, theta, flip = jac_wcs.getDecomposition()

will return the equivalent expansion, shear, rotation and possible flip corresponding to this transformation. See the docstring for that method for more information.

A JacobianWCS is initialized with the command:

>>> wcs = galsim.JacobianWCS(dudx, dudy, dvdx, dvdy)
Parameters:

  • dudx – du/dx

  • dudy – du/dy

  • dvdx – dv/dx

  • dvdy – dv/dy

property dudx

du/dx

property dudy

du/dy

property dvdx

dv/dx

property dvdy

dv/dy

getDecomposition()[source]

Get the equivalent expansion, shear, rotation and possible flip corresponding to this Jacobian transformation.

A non-singular real matrix can always be decomposed into a symmetric positive definite matrix times an orthogonal matrix:

M = P Q

In our case, P includes an overall scale and a shear, and Q is a rotation and possibly a flip of (x,y) -> (y,x).

( dudx dudy ) = scale/sqrt(1-g1^2-g2^2) ( 1+g1 g2 ) ( cos(theta) -sin(theta) ) F ( dvdx dvdy ) ( g2 1-g1 ) ( sin(theta) cos(theta) )

where F is either the identity matrix, ( 1 0 ), or a flip matrix, ( 0 1 ).

( 0 1 ) ( 1 0 )

If there is no flip, then this means that the effect of:

>>> prof.transform(dudx, dudy, dvdx, dvdy)

is equivalent to:

>>> prof.rotate(theta).shear(shear).expand(scale)

in that order. (Rotation and shear do not commute.)

The decomposition is returned as a tuple: (scale, shear, theta, flip), where scale is a float, shear is a Shear, theta is an Angle, and flip is a bool.

getMatrix()[source]

Get the Jacobian as a NumPy matrix:

numpy.array( [[ dudx, dudy ],

[ dvdx, dvdy ]] )

class galsim.AffineTransform(dudx, dudy, dvdx, dvdy, origin=None, world_origin=None)[source]

Bases: UniformWCS

This WCS is the most general linear transformation. It involves a 2x2 Jacobian matrix and an offset. You can provide the offset in terms of either the image_pos (x0,y0) where (u,v) = (0,0), or the world_pos (u0,v0) where (x,y) = (0,0). Or, in fact, you may provide both, in which case the image_pos (x0,y0) corresponds to the world_pos (u0,v0).

The conversion functions are:

u = dudx (x-x0) + dudy (y-y0) + u0 v = dvdx (x-x0) + dvdy (y-y0) + v0

An AffineTransform has attributes dudx, dudy, dvdx, dvdy, x0, y0, u0, v0 that you can access directly if that is convenient.

An AffineTransform is initialized with the command:

>>> wcs = galsim.AffineTransform(dudx, dudy, dvdx, dvdy, origin=None, world_origin=None)
Parameters:

  • dudx – du/dx

  • dudy – du/dy

  • dvdx – dv/dx

  • dvdy – dv/dy

  • origin – Optional origin position for the image coordinate system. If provided, it should be a PositionD or PositionI. [default: PositionD(0., 0.)]

  • world_origin – Optional origin position for the world coordinate system. If provided, it should be a PositionD. [default: PositionD(0., 0.)]

property dudx

du/dx

property dudy

du/dy

property dvdx

dv/dx

property dvdy

dv/dy

property origin

The image coordinate position to use as the origin.

property world_origin

The world coordinate position to use as the origin.

class galsim.UVFunction(ufunc, vfunc, xfunc=None, yfunc=None, origin=None, world_origin=None, uses_color=False)[source]

Bases: EuclideanWCS

This WCS takes two arbitrary functions for u(x,y) and v(x,y).

The ufunc and vfunc parameters may be:

  • python functions that take (x,y) arguments

  • python objects with a __call__ method that takes (x,y) arguments

  • strings which can be parsed with eval(‘lambda x,y: ‘+str)

You may also provide the inverse functions x(u,v) and y(u,v) as xfunc and yfunc. These are not required, but if you do not provide them, then any operation that requires going from world to image coordinates will raise a NotImplementedError.

Note: some internal calculations will be faster if the functions can take NumPy arrays for x,y and output arrays for u,v. Usually this does not require any change to your function, but it is worth keeping in mind. For example, if you want to do a sqrt, you may be better off using numpy.sqrt rather than math.sqrt.

A UVFunction is initialized with the command:

>>> wcs = galsim.UVFunction(ufunc, vfunc, origin=None, world_origin=None)
Parameters:

  • ufunc – The function u(x,y)

  • vfunc – The function v(x,y)

  • xfunc – The function x(u,v) (optional)

  • yfunc – The function y(u,v) (optional)

  • origin – Optional origin position for the image coordinate system. If provided, it should be a PositionD or PositionI. [default: PositionD(0., 0.)]

  • system. (world_origin Optional origin position for the world coordinate) – If provided, it should be a PositionD. [default: PositionD(0., 0.)]

  • uses_color – If True, then the functions take three parameters (x,y,c) or (u,v,c) where the third term is some kind of color value. (The exact meaning of “color” here is user-defined. You just need to be consistent with the color values you use when using the wcs.) [default: False]

property origin

The image coordinate position to use as the origin.

property ufunc

The input ufunc

property vfunc

The input vfunc

property world_origin

The world coordinate position to use as the origin.

property xfunc

The input xfunc

property yfunc

The input yfunc

Celestial WCS’s

class galsim.RaDecFunction(ra_func, dec_func=None, origin=None)[source]

Bases: CelestialWCS

This WCS takes an arbitrary function for the Right Ascension (ra) and Declination (dec).

In many cases, it can be more convenient to calculate both ra and dec in a single function, since there will typically be intermediate values that are common to both, so it may be more efficient to just calculate those once and thence calculate both ra and dec. Thus, we provide the option to provide either a single function or two separate functions.

The function parameters used to initialize an RaDecFunction may be:

  • a python functions that take (x,y) arguments

  • a python object with a __call__ method that takes (x,y) arguments

  • a string which can be parsed with eval(‘lambda x,y: ‘+str)

The return values, ra and dec, should be given in _radians_.

The first argument is called ra_func, but if dec_func is omitted, then it is assumed to calculate both ra and dec. The two values should be returned as a tuple (ra,dec).

We don’t want a function that returns Angle instances, because we want to allow for the possibility of using NumPy arrays as inputs and outputs to speed up some calculations. The function isn’t _required_ to work with NumPy arrays, but it is possible that some things will be faster if it does. If it were expected to return Angle instances, then it definitely couldn’t work with arrays.

An RaDecFunction is initialized with either of the following commands:

>>> wcs = galsim.RaDecFunction(radec_func, origin=None)
>>> wcs = galsim.RaDecFunction(ra_func, dec_func, origin=None)
Parameters:

  • ra_func – If dec_func is also given: A function ra(x,y) returning ra in radians. If dec_func=None: A function returning a tuple (ra,dec), both in radians.

  • dec_func – Either a function dec(x,y) returning dec in radians, or None (in which case ra_func is expected to return both ra and dec. [default: None]

  • origin – Optional origin position for the image coordinate system. If provided, it should be a PositionD or PositionI. [default: PositionD(0., 0.)]

property origin

The image coordinate position to use as the origin.

property radec_func

The input radec_func

class galsim.AstropyWCS(file_name=None, dir=None, hdu=None, header=None, compression='auto', wcs=None, origin=None)[source]

Bases: CelestialWCS

This WCS uses astropy.wcs to read WCS information from a FITS file. It requires the astropy.wcs python module to be installed.

Astropy may be installed using pip, fink, or port:

>>> pip install astropy
>>> fink install astropy-py27
>>> port install py27-astropy

It also comes by default with Enthought and Anaconda. For more information, see their website:

An AstropyWCS is initialized with one of the following commands:

>>> wcs = galsim.AstropyWCS(file_name=file_name)  # Open a file on disk
>>> wcs = galsim.AstropyWCS(header=header)        # Use an existing pyfits header
>>> wcs = galsim.AstropyWCS(wcs=wcs)              # Use an existing astropy.wcs.WCS instance

Exactly one of the parameters file_name, header or wcs is required. Also, since the most common usage will probably be the first, you can also give a file_name without it being named:

>>> wcs = galsim.AstropyWCS(file_name)
Parameters:

  • file_name – The FITS file from which to read the WCS information. This is probably the usual parameter to provide. [default: None]

  • dir – Optional directory to prepend to file_name. [default: None]

  • hdu – Optionally, the number of the HDU to use if reading from a file. The default is to use either the primary or first extension as appropriate for the given compression. (e.g. for rice, the first extension is the one you normally want.) [default: None]

  • header – The header of an open pyfits (or astropy.io) hdu. Or, it can be a FitsHeader object. [default: None]

  • compression – Which decompression scheme to use (if any). See galsim.fits.read() for the available options. [default: ‘auto’]

  • wcs – An existing astropy.wcs.WCS instance [default: None]

  • origin – Optional origin position for the image coordinate system. If provided, it should be a PositionD or PositionI. [default: None]

property origin

The origin in image coordinates of the WCS function.

property wcs

The underlying astropy.wcs.WCS object.

class galsim.PyAstWCS(file_name=None, dir=None, hdu=None, header=None, compression='auto', wcsinfo=None, origin=None)[source]

Bases: CelestialWCS

This WCS uses PyAst (the python front end for the Starlink AST code) to read WCS information from a FITS file. It requires the starlink.Ast python module to be installed.

Starlink may be installed using pip:

>>> pip install starlink-pyast

For more information, see their website:

https://pypi.python.org/pypi/starlink-pyast/

A PyAstWCS is initialized with one of the following commands:

>>> wcs = galsim.PyAstWCS(file_name=file_name)  # Open a file on disk
>>> wcs = galsim.PyAstWCS(header=header)        # Use an existing pyfits header
>>> wcs = galsim.PyAstWCS(wcsinfo=wcsinfo)      # Use an existing starlink.Ast.FrameSet

Exactly one of the parameters file_name, header or wcsinfo is required. Also, since the most common usage will probably be the first, you can also give a file name without it being named:

>>> wcs = galsim.PyAstWCS(file_name)
Parameters:

  • file_name – The FITS file from which to read the WCS information. This is probably the usual parameter to provide. [default: None]

  • dir – Optional directory to prepend to file_name. [default: None]

  • hdu – Optionally, the number of the HDU to use if reading from a file. The default is to use either the primary or first extension as appropriate for the given compression. (e.g. for rice, the first extension is the one you normally want.) [default: None]

  • header – The header of an open pyfits (or astropy.io) hdu. Or, it can be a FitsHeader object. [default: None]

  • compression – Which decompression scheme to use (if any). See galsim.fits.read() for the available options. [default:’auto’]

  • wcsinfo – An existing starlink.Ast.FrameSet [default: None]

  • origin – Optional origin position for the image coordinate system. If provided, it should be a PositionD or PositionI. [default: None]

property origin

The origin in image coordinates of the WCS function.

property wcsinfo

The underlying starlink.Ast.FrameSet for this object.

class galsim.WcsToolsWCS(file_name, dir=None, origin=None)[source]

Bases: CelestialWCS

This WCS uses wcstools executables to perform the appropriate WCS transformations for a given FITS file. It requires wcstools command line functions to be installed.

Note: It uses the wcstools executables xy2sky and sky2xy, so it can be quite a bit less efficient than other options that keep the WCS in memory.

See their website for information on downloading and installing wcstools:

A WcsToolsWCS is initialized with the following command:

>>> wcs = galsim.WcsToolsWCS(file_name)
Parameters:

  • file_name – The FITS file from which to read the WCS information.

  • dir – Optional directory to prepend to file_name. [default: None]

  • origin – Optional origin position for the image coordinate system. If provided, it should be a PositionD or PositionI. [default: None]

property file_name

The file name of the FITS file with the WCS information.

property origin

The origin in image coordinates of the WCS function.

class galsim.GSFitsWCS(file_name=None, dir=None, hdu=None, header=None, compression='auto', origin=None, _data=None, _doiter=True)[source]

Bases: CelestialWCS

This WCS uses a GalSim implementation to read a WCS from a FITS file.

It doesn’t do nearly as many WCS types as the other options, and it does not try to be as rigorous about supporting all possible valid variations in the FITS parameters. However, it does several popular WCS types properly, and it doesn’t require any additional python modules to be installed, which can be helpful.

Currrently, it is able to parse the following WCS types: TAN, STG, ZEA, ARC, TPV, TNX

A GSFitsWCS is initialized with one of the following commands:

>>> wcs = galsim.GSFitsWCS(file_name=file_name)  # Open a file on disk
>>> wcs = galsim.GSFitsWCS(header=header)        # Use an existing pyfits header

Also, since the most common usage will probably be the first, you can also give a file name without it being named:

>>> wcs = galsim.GSFitsWCS(file_name)

In addition to reading from a FITS file, there is also a factory function that builds a GSFitsWCS object implementing a TAN projection. See the docstring of TanWCS for more details.

Parameters:

  • file_name – The FITS file from which to read the WCS information. This is probably the usual parameter to provide. [default: None]

  • dir – Optional directory to prepend to file_name. [default: None]

  • hdu – Optionally, the number of the HDU to use if reading from a file. The default is to use either the primary or first extension as appropriate for the given compression. (e.g. for rice, the first extension is the one you normally want.) [default: None]

  • header – The header of an open pyfits (or astropy.io) hdu. Or, it can be a FitsHeader object. [default: None]

  • compression – Which decompression scheme to use (if any). See galsim.fits.read() for the available options. [default: ‘auto’]

  • origin – Optional origin position for the image coordinate system. If provided, it should be a PositionD or PositionI. [default: None]

property origin

The origin in image coordinates of the WCS function.

galsim.FitsWCS(file_name=None, dir=None, hdu=None, header=None, compression='auto', text_file=False, suppress_warning=False)[source]

This factory function will try to read the WCS from a FITS file and return a WCS that will work. It tries a number of different WCS classes until it finds one that succeeds in reading the file.

If none of them work, then the last class it tries, AffineTransform, is guaranteed to succeed, but it will only model the linear portion of the WCS (the CD matrix, CRPIX, and CRVAL), using reasonable defaults if even these are missing. If you think that you have the right software for one of the WCS types, but FitsWCS still defaults to AffineTransform, it may be helpful to update your installation of astropy and/or starlink (if you don’t already have the latest version).

Note: The list of classes this function will try may be edited, e.g. by an external module that wants to add an additional WCS type. The list is galsim.fitswcs.fits_wcs_types.

Parameters:

  • file_name – The FITS file from which to read the WCS information. This is probably the usual parameter to provide. [default: None]

  • dir – Optional directory to prepend to file_name. [default: None]

  • hdu – Optionally, the number of the HDU to use if reading from a file. The default is to use either the primary or first extension as appropriate for the given compression. (e.g. for rice, the first extension is the one you normally want.) [default: None]

  • header – The header of an open pyfits (or astropy.io) hdu. Or, it can be a FitsHeader object. [default: None]

  • compression – Which decompression scheme to use (if any). See galsim.fits.read() for the available options. [default: ‘auto’]

  • text_file – Normally a file is taken to be a fits file, but you can also give it a text file with the header information (like the .head file output from SCamp). In this case you should set text_file = True to tell GalSim to parse the file this way. [default: False]

  • suppress_warning – Whether to suppress a warning that the WCS could not be read from the FITS header, so the WCS defaulted to either a PixelScale or AffineTransform. [default: False] (Note: this is (by default) set to True when this function is implicitly called from one of the galsim.fits.read* functions.)

galsim.TanWCS(affine, world_origin, units=coord.arcsec)[source]

This is a function that returns a GSFitsWCS object for a TAN WCS projection.

The TAN projection is essentially an affine transformation from image coordinates to Euclidean (u,v) coordinates on a tangent plane, and then a “deprojection” of this plane onto the sphere given a particular RA, Dec for the location of the tangent point. The tangent point will correspond to the location of (u,v) = (0,0) in the intermediate coordinate system.

Parameters:

  • affine – An AffineTransform defining the transformation from image coordinates to the coordinates on the tangent plane.

  • world_origin – A CelestialCoord defining the location on the sphere where the tangent plane is centered.

  • units – The angular units of the (u,v) intermediate coordinate system. [default: galsim.arcsec]

Returns:

a GSFitsWCS describing this WCS.

galsim.FittedSIPWCS(x, y, ra, dec, wcs_type='TAN', order=3, center=None)[source]

A WCS constructed from a list of reference celestial and image coordinates.

Parameters:

  • x – Image x-coordinates of reference stars in pixels

  • y – Image y-coordinates of reference stars in pixels

  • ra – Right ascension of reference stars in radians

  • dec – Declination of reference stars in radians

  • wcs_type – The type of the tangent plane projection to use. Should be one of [‘TAN’, ‘STG’, ‘ZEA’, or ‘ARC’]. [default: ‘TAN’]

  • order – The order of the Simple Imaging Polynomial (SIP) used to describe the WCS distortion. SIP coefficients kick in when order >= 2. If you supply order=1, then just fit a WCS without any SIP coefficients. [default: 3]

  • center – A CelestialCoord defining the location on the sphere where the tangent plane is centered. [default: None, which means use the average position of the list of reference stars]

Celestial Coordinates

Our CelestialCoord class is currently hosted as part of the LSSTDESC.Coord package:

https://github.com/LSSTDESC/Coord

An earlier version of this code was originally implemented in GalSim, so we still import the relevant classes into the galsim namespace. You may therefore use either galsim.CelestialCoord or coord.CelestialCoord as they are equivalent.

class galsim.CelestialCoord(ra, dec=None)[source]

This class defines a position on the celestial sphere, normally given by two angles, ra and dec.

This class can be used to perform various calculations in spherical coordinates, such as finding the angular distance between two points in the sky, calculating the angles in spherical triangles, projecting from sky coordinates onto a Euclidean tangent plane, etc.

Initialization:

A CelestialCoord object is constructed from the right ascension and declination:

coord.CelestialCoord.__init__()

>>> c = CelestialCoord(ra=12*hours, dec=31*degrees)
>>> print(c)
coord.CelestialCoord(3.141592653589793 radians, 0.5410520681182421 radians)

Attributes:

A CelestialCoord has the following (read-only) attributes:

ra:

The right ascension (an Angle instance)

dec:

The declination (an Angle instance)

>>> print(c.ra / degrees, c.dec / degrees)
180.0 31.0

In addition there is a convenience access property that returns ra and dec in radians.

rad:

A tuple (ra.rad, dec.rad)

>>> print(c.rad)
(3.141592653589793, 0.5410520681182421)

Spherical Geometry:

The basic spherical geometry operations are available to work with spherical triangles

For three coordinates cA, cB, cC making a spherical triangle, one can calculate the sides and angles via:

coord.CelestialCoord.distanceTo()
coord.CelestialCoord.angleBetween()
>>> cA = CelestialCoord(0 * degrees, 0 * degrees)
>>> cB = CelestialCoord(0 * degrees, 10 * degrees)
>>> cC = CelestialCoord(10 * degrees, 0 * degrees)
>>> a = cB.distanceTo(cC)
>>> b = cC.distanceTo(cA)
>>> c = cA.distanceTo(cB)
>>> print(a.deg, b.deg, c.deg)
14.106044260566366 10.0 10.0
>>> A = cA.angleBetween(cB, cC)
>>> B = cB.angleBetween(cC, cA)
>>> C = cC.angleBetween(cA, cB)
>>> print(A.deg, B.deg, C.deg)
90.0 45.43854858674231 45.43854858674231

Projections:

Local tangent plane projections of an area of the sky can be performed using the project method:

coord.CelestialCoord.project()

>>> center = CelestialCoord(ra=10*hours, dec=30*degrees)
>>> sky_coord = CelestialCoord(ra=10.5*hours, dec=31*degrees)
>>> print(sky_coord)
coord.CelestialCoord(2.748893571891069 radians, 0.5410520681182421 radians)
>>> u, v = center.project(sky_coord)
>>> print(u.deg, v.deg)
-6.452371275343261 1.21794987288635

and back:

coord.CelestialCoord.deproject()

>>> sky_coord = center.deproject(u,v)
>>> print(sky_coord)
coord.CelestialCoord(2.748893571891069 radians, 0.5410520681182421 radians)

where u and v are Angles and center and sky_coord are CelestialCoords.

angleBetween(coord2, coord3)[source]

Find the open angle at the location of the current coord between coord2 and coord3.

The current coordinate along with the two other coordinates form a spherical triangle on the sky. This function calculates the angle between the two sides at the location of the current coordinate.

Note that this returns a signed angle. The angle is positive if the sweep direction from coord2 to coord3 is counter-clockwise (as observed from Earth). It is negative if the direction is clockwise.

Parameters:

  • coord2 – A second CelestialCoord

  • coord3 – A third CelestialCoord

Returns:

the angle between the great circles joining the other two coordinates to the current coordinate.

area(coord2, coord3)[source]

Find the area of a spherical triangle in steradians.

The current coordinate along with the two other coordinates form a spherical triangle on the sky. This function calculates the area of that spherical triangle, which is measured in steradians (i.e. surface area of the triangle on the unit sphere).

Parameters:

  • coord2 – A second CelestialCoord

  • coord3 – A third CelestialCoord

Returns:

the area in steradians of the given spherical triangle.

property dec

A read-only attribute, giving the Declination as an Angle

deproject(u, v, projection=None)[source]

Do the reverse process from the project() function.

i.e. This takes in a position (u,v) and returns the corresponding celestial coordinate, using the current coordinate as the center point of the tangent plane projection.

Parameters:

  • u – The u position on the tangent plane to deproject (must be an Angle instance)

  • v – The v position on the tangent plane to deproject (must be an Angle instance)

  • projection – The name of the projection to be used. [default: gnomonic, see project docstring for other options]

Returns:

the corresponding CelestialCoord for that position.

deproject_rad(u, v, projection=None)[source]

This is basically identical to the deproject() function except that the output ra, dec are returned as a tuple (ra, dec) in radians rather than packaged as a CelestialCoord object and u and v are in radians rather than Angle instances.

The main advantage to this is that it will work if u and v are NumPy arrays, in which case the output ra, dec will also be NumPy arrays.

Parameters:

  • u – The u position in radians on the tangent plane to deproject

  • v – The v position in radians on the tangent plane to deproject

  • projection – The name of the projection to be used. [default: gnomonic, see project docstring for other options]

Returns:

the corresponding RA, Dec in radians

distanceTo(coord2)[source]

Returns the great circle distance between this coord and another one. The return value is an Angle object

Parameters:

coord2 – The CelestialCoord to calculate the distance to.

Returns:

the great circle distance to coord2.

ecliptic(epoch=2000.0, date=None)[source]

Get the longitude and latitude in ecliptic coordinates corresponding to this position.

The epoch parameter is used to get an accurate value for the (time-varying) obliquity of the ecliptic. The formulae for this are quite straightforward. It requires just a single parameter for the transformation, the obliquity of the ecliptic (the Earth’s axial tilt).

Parameters:

  • epoch – The epoch to be used for estimating the obliquity of the ecliptic, if date is None. But if date is given, then use that to determine the epoch. [default: 2000.]

  • date – If a date is given as a python datetime object, then return the position in ecliptic coordinates with respect to the sun position at that date. If None, then return the true ecliptic coordiantes. [default: None]

Returns:

the longitude and latitude as a tuple (lambda, beta), given as Angle instances.

static from_ecliptic(lam, beta, epoch=2000.0, date=None)[source]

Create a CelestialCoord from the given ecliptic coordinates

Parameters:

  • lam – The longitude in ecliptic coordinates (an Angle instance)

  • beta – The latitude in ecliptic coordinates (an Angle instance)

  • epoch – The epoch to be used for estimating the obliquity of the ecliptic, if date is None. But if date is given, then use that to determine the epoch. [default: 2000.]

  • date – If a date is given as a python datetime object, then return the position in ecliptic coordinates with respect to the sun position at that date. If None, then return the true ecliptic coordiantes. [default: None]

Returns:

the CelestialCoord corresponding to these ecliptic coordinates.

static from_galactic(el, b, epoch=2000.0)[source]

Create a CelestialCoord from the given galactic coordinates

Parameters:

  • el – The longitude in galactic coordinates (an Angle instance)

  • b – The latitude in galactic coordinates (an Angle instance)

  • epoch – The epoch of the returned coordinate. [default: 2000.]

Returns:

the CelestialCoord corresponding to these galactic coordinates.

static from_xyz(x, y, z)[source]

Construct a CelestialCoord from a given (x,y,z) position in three dimensions.

The 3D (x,y,z) position does not need to fall on the unit sphere. The RA, Dec will be inferred from the relations:

\[\begin{split}x &= r \cos(dec) \cos(ra) \\ y &= r \cos(dec) \sin(ra) \\ z &= r \sin(dec)\end{split}\]

where \(r\) is arbitrary.

Parameters:

  • x – The x position in 3 dimensions. Corresponds to r cos(dec) cos(ra)

  • y – The y position in 3 dimensions. Corresponds to r cos(dec) sin(ra)

  • z – The z position in 3 dimensions. Corresponds to r sin(dec)

Returns:

a CelestialCoord instance

galactic(epoch=2000.0)[source]

Get the longitude and latitude in galactic coordinates corresponding to this position.

Parameters:

epoch – The epoch of the current coordinate. [default: 2000.]

Returns:

the longitude and latitude as a tuple (el, b), given as Angle instances.

get_xyz()[source]

Get the (x,y,z) coordinates on the unit sphere corresponding to this (RA, Dec).

\[\begin{split}x &= \cos(dec) \cos(ra) \\ y &= \cos(dec) \sin(ra) \\ z &= \sin(dec)\end{split}\]
Returns:

a tuple (x,y,z)

greatCirclePoint(coord2, theta)[source]

Returns a point on the great circle connecting self and coord2.

Two points, c1 and c2, on the unit sphere define a great circle (so long as the two points are not either coincident or antipodal). We can define points on this great circle by their angle from c1, such that the angle for c2 has 0 < theta2 < pi. I.e. theta increases from 0 as the points move from c1 towards c2.

This function then returns the coordinate on this great circle (where c1 is self and c2 is coord2) that corresponds to the given angle theta.

Parameters:

  • coord2 – Another CelestialCoord defining the great circle to use.

  • theta – The Angle along the great circle corresponding to the desired point.

Returns:

the corresponding CelestialCoord

jac_deproject(u, v, projection=None)[source]

Return the jacobian of the deprojection.

i.e. if the input position is (u,v) then the return matrix is

\[\begin{split}J \equiv \begin{bmatrix} J00 & J01 \\ J10 & J11 \end{bmatrix} = \begin{bmatrix} d\textrm{ra}/du \cos(\textrm{dec}) & d\textrm{ra}/dv \cos(\textrm{dec}) \\ d\textrm{dec}/du & d\textrm{dec}/dv \end{bmatrix}\end{split}\]
Parameters:

  • u – The u position (as an Angle instance) on the tangent plane

  • v – The v position (as an Angle instance) on the tangent plane

  • projection – The name of the projection to be used. [default: gnomonic, see project docstring for other options]

Returns:

the Jacobian as a 2x2 numpy array [[J00, J01], [J10, J11]]

jac_deproject_rad(u, v, projection=None)[source]

Equivalent to jac_deproject, but the inputs are in radians and may be numpy arrays.

Parameters:

  • u – The u position (in radians) on the tangent plane

  • v – The v position (in radians) on the tangent plane

  • projection – The name of the projection to be used. [default: gnomonic, see project docstring for other options]

Returns:

the Jacobian as a 2x2 numpy array [[J00, J01], [J10, J11]]

normal()[source]

Return the coordinate in the “normal” convention of having 0 <= ra < 24 hours.

This convention is not enforced on construction, so this function exists to make it easy to convert if desired.

Functions such as from_galactic and from_xyz will return normal coordinates.

precess(from_epoch, to_epoch)[source]
This function precesses equatorial ra and dec from one epoch to another.

It is adapted from a set of fortran subroutines based on (a) pages 30-34 of the Explanatory Supplement to the AE, (b) Lieske, et al. (1977) A&A 58, 1-16, and (c) Lieske (1979) A&A 73, 282-284.

Parameters:

  • from_epoch – The epoch of the current coordinate

  • to_epoch – The epoch of the returned precessed coordinate

Returns:

a CelestialCoord object corresponding to the precessed position.

project(coord2, projection=None)[source]

Use the currect coord as the center point of a tangent plane projection to project the coord2 coordinate onto that plane.

This function return a tuple (u,v) in the Euclidean coordinate system defined by a tangent plane projection around the current coordinate, with +v pointing north and +u pointing west. (i.e. to the right on the sky if +v is up.)

There are currently four options for the projection, which you can specify with the optional projection keyword argument:

gnomonic:

[default] uses a gnomonic projection (i.e. a projection from the center of the sphere, which has the property that all great circles become straight lines. For more information, see http://mathworld.wolfram.com/GnomonicProjection.html This is the usual TAN projection used by most FITS images.

stereographic:

uses a stereographic proejection, which preserves angles, but not area. For more information, see http://mathworld.wolfram.com/StereographicProjection.html

lambert:

uses a Lambert azimuthal projection, which preserves area, but not angles. For more information, see http://mathworld.wolfram.com/LambertAzimuthalEqual-AreaProjection.html

postel:

uses a Postel equidistant proejection, which preserves distances from the projection point, but not area or angles. For more information, see http://mathworld.wolfram.com/AzimuthalEquidistantProjection.html

The distance or angle errors increase with distance from the projection point of course.

Parameters:

  • coord2 – The coordinate to project onto the tangent plane.

  • projection – The name of the projection to be used. [default: gnomonic, see above for other options]

Returns:

(u,v) as Angle instances

project_rad(ra, dec, projection=None)[source]

This is basically identical to the project() function except that the input ra, dec are given in radians rather than packaged as a CelestialCoord object and the returned u,v are given in radians.

The main advantage to this is that it will work if ra and dec are NumPy arrays, in which case the output u, v will also be NumPy arrays.

Parameters:

  • ra – The right ascension in radians to project onto the tangent plane.

  • dec – The declination in radians to project onto the tangent plane.

  • projection – The name of the projection to be used. [default: gnomonic, see project docstring for other options]

Returns:

(u,v) in radians

property ra

A read-only attribute, giving the Right Ascension as an Angle

property rad

A convenience property, giving a tuple (ra.rad, dec.rad)

static radec_to_xyz(ra, dec, r=1.0)[source]

Convert ra, dec (in radians) to 3D x,y,z coordinates on the unit sphere.

The connection between (ra,dec) and (x,y,z) are given by the following formulae: .. math:

x &= r \cos(dec) \cos(ra)  \\
y &= r \cos(dec) \sin(ra)  \\
z &= r \sin(dec)

For a single ra,dec pair, the following are essentially equivalent:

>>> ra = 12*hours/radians       # May be any angle measured
>>> dec = 31*degrees/radians    # in radians

>>> CelestialCoord.radec_to_xyz(ra, dec)
(-0.8571673007021123, 1.0497271911386187e-16, 0.5150380749100542)

>>> CelestialCoord(ra * radians, dec * radians).get_xyz()
(-0.8571673007021123, 1.0497271911386187e-16, 0.5150380749100542)

However, the advantage of this function is that the input values may be numpy arrays, in which case, the return values will also be numpy arrays.

Parameters:

  • ra – The right ascension(s) in radians. May be a numpy array.

  • dec – The declination(s) in radians. May be a numpy array.

  • r – The distance(s) from Earth (default 1.). May be a numpy array.

Returns:

x, y, z as a tuple.

static xyz_to_radec(x, y, z, return_r=False)[source]

Convert 3D x,y,z coordinates to ra, dec (in radians).

The connection between (ra,dec) and (x,y,z) are given by the following formulae: .. math:

x &= r \cos(dec) \cos(ra)  \\
y &= r \cos(dec) \sin(ra)  \\
z &= r \sin(dec)

For a single (x,y,z) position, the following are essentially equivalent:

>>> x = 0.839       # May be any any 3D location
>>> y = 0.123       # Not necessarily on unit sphere
>>> z = 0.530

>>> CelestialCoord.xyz_to_radec(x, y, z)
(0.14556615088111796, 0.558616191048523)

>>> c = CelestialCoord.from_xyz(x, y, z)
>>> c.ra.rad, c.dec.rad
(0.145566150881118, 0.558616191048523)

However, the advantage of this function is that the input values may be numpy arrays, in which case, the return values will also be numpy arrays.

Parameters:

  • x – The x position(s) in 3 dimensions. May be a numpy array.

  • y – The y position(s) in 3 dimensions. May be a numpy array.

  • z – The z position(s) in 3 dimensions. May be a numpy array.

  • return_r – Whether to return r as well as ra, dec. (default: False)

Returns:

ra, dec as a tuple. Or if return_r is True, (ra, dec, r).

WCS Utilities

galsim.wcs.readFromFitsHeader(header, suppress_warning=True)[source]

Read a WCS function from a FITS header.

This is normally called automatically from within the galsim.fits.read function, but you can also call it directly as:

wcs, origin = galsim.wcs.readFromFitsHeader(header)

If the file was originally written by GalSim using one of the galsim.fits.write() functions, then this should always succeed in reading back in the original WCS. It may not end up as exactly the same class as the original, but the underlying world coordinate system transformation should be preserved.

Note

For UVFunction and RaDecFunction, if the functions that were written to the FITS header were real python functions (rather than a string that is converted to a function), then the mechanism we use to write to the header and read it back in has some limitations:

  1. It apparently only works for cpython implementations.

  2. It probably won’t work to write from one version of python and read from another. (At least for major version differences.)

  3. If the function uses globals, you’ll need to make sure the globals are present when you read it back in as well, or it probably won’t work.

  4. It looks really ugly in the header.

  5. We haven’t thought much about the security implications of this, so beware using GalSim to open FITS files from untrusted sources.

If the file was not written by GalSim, then this code will do its best to read the WCS information in the FITS header. Depending on what kind of WCS is encoded in the header, this may or may not be successful.

If there is no WCS information in the header, then this will default to a pixel scale of 1.

In addition to the wcs, this function will also return the image origin that the WCS is assuming for the image. If the file was originally written by GalSim, this should correspond to the original image origin. If not, it will default to (1,1).

Parameters:

  • header – The fits header with the WCS information.

  • suppress_warning – Whether to suppress a warning that the WCS could not be read from the FITS header, so the WCS defaulted to either a PixelScale or AffineTransform. [default: True]

Returns:

a tuple (wcs, origin) of the wcs from the header and the image origin.

galsim.wcs.compatible(wcs1, wcs2)[source]

A utility to check the compatibility of two WCS. In particular, if two WCS are consistent with each other modulo a shifted origin, we consider them to be compatible, even though they are not equal.