Source code for gavo.base.coords

"""
(Mostly deprecated) code to handle coordinate systems and transform
between them.

Basically all of this should be taken over by stc and astropysics.
"""

#c Copyright 2008-2023, the GAVO project <gavo@ari.uni-heidelberg.de>
#c
#c This program is free software, covered by the GNU GPL.  See the
#c COPYING file in the source distribution.


import math
import types
from math import sin, cos, pi #noflake: exported names
import re

import numpy
from numpy import linalg

from gavo import utils
from gavo.utils import DEG
from gavo.utils import pgsphere
from gavo.utils import pyfits

# temporary measure to shut up astropy's configuration parser.
import builtins
builtins._ASTROPY_SETUP_ = True
wcs = utils.DeferredImport("wcs", "from astropy import wcs")


fitsKwPat = re.compile("[A-Z0-9_-]{1,8}$")

[docs]def makePyfitsFromDict(d): """returns a pyfits header with the cards of d.items(). Only keys "looking like" FITS header keywords are used, i.e. all-uppercase and shorter than 9 characters. """ res = pyfits.Header() res.update(dict((str(key), val) for key, val in d.items() if fitsKwPat.match(key) and val is not None)) return res
_wcsTestDict = { "CRVAL1": 0, "CRVAL2": 0, "CRPIX1": 50, "CRPIX2": 50, "CD1_1": 0.01, "CD1_2": 0, "CD2_1": 0, "CD2_2": 0.01, "NAXIS1": 100, "NAXIS2": 100, "CUNIT1": "deg", "CUNIT2": "deg", "CTYPE1": 'RA---TAN-SIP', "CTYPE2": 'DEC--TAN-SIP', "LONPOLE": 180., }
[docs]class Box(object): """is a 2D box. The can be constructed either with two tuples, giving two points delimiting the box, or with four arguments x0, x1, y0, y1. To access the thing, you can either access the x[01], y[01] attributes or use getitem to retrieve the upper right and lower left corner. The slightly silly ordering of the bounding points (larger values first) is for consistency with Postgresql. """ def __init__(self, x0, x1, y0=None, y1=None): if y0 is None: x0, y0 = x0 x1, y1 = x1 lowerLeft = (min(x0, x1), min(y0, y1)) upperRight = (max(x0, x1), max(y0, y1)) self.x0, self.y0 = upperRight self.x1, self.y1 = lowerLeft def __getitem__(self, index): if index==0 or index==-2: return (self.x0, self.y0) elif index==1 or index==-1: return (self.x1, self.y1) else: raise IndexError("len(box) is always 2") def __str__(self): return "((%.4g,%.4g), (%.4g,%.4g))"%(self.x0, self.y0, self.x1, self.y1) def __repr__(self): return "Box((%g,%g), (%g,%g))"%(self.x0, self.y0, self.x1, self.y1)
[docs]def getBbox(points): """returns a bounding box for the sequence of 2-sequences points. The thing returned is a coords.Box. >>> getBbox([(0.25, 1), (-3.75, 1), (-2, 4)]) Box((0.25,4), (-3.75,1)) """ xCoos, yCoos = [[p[i] for p in points] for i in range(2)] return Box(min(xCoos), max(xCoos), min(yCoos), max(yCoos))
[docs]def clampAlpha(alpha): while alpha>360: alpha -= 360 while alpha<0: alpha += 360 return alpha
[docs]def clampDelta(delta): return max(-90, min(90, delta))
[docs]def straddlesStitchingLine(minRA, maxRA): """returns true if something bordered by minRA and maxRA presumably straddles the stitching line. This assumes minRA<maxRA, and that "something" is less than 180 degrees in longitude. Angles are in degrees here. """ return maxRA>270 and minRA<90
def _calcFootprintMonkeypatch(self, hdr=None, undistort=True): """returns the coordinates of the four corners of an image. This is for monkeypatching astropy.wcs. pywcs, at least up to 1.11 did really badly when non-spatial coordinates are present. TODO: see if there's astropy code that does this better. This method relies on the _monkey_naxis_lengths attribute left by getWCS to figure out the axis lengths. pywcs' hdr argument is always ignored here. """ naxis1, naxis2 = self._monkey_naxis_lengths corners = [[1,1],[1,naxis2], [naxis1,naxis2], [naxis1, 1]] # we want our polygons to have a counterclockwise winding direction, # and until I think about this a bit more deeply or get a sufficiently # capable library, I just revert my sense if the determinant of the # CD matrix is negative. That's gut feeling, not proven right. if linalg.det(self.pixel_scale_matrix)<0: corners.reverse() if undistort: return self.all_pix2world(corners, 1) else: return self.wcs_pix2world(corners, 1) def _monkeypatchWCS(wcsObj, naxis, wcsFields): """monkeypatches wcs.WCS instances for DaCHS' purposes. """ wcsObj._dachs_header = wcsFields wcsObj.longAxis = naxis[0] if len(naxis)>1: wcsObj.latAxis = naxis[1] wcsObj._monkey_naxis_lengths = [wcsFields.get("NAXIS%d"%i) for i in naxis] wcsObj.calcFootprint = types.MethodType( _calcFootprintMonkeypatch, wcsObj)
[docs]def getWCS(wcsFields, naxis=(1,2), relax=True): """returns a WCS instance from wcsFields wcsFields can be either a dictionary or a pyfits header giving some kind of WCS information, or an wcs.WCS instance that is returned verbatim. This will return None if no (usable) WCS information is found in the header. We monkeypatch the resulting WCS structure quite a bit. Among others: * there's longAxis and latAxis attributes taken from naxis * there's _dachs_header, containing the incoming k-v pairs * there's _monkey_naxis_length, the lengths along the WCS axes. """ if isinstance(wcsFields, wcs.WCS): return wcsFields if isinstance(wcsFields, dict): wcsFields = makePyfitsFromDict(wcsFields) # pywcs used to invent identity transforms if no WCS keys are present. # Hence. we do some sanity checking up front to weed those out. if ("CD1_1" not in wcsFields and "CDELT1" not in wcsFields and "PC1_1" not in wcsFields and "CNPIX1" not in wcsFields): return None # workaround for a bug in pywcs 1.11: .*_ORDER=0 must not happen for key in ["AP_ORDER", "BP_ORDER", "A_ORDER", "B_ORDER"]: if wcsFields.get(key)==0: del wcsFields[key] # wcslib (at least in 5.16) wants to do DSS corrections if they see # either XPIXELSZ and YPIXELSZ; hide them, because that almost # never works. for key in ["XPIXELSZ", "YPIXELSZ"]: if key in wcsFields: del wcsFields[key] wcsObj = wcs.WCS(wcsFields, relax=relax, naxis=naxis) _monkeypatchWCS(wcsObj, naxis, wcsFields) return wcsObj
[docs]def pix2foc(wcsFields, pixels): """returns the focal plane coordinates for the 2-sequence pixels. (this is a thin wrapper intended to abstract for pix2world's funky calling convention; also, we fix on the silly "0 pixel is 1 convention") """ wcsObj = getWCS(wcsFields) val = wcsObj.pix2foc((pixels[0],), (pixels[1],), 1) return val[0][0], val[1][0]
[docs]def pix2sky(wcsFields, pixels): """returns the sky coordinates for the 2-sequence pixels. (this is a thin wrapper intended to abstract for pix2world's funky calling convention; also, we fix on the silly "0 pixel is 1 convention") """ wcsObj = getWCS(wcsFields) val = wcsObj.all_pix2world((pixels[0],), (pixels[1],), 1) return val[0][0], val[1][0]
[docs]def sky2pix(wcsFields, longLat): """returns the pixel coordinates for the 2-sequence longLad. (this is a thin wrapper intended to abstract for world2pix's funky calling convention; also, we fix on the silly "0 pixel is 1 convention") """ val = getWCS(wcsFields).wcs_world2pix((longLat[0],), (longLat[1],), 1) return val[0][0], val[1][0]
[docs]def getPixelSizeDeg(wcsFields): """returns the sizes of a pixel at roughly the center of the image for wcsFields. Near the pole, this gets a bit weird; we do some limitation of the width of RA pixels there. """ wcs = getWCS(wcsFields) width, height = wcs._dachs_header["NAXIS1"], wcs._dachs_header["NAXIS2"] cosDelta = max(0.01, math.cos(pix2sky(wcs, (width/2, height/2))[1]*DEG)) p0 = pix2sky(wcs, (width/2, height/2)) p1 = pix2sky(wcs, (width/2+1, height/2)) p2 = pix2sky(wcs, (width/2, height/2+1)) return abs(p1[0]-p0[0])*cosDelta, abs(p2[1]-p0[1])
[docs]def getWCSTrafo(wcsFields): """returns a callable transforming pixel to physical coordinates. wcsFields is passed to getWCS, see there for legal types. """ wcs = getWCS(wcsFields) return lambda x, y: pix2sky(wcs, (x, y))
[docs]def getInvWCSTrafo(wcsFields): """returns a callable transforming physical to pixel coordinates. wcsFields is passed to getWCS, see there for legal types. """ wcs = getWCS(wcsFields) return lambda ra, dec: sky2pix(wcs, (ra,dec))
[docs]def getSpolyFromWCSFields(wcsFields): """returns a pgsphere spoly corresponding to wcsFields wcsFields is passed to getWCS, see there for legal types. The polygon returned is computed by using the four corner points assuming a rectangular image. This typically is only loosely related to a proper spherical polygon describing the shape, as image boundaries in the usual projects are not great circles. Also, the spoly will be in the coordinate system of the WCS. If that is not ICRS, you'll probably get something incompatible with most of the VO. """ wcs = getWCS(wcsFields) return pgsphere.SPoly([pgsphere.SPoint.fromDegrees(*p) for p in wcs.calcFootprint(wcs._dachs_header)])
[docs]def getCenterFromWCSFields(wcsFields, spatialAxes=(1,2)): """returns RA and Dec of the center of an image described by wcsFields. This will use the 1-based axes given by spatialAxes to figure out the pixel lengths of the axes. """ wcs = getWCS(wcsFields) if wcs is None: raise ValueError("No WCS solution found") center1 = wcs._dachs_header["NAXIS%s"%spatialAxes[0]]/2. center2 = wcs._dachs_header["NAXIS%s"%spatialAxes[1]]/2. return pix2sky(wcs, (center1, center2))
[docs]def getCoveringCircle(wcsFields, spatialAxes=(1,2)): """returns a pgsphere.scircle large enough to cover the image described by wcsFields. """ wcs = getWCS(wcsFields) center = getCenterFromWCSFields(wcs) height, width = (wcs._dachs_header["NAXIS%s"%spatialAxes[0]], wcs._dachs_header["NAXIS%s"%spatialAxes[1]]) radius = max( getGCDist(center, pix2sky(wcs, corner)) for corner in [(0, 0), (height, 0), (0, width), (height, width)]) return pgsphere.SCircle(pgsphere.SPoint.fromDegrees(*center), radius*DEG)
[docs]def getSkyWCS(hdr): """returns a pair of a wcs.WCS instance and a sequence of the spatial axes. This will be None, () if no WCS could be discerned. There's some heuristics involved in the identification of the spatial coordinates that will probably fail for unconventional datasets. """ wcsAxes = [] # heuristics: iterate through CTYPEn, anything that's got # a - is supposed to be a position (needs some refinement :-) for ind in range(1, hdr["NAXIS"]+1): if "-" in hdr.get("CTYPE%s"%ind, ""): wcsAxes.append(ind) if not wcsAxes: # more heuristics to be inserted here return None, () if len(wcsAxes)!=2: raise utils.ValidationError("This FITS has !=2" " spatial WCS axes. Please contact the DaCHS authors and" " make them support it.", "PUBDID") return getWCS(hdr, naxis=wcsAxes), wcsAxes
[docs]def getPixelLimits(cooPairs, wcsFields): """returns pixel cutout slices for covering cooPairs in an image with wcsFields. cooPairs is a sequence of (ra, dec) tuples. wcsFields is a DaCHS-enhanced wcs.WCS instance. Behaviour if cooPairs use a different coordinate system from wcsFields is undefined at this point. Each cutout slice is a tuple of (FITS axis number, lower limit, upper limit). If cooPairs is off the wcsFields coverage, a null cutout on the longAxis is returned. """ latAxis = wcsFields.latAxis longAxis = wcsFields.longAxis latPixels = wcsFields._dachs_header["NAXIS%d"%latAxis] longPixels = wcsFields._dachs_header["NAXIS%d"%longAxis] # pywcs used to do really funny things when we "wrap around". Therefore, we # clamp values to be within the pywcs-safe region. Note that # due to spherical magic this is not enough to ensure +/-Inf behaviour # for SODA; TODO: re-evaluate this for astropy.wcs at some point. cooPairs = [( min(359.99999, max(-89.9999, ra)), min(89.9999, max(-89.9999, dec))) for ra, dec in cooPairs] slices = [] pixelFootprint = numpy.asarray( numpy.round(wcsFields.wcs_world2pix(cooPairs, 1)), numpy.int32) pixelLimits = [ [min(pixelFootprint[:,0]), max(pixelFootprint[:,0])], [min(pixelFootprint[:,1]), max(pixelFootprint[:,1])]] # see if we're completely off coverage if pixelLimits[0][1]<0 or pixelLimits[1][1]<0: return [[longAxis, 0, 0]] if pixelLimits[0][0]>longPixels or pixelLimits[1][0]>latPixels: return [[longAxis, 0, 0]] # now crop to the actual pixel values pixelLimits = [ [max(pixelLimits[0][0], 1), min(pixelLimits[0][1], longPixels)], [max(pixelLimits[1][0], 1), min(pixelLimits[1][1], latPixels)]] if pixelLimits[0]!=[1, longPixels]: slices.append([longAxis]+pixelLimits[0]) if pixelLimits[1]!=[1, latPixels]: slices.append([latAxis]+pixelLimits[1]) return slices
# let's do a tiny vector type. It's really not worth getting some dependency # for this.
[docs]class Vector3(object): """is a 3d vector that responds to both .x... and [0]... >>> x, y = Vector3(1,2,3), Vector3(2,3,4) >>> x+y Vector3(3.000000,5.000000,7.000000) >>> 4*x Vector3(4.000000,8.000000,12.000000) >>> x*4 Vector3(4.000000,8.000000,12.000000) >>> x*y 20 >>> "%.6f"%abs(x) '3.741657' >>> "%.4f"%abs((x+y).normalized()) '1.0000' """ def __init__(self, x, y=None, z=None): if isinstance(x, tuple): self.coos = x else: self.coos = (x, y, z) def __repr__(self): return "Vector3(%f,%f,%f)"%tuple(self.coos) def __str__(self): def cutoff(c): if abs(c)<1e-10: return 0 else: return c rounded = [cutoff(c) for c in self.coos] return "[%.2g,%.2g,%.2g]"%tuple(rounded) def __getitem__(self, index): return self.coos[index] def __mul__(self, other): """does either scalar multiplication if other is not a Vector3, or a scalar product. """ if isinstance(other, Vector3): return self.x*other.x+self.y*other.y+self.z*other.z else: return Vector3(self.x*other, self.y*other, self.z*other) __rmul__ = __mul__ def __truediv__(self, scalar): return Vector3(self.x/scalar, self.y/scalar, self.z/scalar) def __add__(self, other): return Vector3(self.x+other.x, self.y+other.y, self.z+other.z) def __sub__(self, other): return Vector3(self.x-other.x, self.y-other.y, self.z-other.z) def __abs__(self): return math.sqrt(self.x**2+self.y**2+self.z**2)
[docs] def cross(self, other): return Vector3(self.y*other.z-self.z*other.y, self.z*other.x-self.x*other.z, self.x*other.y-self.y*other.x)
[docs] def normalized(self): return self/abs(self)
[docs] def getx(self): return self.coos[0]
[docs] def setx(self, x): self.coos[0] = x
x = property(getx, setx)
[docs] def gety(self): return self.coos[1]
[docs] def sety(self, y): self.coos[1] = y
y = property(gety, sety)
[docs] def getz(self): return self.coos[2]
[docs] def setz(self, z): self.coos[2] = z
z = property(getz, setz)
[docs]def sgn(a): if a<0: return -1 elif a>0: return 1 else: return 0
[docs]def computeUnitSphereCoords(alpha, delta): # TODO: replaced by mathtricks.spherToCart """returns the 3d coordinates of the intersection of the direction vector given by the spherical coordinates alpha and delta with the unit sphere. alpha and delta are given in degrees. >>> print(computeUnitSphereCoords(0,0)) [1,0,0] >>> print(computeUnitSphereCoords(0, 90)) [0,0,1] >>> print(computeUnitSphereCoords(90, 90)) [0,0,1] >>> print(computeUnitSphereCoords(90, 0)) [0,1,0] >>> print(computeUnitSphereCoords(180, -45)) [-0.71,0,-0.71] """ return Vector3(*utils.spherToCart(alpha*DEG, delta*DEG))
[docs]def dirVecToCelCoos(dirVec): """returns alpha, delta in degrees for the direction vector dirVec. >>> dirVecToCelCoos(computeUnitSphereCoords(25.25, 12.125)) (25.25, 12.125) >>> dirVecToCelCoos(computeUnitSphereCoords(25.25, 12.125)*16) (25.25, 12.125) >>> "%g,%g"%dirVecToCelCoos(computeUnitSphereCoords(25.25, 12.125)+ ... computeUnitSphereCoords(30.75, 20.0)) '27.9455,16.0801' """ dirVec = dirVec.normalized() alpha = math.atan2(dirVec.y, dirVec.x) if alpha<0: alpha += 2*math.pi return alpha*180./math.pi, math.asin(dirVec.z)*180./math.pi
[docs]def getTangentialUnits(cPos): """returns the unit vectors for RA and Dec at the unit circle position cPos. We compute them by solving u_1*p_1+u_2*p_2=0 (we already know that u_3=0) simultaneously with u_1^2+u_2^2=1 for RA, and by computing the cross product of the RA unit and the radius vector for dec. This becomes degenerate at the poles. If we're exactly on a pole, we *define* the unit vectors as (1,0,0) and (0,1,0). Orientation is a pain -- the convention used here is that unit delta always points to the pole. >>> cPos = computeUnitSphereCoords(45, -45) >>> ua, ud = getTangentialUnits(cPos) >>> print(abs(ua), abs(ud), cPos*ua, cPos*ud) 1.0 1.0 0.0 0.0 >>> print(ua, ud) [-0.71,0.71,0] [-0.5,-0.5,-0.71] >>> ua, ud = getTangentialUnits(computeUnitSphereCoords(180, 60)) >>> print(ua, ud) [0,-1,0] [0.87,0,0.5] >>> ua, ud = getTangentialUnits(computeUnitSphereCoords(0, 60)) >>> print(ua, ud) [0,1,0] [-0.87,0,0.5] >>> ua, ud = getTangentialUnits(computeUnitSphereCoords(0, -60)) >>> print(ua, ud) [0,1,0] [-0.87,0,-0.5] """ try: normalizer = 1/math.sqrt(cPos.x**2+cPos.y**2) except ZeroDivisionError: return Vector3(1,0,0), Vector3(0,1,0) alphaUnit = normalizer*Vector3(cPos.y, -cPos.x, 0) deltaUnit = normalizer*Vector3(cPos.x*cPos.z, cPos.y*cPos.z, -cPos.x**2-cPos.y**2) # now orient the vectors: in delta, we always look towards the pole if sgn(cPos.z)!=sgn(deltaUnit.z): deltaUnit = -1*deltaUnit # XXX this breaks on the equator # The orientation of alphaUnit depends on the hemisphere if cPos.z<0: # south if deltaUnit.cross(alphaUnit)*cPos<0: alphaUnit = -1*alphaUnit else: # north if deltaUnit.cross(alphaUnit)*cPos>0: alphaUnit = -1*alphaUnit return alphaUnit, deltaUnit
[docs]def movePm(alphaDeg, deltaDeg, pmAlpha, pmDelta, timeDiff, foreshort=0): """returns alpha and delta for an object with pm pmAlpha after timeDiff. pmAlpha has to have cos(delta) applied, everything is supposed to be in degrees, the time unit is yours to choose. """ alpha, delta = alphaDeg/180.*math.pi, deltaDeg/180.*math.pi pmAlpha, pmDelta = pmAlpha/180.*math.pi, pmDelta/180.*math.pi sd, cd = math.sin(delta), math.cos(delta) sa, ca = math.sin(alpha), math.cos(alpha) muAbs = math.sqrt(pmAlpha**2+pmDelta**2); muTot = muAbs+0.5*foreshort*timeDiff; if muAbs<1e-20: return alphaDeg, deltaDeg # this is according to Mueller, 115 (4.94) dirA = pmAlpha/muAbs; dirD = pmDelta/muAbs; sinMot = sin(muTot*timeDiff); cosMot = cos(muTot*timeDiff); dirVec = Vector3(-sd*ca*dirD*sinMot - sa*dirA*sinMot + cd*ca*cosMot, -sd*sa*dirD*sinMot + ca*dirA*sinMot + cd*sa*cosMot, +cd*dirD*sinMot + sd*cosMot) return dirVecToCelCoos(dirVec)
[docs]def getGCDist(pos1, pos2): """returns the distance along a great circle between two points. The distance is in degrees, the input positions are in degrees. """ scalarprod = computeUnitSphereCoords(*pos1)*computeUnitSphereCoords(*pos2) # cope with numerical trouble if scalarprod>=1: return 0 return math.acos(scalarprod)/DEG
def _test(): # pragma: no cover import doctest doctest.testmod() if __name__=="__main__": # pragma: no cover _test()