Add ability to handle WCS's more complicated than just a pixel_scale

[rmjarvis]

Currently, we use pixel_scale for our mapping from pixel coordinates to physical coordinates (typically arcsec). However, it would be nice to be able to use more complicated WCS's. At least a CD matrix that is constant across the field would be a minimum addition I think.

I'll probably also add a version where u(x,y), v(x,y) are arbitrary functions, since that's the most general case, so if that works, then it should be fairly straightforward to add in other specific cases as needed.

Finally, one specific case that I'll add is the DES WCS, read in from the image header. As mentioned in #363, this is required for the DES_Shapelet class to produce accurate images, since it works in sky coordinates, which means using the real WCS. So I'll add that to the des module as part of this issue too.

[rmjarvis]

Some thoughts about making Galsim WCS-enabled.

First, Galsim already has a kind of WCS in there, which we call the pixel_scale. Every time we do

x = ix * pixel_scale
y = iy * pixel_scale

(or sometimes dx rather than pixel_scale), what we really mean is:

x,y = wcs.transform(ix,iy)

where the wcs in this case is simply scaling by the pixel scale.

So I think what we ought to do is the following:

• BaseWCS should be a base class that just defines the syntax using mostly pure virtual functions.
• PixelScale will be a derived class that implement just a pixel_scale factor.
• AffineTransform will be a 2x2 transformation with a shift.
• sky_pos = wcs.transform(chip_pos) should take either PositionI or PositionD for chip_pos and return a PositionD.
• chip_pos = wcs.inverseTransform(sky_pos) should return a PositionD.
• wcs.localAffineTransform(chip_pos) should return an AffineTranform.
• The Image class should carry around a wcs item rather than a scale.
• The Image constructor should accept either a pixel_scale (since that will still be the most common use case) or a wcs. For the former case, it would trivially convert to a PixelScale.
• Everywhere we have a pixel_scale or dx parameter, we should allow either a pixel_scale or a wcs object. (I'd like to move away from using dx, since it is not so self-explanatory as to its meaning.
• The draw command should default to using the local affine transform at the center of the object, rather than call wcs.transform() for every position. However, we can have an option, say use_exact_wcs, that would explicitly call it for every location if the WCS changes significantly across the size of a galaxy. (I'm think of the DES glowing edges here!)

While I'm at it, I think I'd like to get rid of the Ellipse class. It's essentially an AffineTransform without the rotation component, so I think any places we currently use Ellipse can be replaced by AffineTransform.

Comments?

[reikonakajima]

So I guess any shear/shift/dilation operation would be defined as an AffineTransform? (Reflection and rotations are affine transforms as well.) The Ellipse class can be useful for defining galaxy shapes, though? It makes a handy catalog entry for individual galaxies; should we really get rid of the Ellipse class?

[rmjarvis]

The Ellipse class can be useful for defining galaxy shapes, though? It makes a handy catalog entry for individual galaxies; should we really get rid of the Ellipse class?

We don't really use it at all in Galsim for this kind of thing. We always treat the shear, magnification and shift separately. The only place we actually use it now as a way to implement applyDilation (and related), which is easily changed to a function that just does that directly.

It seems to me that the fact that we've been carrying around this functionality for a while now and no one has actually used it for something real probably means that it isn't as useful as we originally thought.

[reikonakajima]

Yes, I think you're right, you don't really need a galaxy catalog entry for GalSim.

[barnabytprowe]

Hi Mike,

I think your suggestions all sound sensible. I'd also like to see us providing library functions for converting ellipticity and sizes between coordinate systems, of course assuming local linearity: these are handy for even just catalogue handling but they are not implemented elsewhere as far as I've found in the standard WCS packages. I have code for doing this with ellipticities that I can donate to the branch (actually was going to run this by E. Sheldon for including in his esutil.py package too, so maybe we can share...)

I did have one question though: what is the timeline you have in mind? There is some fairly major surgery here that worries me if it were to mean that a lot of validation needed to be redone right before building the simulation images for GREAT3. I guess the ditching of the Ellipse in particular worries me on that score: the thing we need to get right is the shearing above all else, and it is adding an extra point of possible failure to overhaul this on master before the challenge launch. What do you think?

[rmjarvis]

I'd also like to see us providing library functions for converting ellipticity and sizes between coordinate systems

Good idea. I can do that.

what is the timeline you have in mind?

Soonish, but I don't think my changes would affect any validation tests. My plan wouldn't change the actual calculations that get done for (e.g.) applyShear. I just wouldn't go through CppElliipse anymore -- I'd go straight to SBTransform with the same math. If you look at the Ellipse getMatrix function, you'll see that the heavy lifting is already in the CppShear class, not in CppEllipse. All the latter does is scale the matrix by exp(mu). So it's easy skip the middleman here.

Plus, we have unit tests that test backwards-compatibility to high precision, so if I do mess something up in making the conversion, that will catch it.

[TallJimbo]

On the Ellipse class: if it's not needed, I'd be happy to take over the name and essentially replace it with some variation of my ellipses code from the HSC/LSST pipeline. That would be useful for some HSM changes I'd like to make in GalSIm (#368), and it may also be generally useful for users in helping convert between different ellipse parametrizations (the Shear class does the most important conversions already, of course, and I'd leverage those whenever possible).

[TallJimbo]

Oh, also, on the subject of transforming ellipticities between coordinate systems, I think the simplest approach by far is to convert back to moments (setting radius to unity if necessary), then transform the moments matrix (just matrix multiplications with the affine transform matrix), then convert back. So #367 may be helpful with implementing that.

[rmandelb]

Well, to be honest, I do find it slightly scary to take something that has meant one thing, and replace it with something similar but not quite identical. We'd be both (a) removing the shift and (b) changing the scale parameter in the ellipse, which I thought currently means amount of expansion/dilation to apply to an object, so that in the new version it will refer to the object's absolute size. I'm more worried about (b). For this reason I would kind of prefer to make a new Moments class which could have methods converting between size+e1+e2 and the moments matrix. But I guess it's a weak preference.

[rmjarvis]

I have an update on this issue. I've finished a fairly substantial refactoring of the code. I haven't added any new WCS classes yet, but I changed a lot of things to facilitate the addition of more WCS classes, so I figured it was worth reporting on what I've done so far to get feedback from anyone interested.

The main change I made is that I now have a single python layer Image class that has an image attribute and a wcs attribute. The former is one of the C++ image classes, and the latter is the WCS, which now encapsulates the functionality that had been the scale parameter.

This simplifies the use of Image objects considerably, since there is now just a single type for people to use. The data type can be specified as Image(..., dtype = numpy.float64) for example, or the old letter notation will also still work: ImageD(...). The latter is just an alias for the former now. The returned object is an Image instance.

All of the scale stuff in C++ has been ripped out and moved to the python layer. Specifically the image.wcs attribute. The draw command uses this to convert a profile from sky coordinates to chip coordinates before passing it on to the C++ draw command. For more complicated WCS types (which I haven't written yet), we will be able to use the existing SBTransform functionality to effect the conversion, which should be rather straightforward.

Here is a summary of the API changes so far that I put into the CHANGELOG:

• Changed the name of the dx parameter in the draw, drawShoot, drawK methods of GSObject and the constructors of InterpolatedImage and CorrelatedNoise to the name scale.
• Changed the xw and yw parameters of the Pixel constructor to a single scale parameter.
  • pix = Pixel(xw=scale) should now be either pix = Pixel(scale=scale) or simply pix = Pixel(scale).
• Added new Box class to take up the functionality that had been Pixel with unequal values of xw and yw.
  • box = Pixel(xw=width, yw=height) should now be either
    box = Box(width=width, height=height) or box = Box(width, height).
• Changed the dx_cosmos parameter of getCOSMOSNoise to cosmos_scale.
• Combined the old Image, ImageView and ConstImageView arrays of class names into a single python layer Image class that automatically constructs the appropriate C++ image class as an attribute.
  • im = Image[type](...) should now be Image(..., dtype=type)
  • im = ImageView[type](numpy_array.astype(type)) should now be im = Image(numpy_array.astype(type)). i.e. The data type inherits from the numpy_array argument when appropriate. If it is already the correct type, you do not need the astype(type) part.
  • im = ConstImageView[type](numpy_array.astype(type)) should now be im = Image(numpy_array.astype(type), make_const=True)
  • im = ImageF(...) and similar is still valid.
  • im = ImageViewF(...) and similar should now be im = ImageF(...) (preserving the same type letter S, I, F or D).
  • im = ConstImageViewF(...) and similar should now be im = ImageF(..., make_const=True) (again preserving the type letter).
  • im = ImageF(...) may now be written as im = Image(...). That is, the numpy.float32 type is the default data type if you do not specify something else either through the type letter or the dtype parameter.
• Changed the handling of the scale and init_value parameters of the Image constructor, so that now they have to be named keyword arguments rather than a positional arguments.
  • im = ImageF(nx, ny, scale, init_val) should now be im = ImageF(nx, ny, scale=scale, init_value=init_val).
• Removed the im.at(x,y) syntax. This had been equivalent to im(x,y), so any such code should now be switched to that.

Here are my plans for the WCS classes to implement: (Only PixelScale is implemented yet.)

• PixelScale encapsulates the normal behavior of using just a pixel scale. Most of the time you will not need to explicitly construct this kind of WCS object, since the notation of ImageF(..., scale=pixel_scale) or im.scale = pixel_scale will set this up automatically for you.
• ShearWCS is a sheared coordinate system, where the shear is constant across the field.
• JacobianWCS has an arbitrary 2x2 matrix to convert from (x,y) to (u,v):
  • u = dudx x + dudy y
  • v = dvdx x + dvdy y
• AffineTransform is the most general constant WCS transformation. It involves a jacobian and an arbitrary shift:
  • u = dudx x + dudy y + u0
  • v = dvdx x + dvdy y + v0
• FitsWCS reads WCS information from a FITS header. It requires astropy (or possibly whatever preceded this as a python wrapper to wcstools) to be installed.
• TanWCS reads either 'TAN', 'TNX' or 'TPV' WCS information from a FITS header file. These are probably the most common WCS types, and I have code that can parse them, so I thought it would be useful to let these have a native GalSim implementation that does not require astropy.

And here is the planned functionality for the WCS classes:

• sky_pos = wcs.toSky(chip_pos) converts a position in chip coordinates to sky coordinates.
• chip_pos = wcs.toChip(sky_pos) converts a position in sky coordinates to chip coordinates.
• area = wcs.pixelArea() returns the area (in arcsec**2) of a pixel.
• wcs.applyTo(profile) converts the given profile (a GSObject) from sky coordinates to chip coordinates.
• wcs.applyInverseTo(profile) converts the given profile from chip coordinates to sky coordinates.
• local_wcs = wcs.local(chip_pos=chip_pos) returns the local linear approximation of the WCS at a given point. Alternatively, sky_pos may be provided if that is more convenient.
• local_affine = wcs.localAffine(chip_pos=chip_pos) returns the local linear approximation of the WCS at a given point when you need that returned value to be an AffineTransform type. For variable WCS, this will be the same as wcs.local(...), but local is allowed to return any constant WCS type, which might be more efficient in some cases. e.g. PixelScale can return itself.

If the wcs is not constant, you need to provide either a chip_pos or sky_pos argument to many of the above functions to tell it where you want the pixel area or the profile drawn or whatever. e.g. the following commands are equivalent:

• area = wcs.pixelArea(chip_pos = chip_pos)
• area = wcs.local(chip_pos = chip_pos).pixelArea()
• area = wcs.localAffine(chip_pos = chip_pos).pixelArea()

I think I've converted almost all the places that used the pixel scale to use one of these functions instead. I believe there are only a couple of places where I haven't done so, where I instead use a wcs.linearScale method which returns the pixel scale for PixelScale. For these places, I need to think more to figure out what the right thing to do in those cases is when the pixels are not square.

Also, there are a few places where using a wcs other than PixelScale raises a NotImplementedError. In those cases, I think the implementation will require a separate issue. These are:

• Adding correlated noise to an image with a non-trivial WCS. I think this is possible, but it requires some thinking. (And probably help from @barnabytprowe!)
• Fitting a shapelet profile from an image with a non-trivial WCS. This is definitely possible, since I have code to do it elsewhere, but it's a largish job to incorporate it into GalSim, so I'll do that on another issue.

If you got through this long post, thank you. And I would be grateful for any comments that you have about any of this. :)

[msimet]

Just to check that I'm understanding--you've combined scale and dx into a single scale parameter and are planning to use scale to denote an arbitrary(ish) WCS?

[msimet]

I should say--I did get through the post and that all looks really good--just wanted to check I understood what was going on with the kwargs!

[TallJimbo]

I've just given this a cursory look so far, and I like how you're handling the Image classes in Python. I also think your plans for WCS sound reasonable, but you might want to consider plugging into AstroPy for standard WCS-handling - it's rapidly becoming a standard Python package for astronomers, to the point where if we were starting the LSST stack from this point I think we'd build on top of it for a lot of things. GalSim is young enough that it might still be able to without a lot of refactoring (the LSST stack is not). I realize that potentially adding a new dependency is a much larger question, but it's quite relevant here as I think WCS is one of the main things they've implemented that we might want.

[rmjarvis]

Just to check that I'm understanding--you've combined scale and dx into a single scale parameter and are planning to use scale to denote an arbitrary(ish) WCS?

No. There was never any place that used both dx and scale. Rather, we had three parameter names in different places that all basically meant the pixel scale. In Image, we used scale, in draw and others we used dx, and in Pixel we used xw. Now all three use scale to mean the pixel scale.

However, image.scale cannot be an arbitrary (even -ish) WCS. It is only used when you want the WCS to be a trivial pixel scale conversion. Anything more complicated is set as the image.wcs attribute. It can be set via a wcs parameter to either the Image constructor or the draw command (and similar commands). Or you can explicitly set im.wcs = whatever after an image is constructed the same way you had been able to (and still can) write im.scale = 0.2.

Basically, the scale parameter is shorthand for im.wcs = galsim.PixelScale(scale). In fact, that is what it is doing behind the scenes via a python property (it's not technically an attribute). But you can treat it like an attribute the same way we had used im.scale before.

Hopefully this clears things up for you. :)

[rmjarvis]

you might want to consider plugging into AstroPy for standard WCS-handling

Thanks Jim. I had planned to use AstroPy for the more complicated WCS handling. Especially the WCS read in from FITS headers. (Called FitsWCS above.)

I haven't looking into it too closely yet, so maybe I should do that before I get too far along in implementing some of the other classes I had planned. At the very least, it might be worth trying to imitate some of their style to make it easier for people who already know the syntax in the astropy.wcs module.

[rmjarvis]

I looked at astropy.wcs, and as I suspected, it is pretty much just for dealing with WCS projections that are in FITS headers. So not the simpler kinds of things that I think we also want here. But it will be useful as a backend for the FitsWCS class that I was planning to write.

Also, it still doesn't have the TNX or TPV projections, which have been fairly standard for years now. (At least the former -- I wrote the first version of my TNX parser in grad school.) So we should have a class that can handle those. Any preference on whether we combine them into a single FitsWCS class or have two separate classes?

Advantage of combining them: Users don't have to know that astropy.wcs can't handle the TNX projection. They can just use the same class for all FITS files.

Advantage of not combining them: Users don't have to have astropy installed if they are only going to use FITS file with TAN, TNX or TPV projections. They can just read them with the TanWCS class.

Finally, a syntax question: astropy.wcs uses the terms "world" and "pix" where I used "sky" and "chip". Any preferences between these options? It's easy to switch now if people have a preference.

[rmandelb]

Hi Mike - sorry I've been slow at commenting; I am visiting Princeton and have a lot of meetings. I will try to look it over and comment this week. However, a question came up today that might be related so I wanted to ask you about it on here:

Jim and I were talking about how best to support making output files obeying the geometry of some camera, so that you can make fake HSC / DES / LSST / etc. data (presumably with some code in a galsim.hsc etc. module). Once you have the WCS implemented, would it be sensible to say something like "the camera has N chips, and here is a list of the WCS for each chip" as a way of defining the geometry?

[rmjarvis]

would it be sensible to say something like "the camera has N chips, and here is a list of the WCS for each chip" as a way of defining the geometry?

Clearly the WCS is the main thing you would need for that. Once this is in place, it would not be too hard to have a set of images with compatible WCS solutions that make up an array of CCDs in a focal plane. We could even have some survey-specific modules that know how the chips are arranged for that camera and have them be set up automatically. I'm thinking of a function like image_list = galsim.hsc.buildFocalPlaneImages(ra=ra, dec=dec).

A more interesting capability that we could add along this lines is to have GalSim know that a set of images corresponds to a single focal plane. Then we could have functionality where you just tell GalSim to draw an object in that focal plane and it would automatically figure out in which image the object belongs. That might be a handy addition that would facilitate building full exposures that include a bunch of CCD images. Something to think about...

[barnabytprowe]

Hi Mike,

I've been digesting this and I think the API looks very good. I so far have a couple of randomly scattered thoughts, questions, comments...

Adding correlated noise to an image with a non-trivial WCS. I think this is possible, but it requires some thinking. (And probably help from @barnabytprowe!)

Actually I think this will not be too bad... This only makes sense for the "constant" WCS transformations, as you term them. And in fact the full AffineTransform is also not needed, as the noise correlation functions assume translational symmetry anyway. In which case, we merely need to split up the JacobianWCS part of the WCS into equivalent dilations, rotations and shears... And possible mirror image inversions. All of which bar the last are already implemented in the CorrelatedNoise profile. So I think we just need to work out the right theta, g and mu given a set order of applying each transformation, and make a mirrorFlip or similar method for GSObjects to handle cases where the Jacobian determinant is negative. As you know, I'd been thinking about this stuff recently anyway!

• wcs.applyTo(profile) converts the given profile (a GSObject) from sky coordinates to chip coordinates.
• wcs.applyInverseTo(profile) converts the given profile from chip coordinates to sky coordinates.

I confess, I find those method names rather uninformative. What would you say to

wcs.fromSkytoChip(profile)
wcs.fromChiptoSky(profile)

instead? Saves having to remember which is the fwd/inverse direction.

im = ImageF(...) may now be written as im = Image(...). That is, the numpy.float32 type is the default data type if you do not specify something else either through the type letter or the dtype parameter.

I know that this makes sense for much astro data, and for huge sims like GREAT3, but in general usage doesn't a default float size of single precision seem a bit... I dunno, 1980s??? ;) I admit a preference for doubles for general use, and it's also NumPy's default float size, although there are many NumPy defaults/conventions that are inconvenient I admit.

Anyway, the fact I bring up such minor queries is testament to how well thought-through I think this looks (as well as being achievable).

[barnabytprowe]

Also, it still doesn't have the TNX or TPV projections, which have been fairly standard for years now. (At least the former -- I wrote the first version of my TNX parser in grad school.) So we should have a class that can handle those. Any preference on whether we combine them into a single FitsWCS class or have two separate classes?

Advantage of combining them: Users don't have to know that astropy.wcs can't handle the TNX projection. They can just use the same class for all FITS files.

Advantage of not combining them: Users don't have to have astropy installed if they are only going to use FITS file with TAN, TNX or TPV projections. They can just read them with the TanWCS class.

Surely it's possible to have the best of both worlds, i.e. only import astropy if needed, e.g. using one of the astropy-supported projections that aren't TAN, TNX or TPV?

Finally, a syntax question: astropy.wcs uses the terms "world" and "pix" where I used "sky" and "chip". Any preferences between these options? It's easy to switch now if people have a preference.

I was thinking about that when I read the original post. I don't really have a preference for either, but I do see a benefit in not using different names when we mean the same thing. If it's easy to switch now, I'd vote to do that.

[rmandelb]

Re: world/pix vs. sky/chip - I personally have a slight preference for the former, and I also like the idea of consistency with astropy given its current popularity.

[rmandelb]

Hi Mike -

Thanks for the update on your progress on this issue. Most of it looks good; for example, I worried slightly when you mentioned the refactoring of the Image class but after reading the details, I think that the new version is likely to be more user-friendly. There are just two cases when that’s not true, but I might have misunderstood, so here are my questions about those two cases:

(1) Sometimes we wanted to make sure to have an ImageView, so we would have something (im) that could have been either an Image or an ImageView, and we would ensure we had an ImageView using im.view(). Is it now the case that we’d have to do galsim.Image(im.array)? This seems a bit clumsier / less obvious what’s going on. Or does the refactoring mean it’s no longer necessary to do something like this at all (i.e., the former distinction between Image and ImageView is now unnecessary)?

(2) If an Image now has an image attribute, then for a galsim.Image called im, do all of our lines of code that were im.array have to turn into im.image.array?

No concerns about your WCS plans - it will be good to have that functionality in GalSim.

[rmandelb]

About having survey-specific modules that know how chips are arranged for a particular camera: I wonder if we could make a module that can talk to the LSST stack (for someone who has that installed) and can use the butler to learn about camera geometry, given that the LSST stack has information about a number of different cameras, and could end up with information about more as more people adopt it for regular use.

Hopefully Jim can comment on whether this is an insane idea.

Of course, people may want to have their own galsim.surveyX modules for their surveys that have their own buildFocalPlaneImages methods without requiring the LSST stack.

And I like your idea of being able to simulate a single focal plane with many chips, and have GalSim figure out which image a particular object lives in. Sounds like it could be a fairly complicated thing to implement in a general way.

[rmjarvis]

Thanks for the comments. I'll try to respond to them all here.

Or does the refactoring mean ... the former distinction between Image and ImageView is now unnecessary?

Yes and no. The end user should never need to know about the difference. They can always just use the Image class as is and never worry about the fact that there are actually 12 different C++ image classes behind the scenes.

But in the C++ layer, we still do have the distinction, because C++ needs to treat them differently. So in our python functions that do call down to a C++ function, you still might need to make a conversion to get the right type. c.f. the code in hsm.py where I create image_view objects to pass to C++ functions.

do all of our lines of code that were im.array have to turn into im.image.array?

No. im.array is a property which returns im.image.array. A bit of python decorator magic. :)

Surely it's possible to have the best of both worlds, i.e. only import astropy if needed, e.g. using one of the astropy-supported projections that aren't TAN, TNX or TPV?

I suppose that could work. It would be a bit cleaner to have them be separate classes, since the implementation of all the functions will be quite different. But I could have FitsWCS be a function rather than a class that would check the projection type and return either an AstropyWCS or a TanWCS as appropriate. Then subsequent operations would work correctly using the appropriate class methods.

And the end user would never need to know that this was the way it worked unless they investigated more closely. As far as they were concerned, FitsWCS would work just as though it were another WCS class.

Re: world/pix vs. sky/chip

I will point out one difference in meaning between what wcstools calls "world" and what I have been calling "sky". The "world" coordinates are RA, Dec. But "sky" coordinates are really a rectilinear coordinate system on a tangent plane projection.

These are subtly different things, and it is important, because the galaxy profiles are really defined with respect to a rectilinear coordinate system where arcsec mean the same thing in both directions (as opposed to having to put in a cos(dec) term in the ra direction). So it might be confusing to use the term "world" when we mean something slightly different from what wcstools (and astropy.wcs) means by "world".

To be honest with myself, I suppose I am resistant partly because I don't really like the term world (a synonym of Earth) to mean the coordinates in the sky (not on Earth). But I do see the benefits of conforming to external standards, so I'll try not to be quixotic about this. :)

wcs.fromSkytoChip(profile) wcs.fromChiptoSky(profile)

I agree with your objection. I also didn't really like the names applyTo and applyInverseTo, so I'm glad you called me on it. How about this: I already have chip_pos = wcs.toChip(sky_pos) (and the converse). I could use the same name to convert the profile as well: profile_in_chip_coords = wcs.toChip(profile_in_sky_coords).

Would that be confusing to have one name do basically the same thing to two different kinds of objects? In C++ I would use overloading, but in python I would implement this with an isinstance check for the two cases.

in general usage doesn't a default float size of single precision seem a bit... I dunno, 1980s??? ;)

I don't think so. I normally also have a preference for double most of the time, but I think this is an exception. It is very rare for an image to require double precision. The noise level that we add to the image is usually quite a lot higher than the numerical noise of float32. So all the extra digits for float64 would be wasted bits.

And as you pointed out, it really is quite a bit of savings when building a large suite of images for testing as we did for Great3, which I think/hope will be a common use of GalSim. There is a reason ImageF has been by far the more common choice in our tests and demos. In fact, none of our demos even use ImageD. We mention it once in demo3 to let people know about it, but we don't actually use it anywhere.

[rmjarvis]

I've been thinking about the names to use for the different coordinates, and another package that might be worth conforming to is Sextractor. It uses world and image. Again with the world name, so I think (despite my misgivings) we should use that for what I have been calling sky coordinates.

Furthermore, sextractor has things like X2WORLD and such which are measured in the local tangent plane coordinate system, not in the non-Euclidean distorted RA/Dec system. So I think it will be ok to use world for that same meaning in GalSim when we talk about the profile in world coordinates.

But I think I like "image" better than "pix" for what I have been calling chip coordinates. Mostly because I want to use it as a prefix to _pos, and image_pos sounds better than pix_pos. (We already use image_pos in the config to mean the position on an image, so this would also conform with our own prior nomenclature.) So would people be ok with: world and image for the names of the two coordinate systems?

The other thing I've been thinking about here is how best to handle WCS projections which are inherently in terms of RA/Dec. I think it would be ok to have the positions returned by toWorld(image_pos) be (ra, dec), but have the effect of toWorld(image_profile) be to return a profile in the local tangent plane system (which I call (u,v)).

Even though that sounds like a mismatch, I think it is actually the most intuitive of the various options I've considered, and it is what sextractor does for its second moment (and other) measurements. You never really want the profile to be in terms of ra, dec, since that doesn't really make sense. And for the world coordinates, it's sub-optimal to have a single tangent plane projection for the whole image -- you really want to always take the local tangent plane at the location of each object you consider. So that means there isn't really a well-defined (u,v) position for the toWorld(image_profile) to return.

However, I was thinking that to make this work well, I might want to write a CelestialCoord class to hold the ra, dec values rather than return them as a Position2D, since the arithmetic operations on a Position2D don't really work correctly on the sphere. Then when a WCS class is really returning an ra, dec value, or when it needs one as an input, it would use this class. I think this would help the user keep track of the coordinate systems in the two cases.

[rmjarvis]

I'm up to the point where I am writing the code to read and write all my different WCS classes to FITS headers, and I have a question, probably for @TallJimbo (since I think he is the original author of this bit), but anyone else can weigh in too.

Our galsim.fits.write function has an optional add_wcs parameter with the following documentation:

    @param add_wcs      If `add_wcs` evaluates to `True`, a 'LINEAR' WCS will be added using the 
                        Image's bounding box.  This is not necessary to ensure an Image can be 
                        round-tripped through FITS, as the bounding box (and scale) are always 
                        saved in custom header keys.  If `add_wcs` is a string, this will be used 
                        as the WCS name. (Default `add_wcs = True`.)

I'm guessing this was written with the idea that it might be useful, but in practice there is no place in any example or unit test where we set add_wcs to anything. Either False to turn it off, or to some string so the WCS keys have non-default names.

So do we need this? Can I get rid of this option and just always write the CTYPE, CD, CRPIX, etc. keys?

It would make my current job a bit easier, since I'm not sure what this keyword means when the underlying WCS type requires these keys. (e.g. the AstropyWCS class which reads this kind of stuff from a FITS file in the first place.)

[TallJimbo]

On add_wcs: I think you're right, it was just added in the hopes that it would be useful, with no particular use case in mind. If you haven't seen anything using it, I have no problem with removing it.

[rmjarvis]

On add_wcs: I think you're right, it was just added in the hopes that it would be useful, with no particular use case in mind. If you haven't seen anything using it, I have no problem with removing it.

Thanks. That's what I suspected, but I wanted to confirm.

[rmjarvis]

Update: I think I've pretty much finished all the python-level functionality that I want to do. I still have to put this into the config apparatus, but I'm going to take a break from this for a bit to get to a few other things I need to do. (The great3 PR I promised Rachel, for one.)

But if anyone would like to peruse this, here is a summary of what I've done, which differs a bit from my last summary (due to the helpful comments I received then -- thank you!).

This is copied directly from the docstring of BaseWCS:

    There are two types of WCS classes that we implement.

    1. Local WCS 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 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 local WCS classes:

            PixelScale
            ShearWCS
            JacobianWCS

    2. Non-local WCS classes may have a constant pixel size and shape, but they don't have to.
       They may also have an arbitrary origin in both image coordinates and world coordinates.
       Furthermore, the world coordinates may be either a regular Euclidean coordinate
       system (using galsim.PositionD objects for the world positions) or coordinates on
       the celestial sphere (using galsim.CelestialCoord objects for the world positions).

       Currently we define the following non-local WCS classes:

            OffsetWCS
            OffsetShearWCS
            AffineTransform
            UVFunction
            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 is also a factory function called FitsWCS, which is intended to act like a 
    class initializer.  It 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 it's 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.

    Some things you can do with a WCS object:

    - 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 galsim.PositionD instance.  However, world_pos will
      be a galsim.CelestialCoord if the transformation is in terms of celestial coordinates
      (c.f. wcs.isCelestial()).  Otherwise, it will be a PositionD as well.

    - Convert a GSObject, which is naturally defined in world coordinates, to the equivalent
      profile using image coordinates (or vice versa):

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

      For non-uniform WCS types (c.f. wcs.isUniform()), 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.

    - 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.

    - Shift a transformation to use a new location for what is currently considered
      image_pos = (0,0).  For local WCS types, this also converts to a non-local WCS.

                world_pos1 = wcs.toWorld(PositionD(0,0))
                shifted = wcs.setOrigin(image_origin)
                world_pos2 = shifted.toWorld(image_origin)
                # world_pos1 should be equal to world_pos2

    - 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

      Global WCS objects 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 a PixelScale or OffsetWCS?