@jbosch Sorry for not replying earlier! Essentially I want to use Galsim to simulate galaxies and call the PSF model at their coordinates to use it in correcting for their elipticities etc. So the idea was to make up a single PSF model at all positions by combining all the CCDs from one visit, or is there a better way to approach this?
I’ve been trying to use the wcs.getPixelOrigin() (which seemed to be the only attribute that was not either 0 or the same for every CCD) to do this. But it’s returning the same/overlapping positions for some CCDs in the same visit, which is not what I expected. Am I missing something/is there a better attribute to look for?
@husni, I’m afraid I still don’t quite understand - I get that you want to simulate galaxies convolved with the (modeled) HSC PSF, but I’m not sure what you mean by a “single PSF model at all positions”. If you want a centered image of the PSF model at a position on a CCD image, you can use Psf.computeKernelImage(afw.geom.Point2D(x, y)), with (x, y) given in CCD pixel coordinates.
Combining all of the PSF models for one visit somehow doesn’t seem like a meaningful thing to do, unless you just want some kind of approximate average PSF for the visit. If that’s the case, I’d just recommend just calling the above at a number of randomly-drawn points from randomly-selected CCDs in the visit and averaging them. Because the PSF models are defined and evaluated in pixel coordinates, you shouldn’t need a WCS to evaluate them at all.
If you want to resample the PSF images to sky coordinates, it’s a little trickier (I’d recommend using our WarpedPsf class, which transforms a PSF through a pair of WCSs), but I can walk you through that too.
[quote=“jbosch, post:42, topic:964”] @husni, I’m afraid I still don’t quite understand - I get that you want to simulate galaxies convolved with the (modeled) HSC PSF, but I’m not sure what you mean by a “single PSF model at all positions”. [/quote]
@jbosch Sorry that was not a very clear way to explain it. By that I meant a PSF model that covers the entire part of sky that the HSC covers.
This is actually exactly what I would like to do, only I am trying to construct the Psf. Currently, I’m getting it from the calexp, which I’m in turn getting it from the butler for a single visit and single ccd and looping over them, so it only contains limited x,y positions from one CCD in every iteration. This is why I’m trying to do this addition to get a Psf object that spans the entire sky. Or is there a way to ask the butler for the calexp of the entire sky?
Ok, I understand. I think the core problem here is that a particular PSF is valid only for a particular image, and there is no single image that covers the entire area of sky that HSC covers. There are multiple overlapping epoch-level images at any given point on the sky, of course, so even at a single point no one of those is “the PSF”.
However, when we make coadd images, we coadd the PSFs in a way that is consistent with the way we built the coadd (i.e. it’s the appropriate model for the effective PSF of the coadd). If you get the PSF object from a deepCoadd_calexp image, that’s what you’re getting - but even those are only valid over a small region of sky, called a patch. All of the patches in a larger region called a tract share a common pixel grid (they’re conceptually square subimages of a larger virtual image for the full tract).
You could conceivably make a CoaddPsf that combines PSFs over a tract, but it would be very large and cumbersome to work with, because it would essentially have to store all of the PSF models and WCSs of all of the input images that went into all of those coadd images. I’m afraid it’s not possible to make a CoaddPsf that covers multiple tracts, because the class assumes there’s one WCS that describes the pixel grid of the image to which it corresponds, and different tracts have different WCSs.
If you want one object that can return the PSF at an arbitrary point in the survey, then, I suspect you’ll be better off writing a new class that can look up the tract and patch that correspond to a particular point in celestial coordinates, load the CoaddPsf for that patch, and then evaluate the PSF at that point after converting it to tract coordinates. If you load the deepCoadd_skyMap dataset (this will be an instance of lsst.skymap.RingsSkyMap for HSC), you’ll get an object that can do those tract/patch lookups and coordinate transformatoins.
@jbosch after discussing with @rmandelb I think I’ve made a terminology mistake which has caused the confusion. I should have said the entire field of view instead of the entire sky, since I meant the PSF model from a single visit, rather than deep-coadding all the visits to get the effective PSF. Sorry about that.
So the problem is just that when I try to e.g. plot the positions of the pixels from each of the ccds in a single visit on the same plot I’m finding significant overlap between the CCDs (which if I understand correctly shouldn’t happen within the same visit?). I’m getting the position of each pixel using wcs.getPixelOrigin()[0] + x, and wcs.getPixelOrigin()[1] + y, where wcs = calexp.getWcs(), x and y are between 0 and calexp.getWidth() and calexp.getHeight() respectively, and calexp is being outputted from the butler for each CCD. Is this the right thing to do/is there a better way to do this?
In order to put all of the CCDs in the focal plane on a single plot, you can either make the plot in focal plane coordinates and convert to CCD coordinates using our “cameraGeom” library, or, for a particular visit, you can make the plot in sky coordinates and use the Wcs.
If you want to use the Wcs, you should just use skyToPixel and pixelToSky to transform between sky coordinates and pixel coordinates. Using getPixelOrigin the way you have is only valid when all of the CCDs have the same Wcs parameters aside from their origin parameters. That can never be precisely true for a quality Wcs, because there are at least small differences in rotation and perhaps scale between CCDs (and for HSC, fairly significant optical distortion as well).
