LSST module with survey specific information

[rmjarvis]

It would be nice to have an lsst module in GalSim that people could access to quickly access some of the information about the camera geometry, bandpass information, etc.

For example, it would be nice to be able to get the appropriate WCS function for a given chip with some simple call like

wcs = galsim.lsst.getWCS(chip = 132, ra = ra, dec = dec)

where ra, dec would be the pointing of the field center.

Likewise the bandpasses for the different filters that will be used (I think these are already the ones Josh loaded into examples/data, so this might already be ready to go. I'm not sure what else would be useful, but this is a good start.

And of course, any other survey that wants to create something similar should do so as well. But getting one implementation done will help others see what they would need to do. Particularly w.r.t. the config interface. (I'll help with that part of course!)

[HironaoMiyatake]

I would like to mention that this will help HSC as well, since the HSC pipeline is based on the LSST pipeline.

[rmjarvis]

Except HSC will have different chips, pixels, etc. So I don't think it would translate.

However, once we solve the functionality issues with getting this set up for one survey (HSC or LSST or DES or ...), I think that will make it much easier to get similar things for the other surveys. So if you want to do this in the context of HSC first, that would be fine too. I picked LSST for the title semi-arbitrarily based on the fact that Jim suggested it in the meeting today.

[HironaoMiyatake]

Both HSC and LSST are based on the same code and just configuration files are different. So once we write a code to read a configuration file, I think it would translate. Jim can confirm this.

About the development, I think we should begin with LSST. According to Jim, the code and config format are now being updated on the LSST side and it will be copied over the HSC pipeline later. So there is no need to work on HSC now.

[mdschneider]

This issue is now relevant for a near-term Deliverable in the LSST DESC Science Roadmap. We would like an LSST module to generate a data set for the DESC 'Data Challenge 1' to perform controlled studies of shear systematics.

Some desired features from the DESC Roadmap include:

1. Model the chip layout and WCS (as mentioned in @rmjarvis comments above)
2. Model the field-dependent optics PSF
3. Model atmosphere PSFs for 15-sec exposures (i.e. with more anisotropy than allowed by an elliptical Kolmogorov profile)
4. Model sensor systematics such as the Brigher-Fatter effect
5. Store a pre-computed spider image for LSST

@danielsf has an implementation of the first issue above.

Pinging also @rmandelb @TallJimbo @jmeyers314

[mdschneider]

This is a place to start for addressing this issue: https://github.com/lsst/sims_GalSimInterface

[mdschneider]

Pinging also @dkirkby

[dkirkby]

It looks like LSST DM is approaching this problem somewhat orthogonally: they are implementing a module lsst.sims.GalSimInterface that interfaces GalSim to LSST rather than providing a lightweight standalone definition of LSST chip geometry, etc, that could be wrapped in a galsim.lsst module. Their primary use case is feeding DM from PhoSim/GalSim/data, whereas the use case here is simulating the same catalog as it would appear in different surveys.

Both approaches are probably needed, and should ideally share as much code as possible. However, it sounds like the current DM stack has the chip geometries, etc, entangled with lots of other stuff in some large packages that we don't want GalSim to depend on. I suspect they would be willing to refactor to clean up the dependencies for our use case, and that would be preferable to GalSim duplicating this code.

[rmjarvis]

Model sensor systematics such as the Brigher-Fatter effect

FYI, this is already in GalSim. #524

[rmjarvis]

Model atmosphere PSFs for 15-sec exposures (i.e. with more anisotropy than allowed by an elliptical Kolmogorov profile)

This one is not yet implemented, but I think it is more properly done as part of #549. This has applicability beyond just LSST.

[rmjarvis]

Model the chip layout and WCS (as mentioned in @rmjarvis comments above)

As for the chip geometry stuff, I think we should try to get LSST DM to create from their code some output files with the necessary information that we could read in from GalSim and create GalSim-usable WCS's. If we can standardize on the format, then it's fine if they improve their geom module. They can just make new output files and we can update them in the GalSim repo.

I doubt we will want to port over large portions of the geom package into GalSim. That just seems like the wrong approach to me.

[rmjarvis]

Model the field-dependent optics PSF

This strikes me as potentially a lot of work to get right. We might be able to get estimates of the mean Zernike coefficients as a function of position on the fov. If so, then that plus the spider would probably suffice for most purposes here.

But to accurately model the ways that this pattern may vary from exposure to exposure would be pretty hard I think. So for now, we should probably limit the scope of this to the static mean pattern of the optical PSF.

[rmandelb]

As for the chip geometry stuff, I think we should try to get LSST DM to create from their code some output files with the necessary information that we could read in from GalSim and create GalSim- usable WCS's.

I agree.

We might be able to get estimates of the mean Zernike coefficients as a function of position on the fov. If so, then that plus the spider would probably suffice for most purposes here.

Yes, I agree. This is at least the place to start for any short-term work.

[rmandelb]

Re: brighter-fatter - true that it's in GalSim, but it would be appropriate for the LSST module to know something about the magnitude of the effect for LSST so that it can automatically initialize an appropriate BaseCDModel.

[rmjarvis]

it would be appropriate for the LSST module to know something about the magnitude of the effect for LSST so that it can automatically initialize an appropriate BaseCDModel.

Good point. Yes.

[danielsf]

Just to offer a perspective from within the LSST Simulations team:

@dkirkby is correct. The current implementation of the CatSim-GalSim interface does use the cameraGeom software provided by the DM afw package to take positions on the sky and transform them into positions on the camera and on individual detectors. There are some advantages and disadvantages to this.

Advantages: 1) The framework is actually LSST-agnostic. afw.cameraGeom is a software package that reads in specifications from a text file (or a FITS file, or any other number of input formats) and creates a self-consistent model of all the coordinate transformations associated with a sky-to-camera system. In addition to detector layout, this also includes pincushion and (I think) Zernike distortions. Right now, the API is a little cumbersome, but we are working on python scripts that will make it fairly easy for users to generate an arbitrary camera layout using a simple text file outlining the positions and sizes of chips and pixes.

2) Because of (1), if the engineering team ever changes the specs on the LSST camera, the software will automatically be in sync (assuming it gets updated with the newest distribution of the relevant packages). We do not need to worry about passing around GalSim-specific text files every time something changes at LSST.

Disadvantages: 1) afw is a C++ package that takes about 45 minutes to compile on a modern laptop and occupies 216M of disk space.

2) (1) isn't that bad, but afw brings with it dependencies on a lot of other, unnecessary DM software packages. The Simulations team would definitely like to split cameraGeom off from all of this baggage and make it much more portable. That is a conversation that will have to involve the DM team, since they are more aware of how feasible such a scheme is.

If continuing down this path does not seem palatable, other options are:

1) Write a new python package that re-implements the functionality of afw.cameraGeom without pulling in all of the unnecessary dependencies. As I said, the LSST Simulations team is not terribly pleased with our own dependence on afw, so we are open to writing a new package which the LSST Simulations stack will use to implement camera geometry. We can hand that package over to GalSim so that, at least, GalSim and LSST Simulations will remain in synch. This has the disadvantage that we will have to make sure that sims.cameraGeom and afw.cameraGeom remain in sync, but that is sort of in our (LSST Simulations) job description, anyway.

2) GalSim implements its own LSST camera module and we (LSST Simulations) become responsible for updating the inputs to that code whenever the specifications for the LSST camera change. There's a lot of risk for things to get out of sync with this scheme. My preference would be to come up with a single software package that both GalSim and LSST Simulations depend on for simulating the LSST camera.

We can talk about this in person next week during the WL2 session. This is where my/our (I'm not 100% sure who all I speak for) thinking is at the moment.

[rmjarvis]

Unsurprisingly, I prefer option (1), since it's less work on our end. Especially if this is truly a standalone package that people could just pip install. Then rather than hand it over to GalSim, we just say that the galsim.lsst module depends on the python package lsst_camera_geom (or whatever) and we just import it and use it within galsim.lsst.

[TallJimbo]

I think it ought to be possible to refactor the core of cameraGeom into a lightweight package that GalSim, LSST Sims, and DM could all depend on. DM's cameraGeom in afw would then just provide interoperability with other afw classes.

The main difficulty is that DM will need access to at least some functionality in C++, and we don't presently have any C++ code in the stack that delegates to Python. We've been thinking about doing that, so it's not out of the question, but it might make this project a bit harder if it's the first such case. Another option would be to implement at least some of the lightweight camera package in C++ as well, but not necessarily use Swig or the any of the DM approach to providing Python bindings (e.g. you could use Boost.Python, Cython, or the raw Python C API - whichever was easiest for this problem). DM could then depend on the C++ interface from the lightweight package and then either add its own Python bindings or find a way to make the lightweight package's bindings work with Swig.

[rmjarvis]

You could have the python package call C functions, which are also accessible from your C++ layer. It's very easy to wrap C functions using ctypes. Then they would be callable from both python and C++.

[jmeyers314]

Likewise the bandpasses for the different filters that will be used (I think these are already the ones Josh loaded into examples/data, so this might already be ready to go.

I believe the throughputs I added a while ago are specifically for airmass=1.2. I think it's worth putting in a basic atmospheric extinction function (even if we neglect the specifics of water vapor, aerosols, etc.) to get this right as a function of airmass.

Relatedly, I think it makes sense to specify the atmospheric PSF FWHM specifically for observations at zenith and at 500 nm, and then scale to whatever elevation and wavelength (if using the chromatic code) or filter effective wavelength (if rendering achromatically; this should probably be an option for speed) is requested.

[rmjarvis]

We decided that for now we will make an LSSTWCS class that will try to import lsst.afw.cameraGeom and give a reasonable error message if the stack is not installed.

It will then implement all the required GalSim WCS functions as thin wrappers to the corresponding functionality within cameraGeom:

• _radec(x,y) returns the corresponding ra,dec for a given (possibly lists) x,y
• _xy(ra,dec) does the reverse
• _writeHeader should write something to a fits header that would allow the object to be correctly reconstructed from reading the header. Think serialization. This could probably write the TAN-SIP approximation and also some proprietary keys that have an exact serialization.
• _readHeader should parse the serialization version stored in the header to reconstruct the original wcs.

In addition, some boiler plate:

• _newOrigin(origin) sets the internal origin parameter, which gets returned by self.origin. This doesn't need to interact with cameraGeom.
• copy
• __eq__
• __repr__ should produce a unique serialization
• __str__ should provide something more readable
• __hash__ is typically just hash(repr(self))
• __setstate__, __getstate__ if necessary to make it picklable.

[danielsf]

I am trying to decide what the unit test should be that verifies the

(RA, Dec) -> name of chip that sees it

transformation. In the LSST stack, we create a cartoon camera with no distortions so that we can analytically compute how RA, Dec map to position on the camera and test on that. That seems like it would involve porting in more LSST code than we want for this use case.

I can generate a bunch of examples using the LSST stack, store them as a text file, and just have the unittest validate those examples. Then, the onus would be on us (LSST) to update that text file if our camera model ever changes. Does anyone object to this idea (or have a better one)?

[rmandelb]

That sounds like what I would do.

As a reminder, you could have comments in the camera model routine in GalSim (for developers - not in the docstring but in the code) saying things like "WARNING: IF YOU CHANGE THESE LINES, YOU HAVE TO UPDATE THE UNIT TEST CALLED BLAH BLAH IN TEST_LSST.PY".

[rmjarvis]

This works as a regression test. So basically, the test will check that the galsim.lsst module gets the same answer that the lsst dm stack gets. We'll be relying on the DM code being right for this, but we're already counting on that anyway, so I'm fine with that. :)

Also, a good place to put the code that runs the DM stack and generates this file would be in devel/external. Then reference that script in the test suite so developers will know how to recreate the text file if necessary.

[danielsf]

@rmjarvis are _xy(ra, dec) and _radec(x, y) supposed to accept/return plain floats for ra, dec, or galsim angle objects?

[rmjarvis]

@rmjarvis are _xy(ra, dec) and _radec(x, y) supposed to accept/return plain floats for ra, dec, or galsim angle objects?

floats in all cases. They will have bee converted from Angles to radians upstream. (cf. _posToImage in wcs.py, which is the function that will call _xy)

[danielsf]

@rmjarvis Does the origin of the pixel coordinate system represented by _xy() need to correspond to the CRVAL1, CRVAL2 sky coordinates from the FITS header?

[rmjarvis]

Not normally. In a regular FITS header, there are also CRPIX1, CRPIX2. IIRC, these are the (x,y) coordinates that correspond to the sky coordinates given by CRVAL1, CRVAL2.

The origin (x,y) = (0,0) is normally the origin of the chip coordinate system. Typically, slightly off the lower left corner, since the center of the lower left pixel is normally defined as (1,1).

If you have an image which has both a normal fits header and also a description using cameraGeom, you can probably figure out the appropriate conventions by comparing what values for (x,y) cameraGeom uses vs what ds9 shows. I did a lot of that kind of comparison in writing the other wcs classes. The unit tests in test_wcs.py reflect this, comparing the values we get for various wcs back ends with what ds9 showed for a few choice (x,y) positions for a number of different WCS types.

[danielsf]

@rmjarvis What is eq used for? More specifically, should it return 'True' if two WCSs are equal to within some rounding error? Currently, outputting the WCS to a fits image and reading it back in produces WCSs that are inequal due to rounding error. Is this acceptable behavior, or not?

[rmjarvis]

__eq__ should return whether two objects are precisely equal, not just close. It's probably not "used for" much -- we define it for all GalSim classes to override the default python implementation, which is equivalent to is.

The round trip through the fits header doesn't need to be perfect. It should be possible to make it accurate to pretty good accuracy though, possibly by writing extra (proprietary) header keys. I've been starting them with GS_ when I do this. But even then, we're always limited by the finite number of digits stored in fits header values, so there's only so much we can do.

Once you get most of the functionality set, I can take a look and see if I want to tighten that up. But don't worry about it too much now if you've already got it reasonably close for a round trip.

[rmandelb]

An update that is relevant to this issue:

@mdschneider has been working in #716 on the field-dependent optical PSF aberrations, and I've started looking into a different optical PSF-related issue -- representing the obscuration and struts for LSST. There are three parts to this:

• Just putting in the right obscuration value (we can start the lsst module and tell it this information, i.e., there will be a parameter called galsim.lsst.obscuration, similar to what we have for WFIRST).
• Including the right pattern of struts - for this I am waiting on some info from @aaronroodman , but it's possible that this can be done using a simple variant on our current model for struts (without having to load a full pupil plane image).
• Incorporating field-dependence, i.e., taking into account the fact that the obscuration disk is only perfectly centered if you're at the center of the FOV. I was thinking we might ignore this as a higher-order effect, and begin by only including field-dependence of the aberrations.

[rmandelb]

Aaron sent me an image that shows the geometry of the spiders:

So I'm thinking that these could actually be handled directly using kwargs (instead of by loading a pupil plane image, which is more expensive). The way to do this would be to use the current code with 4 equally-spaced radial struts, but allow for the inner ~90% of the struts to be missing and hence that area is unobscured. We'd have to add a kwarg for what fraction of the strut is missing, but I think I could code that up pretty easily.

@rmjarvis , @jmeyers314 , any thoughts on this idea?

[rmandelb]

@aaronroodman , do you have any updates on this (the info I would need to understand the width / spacing of the spiders)? Do you prefer to wait until after the meeting in May?

[rmjarvis]

Leaving this open, since the PSF stuff hasn't been done yet. (E.g. the spider stuff Rachel mentioned above.)

But PR #690 is merged now with the WCS functionality. (Albeit requiring the LSST stack to be installed in order to use.)

[rmandelb]

For getting the spider right, I found that Aaron did send me the details needed to do the treatment I had proposed:

So I'm thinking that these could actually be handled directly using kwargs (instead of by loading a pupil plane image, which is more expensive). The way to do this would be to use the current code with 4 equally-spaced radial struts, but allow for the inner ~90% of the struts to be missing and hence that area is unobscured.

For reference, the info he sent was as follows:

"The spiders are each 50mm wide and the distance between the spiders is 80cm — that is the distance from the midline of one to the midline of the other along the 45deg line.

Also note that the central obscuration is 61%, that is a radius of 2.55 meters compared to the primary’s radius of 4.18 meters (thats clear area, so little under the 4.2 m glass radius)"

So in principle, I could make it possible to just remove the inner X% of each strut (i.e. make that area unobscured) using the figures given by Aaron. Or use a high-res FITS file, which has a bit of work needed to make compliant with our format, but not much.

Anybody have any thoughts on these options? My only concern with the former is that since the needed resolution is calculated based on the aperture + obscuration ignoring the demands of representing struts, and the spider is very narrow, the first option will require careful manual tuning of the oversampling to make sure we're getting it right.

[jmeyers314]

My only concern with the former is that since the needed resolution is calculated based on the aperture + obscuration ignoring the demands of representing struts, and the spider is very narrow, the first option will require careful manual tuning of the oversampling to make sure we're getting it right.

I think that's equivalent to manually tuning the size of the FITS image for the latter option, though.

Maybe the clearest thing to do would be to use the pupil_plane_scale and pupil_plane_size kwargs to explicitly set the geometry of the pupil, rather than through the aperture + obscuration calculation.

[rmandelb]

I think that's equivalent to manually tuning the size of the FITS image for the latter option, though.

Yes, I agree it is. But we can make the FITS image that we distribute a reasonable default size, so that end users should not have to do much/any tuning.

I think you might be right about the cleanest way to do this.

[aaronroodman]

Rachel & Josh, I think the question is: is there enough resolution to give a nice set of diffraction spikes? and will the calculation be fast enough with that resolution?

If the diffraction spikes don’t show up then there really isn’t much point to including the struts in the first place.

[rmandelb]

I'm making a new branch #556 for work on this.

[rmandelb]

@aaronroodman - do you have your own independent calculation of what the PSF should look like with these struts, that I could use as a check on our calculation in GalSim? The issue I'm having is that I'm used to looking at PSFs with filled in struts, but I don't have great intuition for what the PSF should look like with these spiders. So yes, there are diffraction spikes (I confirmed convergence of the calculation in terms of sampling) but I would like a way to test that what I'm getting is right.

In case you have time to quickly eyeball some images, I made a few images on a common logarithmic color map (i.e., color scale is the same in both plots, and goes from the peak pixel value in the first image down to 10^-5 of that value). For both, the pixel scale is 0.02 arcsec, and the image has been convolved with a top hat of that size. So these are images of the LSST PSF at the field center, without atmosphere or aberrations, with teeny pixels.

First, with LSST aperture size and obscuration, and four filled-in struts (with about the right size to be the LSST spider if you filled in the regions between):

Second, with the FITS image you gave me of the pupil plane, so including the spider properly in principle (assuming that all the sampling is done right etc.):

We can see that the diffraction features change in two ways: (1) In the first case, there is some area around the diffraction spikes that is sort of blocked out, whereas in the second case, there is not.

(2) The location of the bumps in the diffraction spikes changes.

You can see some low-level numerical issues where the results go slightly negative and that gives the little boxes out near the array edge in the second case. If I increase the level of oversampling in GalSim, the resulting PSF does not change, just this low-level numerical noise moves around. So it seems like the calculation has converged.

Finally, use of the pupil plane image changes the image rendering time (compared to the filled-in struts case) by a factor of 40. This is without me making any attempt at speeding anything up.

[rmandelb]

... and by "changes" I mean "increases", of course, sigh. Nothing I do ever goes in the other direction. :)

[rmandelb]

Some more investigations while we wait to hear what Aaron thinks of these images:

Building intuition on the effect of strut configurations

In talking with Josh I realized I don't have a great intuition even for what happens if you have naive rectangular struts (i.e. not a proper spider) and you change the width, and that maybe learning about this would help me understand the results with struts vs. spider. So I made the version with filled in struts of widths varying from 2% to 14% of the pupil diameter. The one with narrow struts gives a PSF that looks like this: The diffraction spikes are sort of approaching an unbroken line. With wide struts we get this: This has the bump-like features in the diffraction spikes, and the feature we saw above with struts that there is some modification of the obscured Airy features around the diffraction spikes.

So I think this kind of makes the above plots make more sense.

[rmandelb]

Understanding timing

These calculations with naive struts take ~1.5 seconds to initialize and render, while with the pupil plane image, it takes ~60 seconds (hence my "factor of 40" comment). I tried to understand where this was coming from, and established the following:

• The maxK and stepK for the two cases are quite similar. So it's not like one case is handling much higher resolution images internally.
• The actual loading of the very large pupil plane file takes ~1.5 sec, so it's not like the 60 seconds is all loading images somehow.
• The ~60 seconds is nearly all in the initialization, not in the rendering. So it's some part of the initialization but it's not the "loading the file" part that is dominating.
• For reference, the stepK and maxK for the case where we use the pupil plane image are (1.7 and 506), while a Kolmogorov with FWHM=0.65 arcsec has maxK~20. So basically, it takes a lot of time to build a very high-res PSF and then when convolving with a typical ground-based PSF it's going to chuck out a lot of the information it has gone to great pains to calculate.

[rmandelb]

Thanks to some hints from Josh, I found that 45 of the 60 seconds come from the calculate_maxk step when it initializes the final InterpolatedImage. This would suggest allowing users to pass in a maxk value in cases where they know that they are going to convolve with something else that has a smaller maxk than the optical PSF might have.

So I tried that. Good news: the time to initialize shrank to ~15 seconds as expected. (Next job: figure out what that 15 seconds is.) Bad news: the time to draw the convolution of the optical and atmospheric PSF inflated from well below a second to 40 seconds. That's because the FFT of the InterpolatedImage is only done once, either when you calculate maxK or otherwise, so it has to happen sometime. This plan doesn't work.

Plan B is that we should downsample the file that Aaron sent; perhaps we can get nearly the same results without having ~16 pixels across each spider. I tried downsampling from 2048^2 to 512^2 and even 256^2, in the spirit of living dangerously. 256^2 raised a warning that it wasn't sampled enough even for a simple unobscured Airy, so forget that. 512^2 is nominally fine - at first - and time to initialize drops to 3 seconds from 60 seconds. For reference, the high resolution optical PSF image (2048^2 pupil plane, 0.02 arcsec pixel scale): While for 512^2 it's I was going to suggest that we go with 512^2 as the default, since it clearly gives almost exactly correct diffraction features much more quickly (and then have a "highres" option for those who really care about all the details).

But... when we convolve with a Kolmogorov with FWHM=0.65 arcsec and draw into an image with 0.2 arcsec pixels, something funny happens. After using 2048^2 pupil plane image: After using 512^2 pupil plane image: You see a funny box feature. The origin of the feature is that when we use the lower res pupil plane image, that sets a maximum physical extent of the interpolated image stored inside the OpticalPSF - 3 arcsec or so. But if you are convolving with a reasonable atmospheric PSF, the natural image size is larger than that. So you get a funny box feature that is not easily removed by changing the OpticalPSF kwargs like gsparams, oversampling, or pad factor.

Then I came up with plan C: If the high-res optical PSF is very expensive to initialize but not nearly as expensive to use, we could cache the results of the initialization step in share/ and then the LSST PSF routine could use the cached things (unless different oversampling etc. are requested).

That also doesn't work for two reasons: 1 - the stuff to be cached would also depend on the aberrations, which presumably will differ for different calls 2 - Josh pointed out that the stuff to be cached depends on wavelength

I'm stumped. Anyone else have suggestions?

[rmandelb]

Note, comment updated on github, please read there rather than by e-mail.

[aaronroodman]

Rachel, what wavelength is this done at?

are the plots you show smoothed or interpolated, or are these the raw outputs of the FFTs?

I’m just going to compare against my own Fraunhofer code, especially for time.

thanks, Aaron

[rmandelb]

Aaron - sorry, should have said, lambda=500 nm in all plots.

These are not the raw FFT outputs because the pixel scale that I adopted for drawing the images did not correspond to the natural pixel scale that you get based on the pupil plane image. And I convolved with the pixel response for the pixel scale I was using.

Can you tell me what you'd prefer to do with respect to pixel scale and pixel response for this comparison, and I can remake the images that way?

[rmandelb]

To be clear, by "natural pixel scale" I mean this:

• GalSim takes the input pupil plane images and pads them by a factor of 3, so the diameter becomes 3*aperture.
• The natural pixel scale then becomes (wavelength)/(3*aperture) which for 500 nm is 0.004112 arcsec.

If we weren't padding then the natural scale would just be wavelength/aperture but then there are some potential numerical issues. Anyway, just wanted to clarify what I meant by "natural pixel scale".

[rmjarvis]

Propose closing this, as all of the work related to this is now in imSim, not GalSim.

[rmandelb]

I agree.