DM-26485: Add skyToPixelArray and pixelToSkyArray convenience functions. #538
Conversation
There was a problem hiding this comment.
This is convenient, and appears to be faster than using the pixelToSky and skyToPixel for more than 50-ish points. I'm wondering if this should allow the input arrays to be 2-D as well. I could imagine getting a set of (x,y) points from meshgrid and wanting to convert them to (RA,DEC). Either, this function could reshape the input arrays (if 2D) internally and return a pair of arrays of the same shape, or explicitly mention in the docstring that the input arrays have to be one dimensional.
| y : `np.ndarray` | ||
| Array of y points. | ||
| degrees : `bool`, optional | ||
| Return ra/dec arrays in degrees if True. |
There was a problem hiding this comment.
Default value should be specified here
There was a problem hiding this comment.
So I used to specify the default, but got barked at by somebody because it is redundant because the default is in the definition...
There was a problem hiding this comment.
Yes, sphinx should include the default value in the output so including it here risks getting out of sync.
There was a problem hiding this comment.
Ah, if it reads out from the function signature, then that should suffice.
| Array of Declination. Units are radians unless | ||
| degrees=True. | ||
| degrees : `bool`, optional | ||
| Input ra/dec arrays are degrees if True. |
There was a problem hiding this comment.
Default value should be specified here as well
| Array of y points. | ||
| """ | ||
| if degrees: | ||
| rd = np.deg2rad(np.vstack((ra, dec))) |
There was a problem hiding this comment.
I think it would make sense to first denote the output from vstack as rd, and then convert it to degrees if degrees=True, rather than having an else statement.
| Parameters | ||
| ---------- | ||
| x : `np.ndarray` | ||
| Array of x points. |
There was a problem hiding this comment.
Here and elsewhere, Array of x values instead of points?
| class SkyWcs: | ||
| def pixelToSkyArray(self, x, y, degrees=False): | ||
| """ | ||
| Convert numpy array pixels (x/y) to numpy array sky (ra/dec) |
There was a problem hiding this comment.
It'd be clearer to use (x,y) instead of (x/y) and similarly for (ra/dec). Similarly for the inverse function.
| rd[0, :] = rd[0, :] % (2.*np.pi) | ||
|
|
||
| if degrees: | ||
| return np.rad2deg(rd[0, :]), np.rad2deg(rd[1, :]) |
There was a problem hiding this comment.
I'm not sure. but it may be more efficient to call rad2deg with rd first and then split them, and to me that seems a bit more elegant.
There was a problem hiding this comment.
Also, since rd[0,:] and rd[1,:] are used quite a few times, it may be better to name them ra and dec. Alternatively, you can use np.vsplit in the previous line and not define rd at all and ignore the comment above.
|
|
||
| xy = self.getTransform().getMapping().applyInverse(rd) | ||
|
|
||
| return xy[0, :], xy[1, :] |
There was a problem hiding this comment.
This does the job, but I'm wondering if using np.vsplit(xy,2) would be more natural in reversing the action of np.vstack earlier.
There was a problem hiding this comment.
I did not know about vsplit. That is indeed the perfect function right here!
|
|
||
| pixPoints = wcs.skyToPixel(spherePointList) | ||
|
|
||
| x, y = wcs.skyToPixelArray(np.deg2rad(raPoints), np.deg2rad(decPoints)) |
There was a problem hiding this comment.
Could you explicitly set degrees=False, future-proofing a change in the default value?
There was a problem hiding this comment.
This is an interesting question. I see the point of explicitly setting this, but at the same time I think I'd like the test to fail if somebody inadvertently changes the default behavior and then they can explicitly update the tests if desired? E.g., it can be argued that tests should not be "future proofed" in certain ways.
There was a problem hiding this comment.
Hmm, interesting. If someone, as a developer, changed the default behavior intentionally, they'd also fix the unit tests, but someone else, as a user, will find their code using these routines breaking. But if this is to protect against unintentional change of default behavior, I see your point. I'm okay with this staying as it is then.
|
Regarding a grid of points, AST already has a super optimized function for returning all the transformed WCS for a grid of points. You wouldn't want to use the implementation of this routine to do a regular grid. |
|
It's implemented in astshim as |
|
Exactly. And the documentation is somewhat sparse and it is tricky to use, but I did figure it out: https://github.com/LSSTDESC/supreme/blob/75a2008d203e28b3d42738374823f1f5e42d9db1/supreme/utils.py#L203 I don't know if it's worth porting that over as well? But maybe on a separate ticket. |
erykoff
force-pushed
the
tickets/DM-26485
branch
from
be33dbe
to
7e4ea93
Compare
| ra %= (2.*np.pi) | ||
|
|
||
| if degrees: | ||
| return np.rad2deg(ra.ravel()), np.rad2deg(dec.ravel()) |
There was a problem hiding this comment.
I don't think ravel is useful here. If xy does not have exactly two rows, the call to applyForward will throw an error. And when xy has exactly two rows, ra and dec will be flat. Flattening isn't incorrect, but may just add to unnecessary calls without actually covering any edge cases.
There was a problem hiding this comment.
Unfortunately, it seems that vsplit does not return 1d arrays but returns 2 2d arrays with shape [1, npoints]. The ravel call will flatten that without copying memory (unlike flatten which does a copy). In addition, a ravel on a true 1d array is essentially a no-op.
No description provided.