GalSim-developers/GalSim — pull request #500 by rmjarvis (2014-01-14 16:30:42+00:00)
Title: #364 New WCS classes.

This is a really huge pull request, so I apologize in advance. It adds in a lot of new functionality with respect to the WCS classes. Here is the summary from the CHANGELOG:

New WCS classes:

* Every place in the code that used to need a pixel scale item (e.g. `Image` constructor, `GSObject.draw()`, `InterpolatedImage`, etc.) now can take a `wcs` item. The `scale` parameter is still an option, but now it is just
shorthand for `wcs = PixelScale(scale)`.
* There are three classes that we call "local WCS classes":
* `PixelScale` describes a simple scale conversion from pixels to arcsec.
* `ShearWCS` describes a uniformly sheared coordinate system.
* `JacobianWCS` describes an arbitrary 2x2 Jacobian matrix.
* There are four non-local WCS classes that use Euclidean coordinates for the world coordinate system:
* `OffsetWCS` is a `PixelScale` with the world (0,0) location offset from the (0,0) position in image coordinates.
* `OffsetShearWCS` is a `ShearWCS` with a similar offset.
* `AffineTransform` is a `JacobianWCS` with an offset. It is the most general possible _uniform_ WCS transformation. i.e. one where the pixel shape is uniform across the image.
* `UVFunction` is an arbitrary transformation from (x,y) coordinates to Euclidean (u,v) coordinates. It takes arbitrary function u(x,y) and v(x,y) as inputs. (And optionally x(u,v) and y(u,v) for the inverse transformations.)
* There are five WCS classes that use celestial coordinates for the world coordinate system. i.e. the world coordinates are in terms of right ascension and declination (RA, Dec). There is a new CelestialCoord class that encapsulates this kind of position on the sphere.
* `RaDecFunction` takes an arbitrary function radec_func(x,y) that returns the RA and Dec.
* `AstropyWCS` uses the astropy.wcs package to read in a given FITS file.
* `PyAstWCS` uses the starlink.Ast package to read in a given FITS file.
* `WcsToolsWCS` uses wcstools commands for a given FITS file.
* `GSFitsWCS` is GalSim code to read FITS files that use TAN and TPV WCS types. Less flexible than the others, but still useful since these are probably the most common WCS types for optical astronomical images. Plus it is quite a bit faster than the others.
* Finally, there is a factory function called `FitsWCS` that will try the various classes that can read FITS files until it finds one that works. It will revert to `AffineTransform` if it cannot find anything better.
* Another function, `TanWCS`, acts like a WCS class. It builds a WCS using TAN projection and returns a `GSFitsWCS` implementing it.
* When reading in an image from a FITS file, the image will automatically try to read the WCS information from the header with the `FitsWCS` function.
* See the docstring for `BaseWCS` (the base class for all of these WCS classes) for information about how to use these classes.
* Also, check out demo3, demo9, demo10, and demo11 for example usage.

To facilitate the celestial classes, there is a new `CelestialCoord` class to define a position on the sky. It does a number of useful things that I use in the code, plus has a few other spherical trig things that I didn't need here, but that I had coded up already in a C++ class I had written a while back for other purposes. So I put them in as well for completeness.

Here is the bit from CHANGELOG:

New `CelestialCoord` class:

* This class describes a position on the celestial sphere according to RightAscension (RA) and Declination (Dec). These two values are accessible as coord.ra and coord.dec. So it is used by some of the WCS classes for
the world coordinate positions.
* It has methods to handle a number of spherical trigonometry operations that are sometimes required when dealing with celestial coordinates:
* `coord1.diatanceTo(coord2)` returns the great circle distance between two coordinates (as a `galsim.Angle`).
* `coord1.angleBetween(coord2,coord3)` returns the angle between the two great circles (coord1-coord2) and (coord1-coord3).
* `coord1.project(coord2)` applies a tangent plane projection of coord2 with respect to the tangent point coord1 using one of 4 possible projection schemes specified by the optional keyword `projection`: lambert, stereographic, gnomonic, or postel. See the docstring for this function for details.
* `coord1.deproject(pos)` reverses the projection to go from the position on the tangent plane back to celestial coordinates.
* `coord.precess(from_epoch, to_epoch)` precesses the coordinates to a different epoch.
* `coord.getGalaxyPos()` returns the galaxy longitude and latitude as a tuple (el, b).

Other updates (from CHANGELOG again):

Python layer API changes:

* `Angle.wrap()` now returns the wrapped angle rather than modifying the original.
* `angle.wrap()` should now be `angle = angle.wrap()`.
* Removed the previously deprecated `Ellipse` and `AtmosphericPSF` classes. Also removed `PhotonArray` from the python layer, since it is only used by the C++ layer.
* Changed Bounds methods `addBorder`, `shift`, and `expand` to return new Bounds objects rather than changing the original (in the python layer only).

Updates to config options:

* Changed the previous behavior of the `image.wcs` field to allow several WCS types: `PixelScale`, `Offset`, `Shear`, `OffsetShear`, `UVFunction`, `RaDecFunction`, `Fits`, and `Tan`.
* Changed the name of `sky_pos` to `world_pos`.
* Remove pix top layer in config structure. Add `draw_method=no_pixel` to do what `pix : None` used to do.
* Add `draw_method=real_space` to try to use real-space convolution. This had been an option for the psf draw, but not the main draw. This is only possible if there is only one items being convolved with the pixel.
* Add ability to index Sequences by any running index, rather than just the default. i.e. `obj_num`, `image_num`, or `file_num`.
* Removed `sky_level_pixel` option for specifying the sky. Now the sky level can only be specified in units of ADU/arcsec^2 (or, technically, whatever the world units are if not arcsec) via the `sky_level` value. The alternate `sky_level_pixel` option was never used in our demos, so I doubt many people have been using it.

Other new features:

* New functions to convert an angle to/from DMS strings. Sometimes handy when dealing with RA or Dec.
* `angle.dms()` returns the angle as a string in the form +/-ddmmss.decimal.
* `angle.hms()` returns the angle as a string in the form +/-hhmmss.decimal.
* `angle = DMS_Angle(str)` convert from a dms string to a `galsim.Angle`.
* `angle = HMS_Angle(str)` convert from an hms string to a `galsim.Angle`.
* `profile.applyTransformation(dudx, dudy, dvdx, dvdy)` applies a general (linear) coordinate transformation to a GSObject profile. It is a generalization of `applyShear`, `applyRotation`, etc. There is also the corresponding `createTransformed` as well.
* A new `galsim.fits.readFile()` function reads a FITS file and returns the hdu_list. Normally, this is equivalent to `pyfits.open(file_name)`, but it has a `compression` option that works the same way `compression` works for the other `galsim.fits.read*` functions, so it may be convenient at times.