From 2d062c20696866d65b554c42d4b693d7a201fbba Mon Sep 17 00:00:00 2001 From: Matt Hilton Date: Tue, 26 May 2026 21:55:59 +0200 Subject: [PATCH 1/5] Convert all docstrings from epydoc to Google style Rewrites @type/@param/@rtype/@return tags to Args:/Returns: sections, @note: Raises ... to Raises: sections, general @note: to Note: sections, and strips L{...} cross-reference syntax throughout all seven astLib modules. Co-Authored-By: Claude Sonnet 4.6 --- astLib/astCalc.py | 222 ++++++----- astLib/astCoords.py | 142 ++++--- astLib/astImages.py | 812 +++++++++++++++++++-------------------- astLib/astPlots.py | 548 +++++++++++++------------- astLib/astSED.py | 910 +++++++++++++++++++++----------------------- astLib/astStats.py | 573 ++++++++++++++-------------- astLib/astWCS.py | 168 ++++---- 7 files changed, 1658 insertions(+), 1717 deletions(-) diff --git a/astLib/astCalc.py b/astLib/astCalc.py index 89e2970..285e84d 100755 --- a/astLib/astCalc.py +++ b/astLib/astCalc.py @@ -17,7 +17,7 @@ """The dark energy density (in the form of a cosmological constant) at z=0.""" OMEGA_R0 = 8.24E-5 -"""The radiation density at z=0 (note this is only used currently in calculation of L{Ez}).""" +"""The radiation density at z=0 (note this is only used currently in calculation of Ez).""" H0 = 70.0 """The Hubble parameter (in km/s/Mpc) at z=0.""" @@ -35,10 +35,11 @@ def dl(z): """Calculates the luminosity distance in Mpc at redshift z. - @type z: float - @param z: redshift - @rtype: float - @return: luminosity distance in Mpc + Args: + z (float): redshift + + Returns: + float: luminosity distance in Mpc """ @@ -51,10 +52,11 @@ def dl(z): def da(z): """Calculates the angular diameter distance in Mpc at redshift z. - @type z: float - @param z: redshift - @rtype: float - @return: angular diameter distance in Mpc + Args: + z (float): redshift + + Returns: + float: angular diameter distance in Mpc """ DM = dm(z) @@ -67,10 +69,11 @@ def dm(z): """Calculates the transverse comoving distance (proper motion distance) in Mpc at redshift z. - @type z: float - @param z: redshift - @rtype: float - @return: transverse comoving distance (proper motion distance) in Mpc + Args: + z (float): redshift + + Returns: + float: transverse comoving distance (proper motion distance) in Mpc """ @@ -101,10 +104,11 @@ def dm(z): def dc(z): """Calculates the line of sight comoving distance in Mpc at redshift z. - @type z: float - @param z: redshift - @rtype: float - @return: transverse comoving distance (proper motion distance) in Mpc + Args: + z (float): redshift + + Returns: + float: transverse comoving distance (proper motion distance) in Mpc """ @@ -129,10 +133,11 @@ def dVcdz(z): """Calculates the line of sight comoving volume element per steradian dV/dz at redshift z. - @type z: float - @param z: redshift - @rtype: float - @return: comoving volume element per steradian + Args: + z (float): redshift + + Returns: + float: comoving volume element per steradian """ @@ -146,10 +151,11 @@ def dl2z(distanceMpc): """Calculates the redshift z corresponding to the luminosity distance given in Mpc. - @type distanceMpc: float - @param distanceMpc: distance in Mpc - @rtype: float - @return: redshift + Args: + distanceMpc (float): distance in Mpc + + Returns: + float: redshift """ @@ -187,10 +193,11 @@ def dc2z(distanceMpc): """Calculates the redshift z corresponding to the comoving distance given in Mpc. - @type distanceMpc: float - @param distanceMpc: distance in Mpc - @rtype: float - @return: redshift + Args: + distanceMpc (float): distance in Mpc + + Returns: + float: redshift """ @@ -228,8 +235,8 @@ def t0(): """Calculates the age of the universe in Gyr at z=0 for the current set of cosmological parameters. - @rtype: float - @return: age of the universe in Gyr at z=0 + Returns: + float: age of the universe in Gyr at z=0 """ @@ -251,13 +258,14 @@ def t0(): #------------------------------------------------------------------------------ def tl(z): - """ Calculates the lookback time in Gyr to redshift z for the current set + """Calculates the lookback time in Gyr to redshift z for the current set of cosmological parameters. - @type z: float - @param z: redshift - @rtype: float - @return: lookback time in Gyr to redshift z + Args: + z (float): redshift + + Returns: + float: lookback time in Gyr to redshift z """ OMEGA_K = 1.0 - OMEGA_M0 - OMEGA_L0 @@ -281,10 +289,11 @@ def tz(z): """Calculates the age of the universe at redshift z for the current set of cosmological parameters. - @type z: float - @param z: redshift - @rtype: float - @return: age of the universe in Gyr at redshift z + Args: + z (float): redshift + + Returns: + float: age of the universe in Gyr at redshift z """ @@ -297,13 +306,15 @@ def tl2z(tlGyr): """Calculates the redshift z corresponding to lookback time tlGyr given in Gyr. - @type tlGyr: float - @param tlGyr: lookback time in Gyr - @rtype: float - @return: redshift - - @note: Raises ValueError if tlGyr is not positive. - + Args: + tlGyr (float): lookback time in Gyr + + Returns: + float: redshift + + Raises: + ValueError: if tlGyr is not positive. + """ if tlGyr < 0.: raise ValueError('Lookback time must be positive') @@ -342,12 +353,14 @@ def tz2z(tzGyr): """Calculates the redshift z corresponding to age of the universe tzGyr given in Gyr. - @type tzGyr: float - @param tzGyr: age of the universe in Gyr - @rtype: float - @return: redshift - - @note: Raises ValueError if Universe age not positive + Args: + tzGyr (float): age of the universe in Gyr + + Returns: + float: redshift + + Raises: + ValueError: if Universe age is not positive. """ if tzGyr <= 0: @@ -362,12 +375,12 @@ def absMag(appMag, distMpc): """Calculates the absolute magnitude of an object at given luminosity distance in Mpc. - @type appMag: float - @param appMag: apparent magnitude of object - @type distMpc: float - @param distMpc: distance to object in Mpc - @rtype: float - @return: absolute magnitude of object + Args: + appMag (float): apparent magnitude of object + distMpc (float): distance to object in Mpc + + Returns: + float: absolute magnitude of object """ absMag = appMag - (5.0*math.log10(distMpc*1.0e5)) @@ -380,10 +393,11 @@ def Ez(z): parameter with redshift, at redshift z for the current set of cosmological parameters. See, e.g., Bryan & Norman 1998 (ApJ, 495, 80). - @type z: float - @param z: redshift - @rtype: float - @return: value of E(z) at redshift z + Args: + z (float): redshift + + Returns: + float: value of E(z) at redshift z """ @@ -397,10 +411,11 @@ def Ez2(z): parameter with redshift, at redshift z for the current set of cosmological parameters. See, e.g., Bryan & Norman 1998 (ApJ, 495, 80). - @type z: float - @param z: redshift - @rtype: float - @return: value of E(z)^2 at redshift z + Args: + z (float): redshift + + Returns: + float: value of E(z)^2 at redshift z """ # This form of E(z) is more reliable at high redshift. It is basically the @@ -419,10 +434,11 @@ def OmegaMz(z): """Calculates the matter density of the universe at redshift z. See, e.g., Bryan & Norman 1998 (ApJ, 495, 80). - @type z: float - @param z: redshift - @rtype: float - @return: matter density of universe at redshift z + Args: + z (float): redshift + + Returns: + float: matter density of universe at redshift z """ ez2 = Ez2(z) @@ -433,12 +449,13 @@ def OmegaMz(z): #------------------------------------------------------------------------------ def OmegaLz(z): - """ Calculates the dark energy density of the universe at redshift z. + """Calculates the dark energy density of the universe at redshift z. - @type z: float - @param z: redshift - @rtype: float - @return: dark energy density of universe at redshift z + Args: + z (float): redshift + + Returns: + float: dark energy density of universe at redshift z """ ez2 = Ez2(z) @@ -447,12 +464,13 @@ def OmegaLz(z): #------------------------------------------------------------------------------ def OmegaRz(z): - """ Calculates the radiation density of the universe at redshift z. + """Calculates the radiation density of the universe at redshift z. + + Args: + z (float): redshift - @type z: float - @param z: redshift - @rtype: float - @return: radiation density of universe at redshift z + Returns: + float: radiation density of universe at redshift z """ ez2 = Ez2(z) @@ -461,18 +479,18 @@ def OmegaRz(z): #------------------------------------------------------------------------------ def DeltaVz(z): - """Calculates the density contrast of a virialised region S{Delta}V(z), - assuming a S{Lambda}CDM-type flat cosmology. See, e.g., Bryan & Norman + """Calculates the density contrast of a virialised region DeltaV(z), + assuming a LambdaCDM-type flat cosmology. See, e.g., Bryan & Norman 1998 (ApJ, 495, 80). - @type z: float - @param z: redshift - @rtype: float - @return: density contrast of a virialised region at redshift z + Args: + z (float): redshift + + Returns: + float: density contrast of a virialised region at redshift z - @note: If OMEGA_M0+OMEGA_L0 is not equal to 1, this routine exits and - prints an error - message to the console. + Raises: + Exception: if OMEGA_M0+OMEGA_L0 is not equal to 1 (non-flat cosmology). """ @@ -491,20 +509,19 @@ def RVirialXRayCluster(kT, z, betaT): """Calculates the virial radius (in Mpc) of a galaxy cluster at redshift z with X-ray temperature kT, assuming self-similar evolution and a flat cosmology. See Arnaud et al. 2002 (A&A, 389, 1) and Bryan & Norman 1998 - (ApJ, 495, 80). A flat S{Lambda}CDM-type flat cosmology is assumed. + (ApJ, 495, 80). A flat LambdaCDM-type flat cosmology is assumed. - @type kT: float - @param kT: cluster X-ray temperature in keV - @type z: float - @param z: redshift - @type betaT: float - @param betaT: the normalisation of the virial relation, for which Evrard et - al. 1996 (ApJ,469, 494) find a value of 1.05 - @rtype: float - @return: virial radius of cluster in Mpc + Args: + kT (float): cluster X-ray temperature in keV + z (float): redshift + betaT (float): the normalisation of the virial relation, for which + Evrard et al. 1996 (ApJ, 469, 494) find a value of 1.05 - @note: If OMEGA_M0+OMEGA_L0 is not equal to 1, this routine exits and - prints an error message to the console. + Returns: + float: virial radius of cluster in Mpc + + Raises: + Exception: if OMEGA_M0+OMEGA_L0 is not equal to 1 (non-flat cosmology). """ @@ -528,4 +545,3 @@ def RVirialXRayCluster(kT, z, betaT): raise Exception("cosmology is NOT flat.") #------------------------------------------------------------------------------ - diff --git a/astLib/astCoords.py b/astLib/astCoords.py index f95383f..d750983 100755 --- a/astLib/astCoords.py +++ b/astLib/astCoords.py @@ -18,12 +18,12 @@ def hms2decimal(RAString, delimiter): """Converts a delimited string of Hours:Minutes:Seconds format into decimal degrees. - @type RAString: string - @param RAString: coordinate string in H:M:S format - @type delimiter: string - @param delimiter: delimiter character in RAString - @rtype: float - @return: coordinate in decimal degrees + Args: + RAString (str): coordinate string in H:M:S format + delimiter (str): delimiter character in RAString + + Returns: + float: coordinate in decimal degrees """ # is it in HH:MM:SS format? @@ -48,12 +48,12 @@ def dms2decimal(decString, delimiter): """Converts a delimited string of Degrees:Minutes:Seconds format into decimal degrees. - @type decString: string - @param decString: coordinate string in D:M:S format - @type delimiter: string - @param delimiter: delimiter character in decString - @rtype: float - @return: coordinate in decimal degrees + Args: + decString (str): coordinate string in D:M:S format + delimiter (str): delimiter character in decString + + Returns: + float: coordinate in decimal degrees """ # is it in DD:MM:SS format? @@ -83,12 +83,12 @@ def decimal2hms(RADeg, delimiter): """Converts decimal degrees to string in Hours:Minutes:Seconds format with user specified delimiter. - @type RADeg: float - @param RADeg: coordinate in decimal degrees - @type delimiter: string - @param delimiter: delimiter character in returned string - @rtype: string - @return: coordinate string in H:M:S format + Args: + RADeg (float): coordinate in decimal degrees + delimiter (str): delimiter character in returned string + + Returns: + str: coordinate string in H:M:S format """ hours = (RADeg/360.0)*24 @@ -137,12 +137,12 @@ def decimal2dms(decDeg, delimiter): """Converts decimal degrees to string in Degrees:Minutes:Seconds format with user specified delimiter. - @type decDeg: float - @param decDeg: coordinate in decimal degrees - @type delimiter: string - @param delimiter: delimiter character in returned string - @rtype: string - @return: coordinate string in D:M:S format + Args: + decDeg (float): coordinate in decimal degrees + delimiter (str): delimiter character in returned string + + Returns: + str: coordinate string in D:M:S format """ # Positive @@ -235,19 +235,18 @@ def calcAngSepDeg(RADeg1, decDeg1, RADeg2, decDeg2): in decimal degrees) in decimal degrees. Note that RADeg2, decDeg2 can be numpy arrays. - @type RADeg1: float - @param RADeg1: R.A. in decimal degrees for position 1 - @type decDeg1: float - @param decDeg1: dec. in decimal degrees for position 1 - @type RADeg2: float or numpy array - @param RADeg2: R.A. in decimal degrees for position 2 - @type decDeg2: float or numpy array - @param decDeg2: dec. in decimal degrees for position 2 - @rtype: float or numpy array, depending upon type of RADeg2, decDeg2 - @return: angular separation in decimal degrees + Args: + RADeg1 (float): R.A. in decimal degrees for position 1 + decDeg1 (float): dec. in decimal degrees for position 1 + RADeg2 (float or numpy.ndarray): R.A. in decimal degrees for position 2 + decDeg2 (float or numpy.ndarray): dec. in decimal degrees for position 2 + + Returns: + float or numpy.ndarray: angular separation in decimal degrees; type + matches the type of RADeg2, decDeg2 """ - + a=numpy.sin(numpy.radians(decDeg1))*numpy.sin(numpy.radians(decDeg2))+numpy.cos(numpy.radians(decDeg1))*numpy.cos(numpy.radians(decDeg2))*numpy.cos(numpy.radians(RADeg1-RADeg2)) mask=numpy.greater(a, 1.0) if mask.sum() > 0: @@ -262,7 +261,7 @@ def calcAngSepDeg(RADeg1, decDeg1, RADeg2, decDeg2): else: a=-1.0 r=numpy.degrees(numpy.arccos(a)) - + # Above gives nan when RADeg1, decDeg1 == RADeg1, decDeg2 indexList=numpy.where(numpy.isnan(numpy.atleast_1d(r)) == True)[0] tolerance=1e-6 @@ -275,13 +274,13 @@ def calcAngSepDeg(RADeg1, decDeg1, RADeg2, decDeg2): else: raise Exception("astCoords: calcAngSepDeg - encountered nan not due to equal RADeg, decDeg coords") elif type(RADeg1) == numpy.ndarray: - if abs(RADeg2 - RADeg1[index]) < tolerance and abs(decDeg2 -decDeg1[index]) < tolerance: + if abs(RADeg2 - RADeg1[index]) < tolerance and abs(decDeg2 -decDeg1[index]) < tolerance: r[index]=0.0 else: raise Exception("astCoords: calcAngSepDeg - encountered nan not due to equal RADeg, decDeg coords") else: r=0.0 - + return r #----------------------------------------------------------------------------- @@ -291,16 +290,14 @@ def shiftRADec(ra1, dec1, deltaRA, deltaDec): (deltaRA, deltaDec) are arcseconds, and output is decimal degrees. Based on an IDL routine of the same name. - @param ra1: float - @type ra1: R.A. in decimal degrees - @param dec1: float - @type dec1: dec. in decimal degrees - @param deltaRA: float - @type deltaRA: shift in R.A. in arcseconds - @param deltaDec: float - @type deltaDec: shift in dec. in arcseconds - @rtype: float [newRA, newDec] - @return: shifted R.A. and dec. + Args: + ra1 (float): R.A. in decimal degrees + dec1 (float): dec. in decimal degrees + deltaRA (float): shift in R.A. in arcseconds + deltaDec (float): shift in dec. in arcseconds + + Returns: + tuple: shifted R.A. and dec. as (newRA, newDec) in decimal degrees """ @@ -330,20 +327,17 @@ def convertCoords(inputSystem, outputSystem, coordX, coordY, epoch): """Converts specified coordinates (given in decimal degrees) between J2000, B1950, and Galactic. - @type inputSystem: string - @param inputSystem: system of the input coordinates (either "J2000", - "B1950" or "GALACTIC") - @type outputSystem: string - @param outputSystem: system of the returned coordinates (either "J2000", - "B1950" or "GALACTIC") - @type coordX: float - @param coordX: longitude coordinate in decimal degrees, e.g. R. A. - @type coordY: float - @param coordY: latitude coordinate in decimal degrees, e.g. dec. - @type epoch: float - @param epoch: epoch of the input coordinates - @rtype: list - @return: coordinates in decimal degrees in requested output system + Args: + inputSystem (str): system of the input coordinates; one of "J2000", + "B1950", or "GALACTIC" + outputSystem (str): system of the returned coordinates; one of "J2000", + "B1950", or "GALACTIC" + coordX (float): longitude coordinate in decimal degrees, e.g. R.A. + coordY (float): latitude coordinate in decimal degrees, e.g. dec. + epoch (float): epoch of the input coordinates + + Returns: + list: coordinates in decimal degrees in the requested output system """ @@ -367,18 +361,17 @@ def calcRADecSearchBox(RADeg, decDeg, radiusSkyDeg): """Calculates minimum and maximum RA, dec coords needed to define a box enclosing a circle of radius radiusSkyDeg around the given RADeg, decDeg coordinates. Useful for freeform queries of e.g. SDSS, UKIDSS etc.. Uses - L{calcAngSepDeg}, so has the same limitations. - - @type RADeg: float - @param RADeg: RA coordinate of centre of search region - @type decDeg: float - @param decDeg: dec coordinate of centre of search region - @type radiusSkyDeg: float - @param radiusSkyDeg: radius in degrees on the sky used to define search - region - @rtype: list - @return: [RAMin, RAMax, decMin, decMax] - coordinates in decimal degrees - defining search box + calcAngSepDeg, so has the same limitations. + + Args: + RADeg (float): RA coordinate of centre of search region + decDeg (float): dec coordinate of centre of search region + radiusSkyDeg (float): radius in degrees on the sky used to define the + search region + + Returns: + list: [RAMin, RAMax, decMin, decMax] - coordinates in decimal degrees + defining search box """ @@ -439,4 +432,3 @@ def calcRADecSearchBox(RADeg, decDeg, radiusSkyDeg): raise Exception("calcRADecSearchBox failed sanity check") return [RAMin, RAMax, decMin, decMax] - diff --git a/astLib/astImages.py b/astLib/astImages.py index c1eab27..d0b60de 100644 --- a/astLib/astImages.py +++ b/astLib/astImages.py @@ -1,12 +1,12 @@ """Module for simple .fits image tasks (rotation, clipping out sections, making .pngs etc.). -(c) 2007-2018 Matt Hilton +(c) 2007-2018 Matt Hilton Some routines in this module will fail if, e.g., asked to clip a section from a .fits image at a position not found within the image (as determined using the WCS). Where this occurs, the function will return None. An error message will be printed to the console when this happens if astImages.REPORT_ERRORS=True (the default). Testing if an astImages function returns None can be -used to handle errors in scripts. +used to handle errors in scripts. """ @@ -16,7 +16,7 @@ import sys import math from astLib import astWCS -from astropy.io import fits as pyfits +from astropy.io import fits as pyfits try: from scipy import ndimage from scipy import interpolate @@ -32,42 +32,38 @@ #--------------------------------------------------------------------------------------------------- def clipImageSectionWCS(imageData, imageWCS, RADeg, decDeg, clipSizeDeg, returnWCS = True): - """Clips a square or rectangular section from an image array at the given celestial coordinates. - An updated WCS for the clipped section is optionally returned, as well as the x, y pixel + """Clips a square or rectangular section from an image array at the given celestial coordinates. + An updated WCS for the clipped section is optionally returned, as well as the x, y pixel coordinates in the original image corresponding to the clipped section. - + Note that the clip size is specified in degrees on the sky. For projections that have varying - real pixel scale across the map (e.g. CEA), use L{clipUsingRADecCoords} instead. + real pixel scale across the map (e.g. CEA), use clipUsingRADecCoords instead. - Similarly, this routine will not work for a WCS that has polynomial distortion coefficients - in the header (e.g., CTYPE1 = 'RA---TAN-SIP' etc.) - again L{clipUsingRADecCoords} can be used + Similarly, this routine will not work for a WCS that has polynomial distortion coefficients + in the header (e.g., CTYPE1 = 'RA---TAN-SIP' etc.) - again clipUsingRADecCoords can be used in such cases. - - @type imageData: np array - @param imageData: image data array - @type imageWCS: astWCS.WCS - @param imageWCS: astWCS.WCS object - @type RADeg: float - @param RADeg: coordinate in decimal degrees - @type decDeg: float - @param decDeg: coordinate in decimal degrees - @type clipSizeDeg: float or list in format [widthDeg, heightDeg] - @param clipSizeDeg: if float, size of square clipped section in decimal degrees; if list, - size of clipped section in degrees in x, y axes of image respectively - @type returnWCS: bool - @param returnWCS: if True, return an updated WCS for the clipped section - @rtype: dictionary - @return: clipped image section (np array), updated astWCS WCS object for - clipped image section, and coordinates of clipped section in imageData in format - {'data', 'wcs', 'clippedSection'}. - - """ - + + Args: + imageData (numpy.ndarray): image data array + imageWCS (astWCS.WCS): WCS object for the image + RADeg (float): R.A. coordinate in decimal degrees + decDeg (float): dec. coordinate in decimal degrees + clipSizeDeg (float or list): if float, size of square clipped section in decimal degrees; + if list, size of clipped section in degrees as [widthDeg, heightDeg] + returnWCS (bool): if True, return an updated WCS for the clipped section + + Returns: + dict: clipped image section (numpy.ndarray), updated WCS object for the clipped image + section, and coordinates of clipped section in imageData in the format + ``{'data', 'wcs', 'clippedSection'}`` + + """ + imHeight=imageData.shape[0] imWidth=imageData.shape[1] xImScale=imageWCS.getXPixelSizeDeg() yImScale=imageWCS.getYPixelSizeDeg() - + if type(clipSizeDeg) == float: xHalfClipSizeDeg=clipSizeDeg/2.0 yHalfClipSizeDeg=xHalfClipSizeDeg @@ -76,21 +72,21 @@ def clipImageSectionWCS(imageData, imageWCS, RADeg, decDeg, clipSizeDeg, returnW yHalfClipSizeDeg=clipSizeDeg[1]/2.0 else: raise Exception("did not understand clipSizeDeg: should be float, or [widthDeg, heightDeg]") - + xHalfSizePix=xHalfClipSizeDeg/xImScale - yHalfSizePix=yHalfClipSizeDeg/yImScale - + yHalfSizePix=yHalfClipSizeDeg/yImScale + cPixCoords=imageWCS.wcs2pix(RADeg, decDeg) - + cTopLeft=[cPixCoords[0]+xHalfSizePix, cPixCoords[1]+yHalfSizePix] cBottomRight=[cPixCoords[0]-xHalfSizePix, cPixCoords[1]-yHalfSizePix] - + X=[int(round(cTopLeft[0])),int(round(cBottomRight[0]))] Y=[int(round(cTopLeft[1])),int(round(cBottomRight[1]))] - + X.sort() Y.sort() - + if X[0] < 0: X[0]=0 if X[1] > imWidth: @@ -99,7 +95,7 @@ def clipImageSectionWCS(imageData, imageWCS, RADeg, decDeg, clipSizeDeg, returnW Y[0]=0 if Y[1] > imHeight: Y[1]=imHeight - + clippedData=imageData[Y[0]:Y[1],X[0]:X[1]] # Update WCS @@ -113,41 +109,39 @@ def clipImageSectionWCS(imageData, imageWCS, RADeg, decDeg, clipSizeDeg, returnW clippedWCS.header['CRPIX1']=oldCRPIX1-X[0] clippedWCS.header['CRPIX2']=oldCRPIX2-Y[0] clippedWCS.updateFromHeader() - + except KeyError: - + if REPORT_ERRORS == True: - + print("WARNING: astImages.clipImageSectionWCS() : no CRPIX1, CRPIX2 keywords found - not updating clipped image WCS.") - + clippedData=imageData[Y[0]:Y[1],X[0]:X[1]] clippedWCS=imageWCS.copy() else: clippedWCS=None - + return {'data': clippedData, 'wcs': clippedWCS, 'clippedSection': [X[0], X[1], Y[0], Y[1]]} - + #--------------------------------------------------------------------------------------------------- def clipImageSectionPix(imageData, XCoord, YCoord, clipSizePix): """Clips a square or rectangular section from an image array at the given pixel coordinates. - - @type imageData: np array - @param imageData: image data array - @type XCoord: float - @param XCoord: coordinate in pixels - @type YCoord: float - @param YCoord: coordinate in pixels - @type clipSizePix: float or list in format [widthPix, heightPix] - @param clipSizePix: if float, size of square clipped section in pixels; if list, - size of clipped section in pixels in x, y axes of output image respectively - @rtype: np array - @return: clipped image section - - """ - + + Args: + imageData (numpy.ndarray): image data array + XCoord (float): x coordinate in pixels + YCoord (float): y coordinate in pixels + clipSizePix (float or list): if float, size of square clipped section in pixels; if list, + size of clipped section in pixels as [widthPix, heightPix] + + Returns: + numpy.ndarray: clipped image section + + """ + imHeight=imageData.shape[0] imWidth=imageData.shape[1] - + if type(clipSizePix) == float or type(clipSizePix) == int: xHalfClipSizePix=int(round(clipSizePix/2.0)) yHalfClipSizePix=xHalfClipSizePix @@ -156,16 +150,16 @@ def clipImageSectionPix(imageData, XCoord, YCoord, clipSizePix): yHalfClipSizePix=int(round(clipSizePix[1]/2.0)) else: raise Exception("did not understand clipSizePix: should be float, or [widthPix, heightPix]") - + cTopLeft=[XCoord+xHalfClipSizePix, YCoord+yHalfClipSizePix] cBottomRight=[XCoord-xHalfClipSizePix, YCoord-yHalfClipSizePix] - + X=[int(round(cTopLeft[0])),int(round(cBottomRight[0]))] Y=[int(round(cTopLeft[1])),int(round(cBottomRight[1]))] - + X.sort() Y.sort() - + if X[0] < 0: X[0]=0 if X[1] > imWidth: @@ -173,47 +167,45 @@ def clipImageSectionPix(imageData, XCoord, YCoord, clipSizePix): if Y[0] < 0: Y[0]=0 if Y[1] > imHeight: - Y[1]=imHeight - + Y[1]=imHeight + return imageData[Y[0]:Y[1],X[0]:X[1]] - + #--------------------------------------------------------------------------------------------------- def clipRotatedImageSectionWCS(imageData, imageWCS, RADeg, decDeg, clipSizeDeg, returnWCS = True): - """Clips a square or rectangular section from an image array at the given celestial coordinates. + """Clips a square or rectangular section from an image array at the given celestial coordinates. The resulting clip is rotated and/or flipped such that North is at the top, and East appears at the left. An updated WCS for the clipped section is also returned. Note that the alignment of the rotated WCS is currently not perfect - however, it is probably good enough in most - cases for use with L{ImagePlot} for plotting purposes. - + cases for use with ImagePlot for plotting purposes. + Note that the clip size is specified in degrees on the sky. For projections that have varying - real pixel scale across the map (e.g. CEA), use L{clipUsingRADecCoords} instead. - - Similarly, this routine will not work for a WCS that has polynomial distortion coefficients - in the header (e.g., CTYPE1 = 'RA---TAN-SIP' etc.) - again L{clipUsingRADecCoords} can be used + real pixel scale across the map (e.g. CEA), use clipUsingRADecCoords instead. + + Similarly, this routine will not work for a WCS that has polynomial distortion coefficients + in the header (e.g., CTYPE1 = 'RA---TAN-SIP' etc.) - again clipUsingRADecCoords can be used in such cases. - - @type imageData: np array - @param imageData: image data array - @type imageWCS: astWCS.WCS - @param imageWCS: astWCS.WCS object - @type RADeg: float - @param RADeg: coordinate in decimal degrees - @type decDeg: float - @param decDeg: coordinate in decimal degrees - @type clipSizeDeg: float - @param clipSizeDeg: if float, size of square clipped section in decimal degrees; if list, - size of clipped section in degrees in RA, dec. axes of output rotated image respectively - @type returnWCS: bool - @param returnWCS: if True, return an updated WCS for the clipped section - @rtype: dictionary - @return: clipped image section (np array), updated astWCS WCS object for - clipped image section, in format {'data', 'wcs'}. - - @note: Returns 'None' if the requested position is not found within the image. If the image - WCS does not have keywords of the form CD1_1 etc., the output WCS will not be rotated. - + + Args: + imageData (numpy.ndarray): image data array + imageWCS (astWCS.WCS): WCS object for the image + RADeg (float): R.A. coordinate in decimal degrees + decDeg (float): dec. coordinate in decimal degrees + clipSizeDeg (float or list): if float, size of square clipped section in decimal degrees; + if list, size of clipped section in degrees as [widthDeg, heightDeg] in the RA, dec. + axes of the output rotated image respectively + returnWCS (bool): if True, return an updated WCS for the clipped section + + Returns: + dict or None: clipped image section (numpy.ndarray) and updated WCS object in the format + ``{'data', 'wcs'}``, or None if the requested position is not found within the image. + + Note: + If the image WCS does not have keywords of the form CD1_1 etc., the output WCS will not + be rotated. + """ - + halfImageSize=imageWCS.getHalfSizeDeg() imageCentre=imageWCS.getCentreWCSCoords() imScale=imageWCS.getPixelSizeDeg() @@ -226,32 +218,32 @@ def clipRotatedImageSectionWCS(imageData, imageWCS, RADeg, decDeg, clipSizeDeg, yHalfClipSizeDeg=clipSizeDeg[1]/2.0 else: raise Exception("did not understand clipSizeDeg: should be float, or [widthDeg, heightDeg]") - + diagonalHalfSizeDeg=math.sqrt((xHalfClipSizeDeg*xHalfClipSizeDeg) \ +(yHalfClipSizeDeg*yHalfClipSizeDeg)) - + diagonalHalfSizePix=diagonalHalfSizeDeg/imScale - + if RADeg>imageCentre[0]-halfImageSize[0] and RADegimageCentre[1]-halfImageSize[1] and decDeg imWidth: @@ -384,8 +370,8 @@ def clipUsingRADecCoords(imageData, imageWCS, RAMin, RAMax, decMin, decMax, retu if Y[0] < 0: Y[0]=0 if Y[1] > imHeight: - Y[1]=imHeight - + Y[1]=imHeight + clippedData=imageData[Y[0]:Y[1],X[0]:X[1]] # Update WCS @@ -399,45 +385,44 @@ def clipUsingRADecCoords(imageData, imageWCS, RAMin, RAMax, decMin, decMax, retu clippedWCS.header['CRPIX1']=oldCRPIX1-X[0] clippedWCS.header['CRPIX2']=oldCRPIX2-Y[0] clippedWCS.updateFromHeader() - + except KeyError: - + if REPORT_ERRORS == True: - + print("WARNING: astImages.clipUsingRADecCoords() : no CRPIX1, CRPIX2 keywords found - not updating clipped image WCS.") - + clippedData=imageData[Y[0]:Y[1],X[0]:X[1]] clippedWCS=imageWCS.copy() else: clippedWCS=None - + return {'data': clippedData, 'wcs': clippedWCS, 'clippedSection': [X[0], X[1], Y[0], Y[1]]} - + #--------------------------------------------------------------------------------------------------- def scaleImage(imageData, imageWCS, scaleFactor): """Scales image array and WCS by the given scale factor. - - @type imageData: np array - @param imageData: image data array - @type imageWCS: astWCS.WCS - @param imageWCS: astWCS.WCS object - @type scaleFactor: float or list or tuple - @param scaleFactor: factor to resize image by - if tuple or list, in format - [x scale factor, y scale factor] - @rtype: dictionary - @return: image data (np array), updated astWCS WCS object for image, in format {'data', 'wcs'}. - + + Args: + imageData (numpy.ndarray): image data array + imageWCS (astWCS.WCS): WCS object for the image + scaleFactor (float or list or tuple): factor to resize image by; if tuple or list, in + format [x scale factor, y scale factor] + + Returns: + dict: image data (numpy.ndarray) and updated WCS object in format ``{'data', 'wcs'}`` + """ if type(scaleFactor) == int or type(scaleFactor) == float: - scaleFactor=[float(scaleFactor), float(scaleFactor)] + scaleFactor=[float(scaleFactor), float(scaleFactor)] scaledData=ndimage.zoom(imageData, scaleFactor) - + # Changed below because ndimage.zoom now uses round instead of int (since scipy 0.13.0) # NOTE: np axes order flips order compared to scaleFactor trueScaleFactor=np.array(scaledData.shape, dtype = float) / np.array(imageData.shape, dtype = float) offset=0. - + # Rescale WCS try: oldCRPIX1=imageWCS.header['CRPIX1'] @@ -445,7 +430,7 @@ def scaleImage(imageData, imageWCS, scaleFactor): CD11=imageWCS.header['CD1_1'] CD21=imageWCS.header['CD2_1'] CD12=imageWCS.header['CD1_2'] - CD22=imageWCS.header['CD2_2'] + CD22=imageWCS.header['CD2_2'] except KeyError: # Try the older FITS header format try: @@ -476,47 +461,49 @@ def scaleImage(imageData, imageWCS, scaleFactor): scaledWCS.header['CD1_2']=scaledCDMatrix[0][1] scaledWCS.header['CD2_2']=scaledCDMatrix[1][1] scaledWCS.updateFromHeader() - + return {'data': scaledData, 'wcs': scaledWCS} - + #--------------------------------------------------------------------------------------------------- def intensityCutImage(imageData, cutLevels): """Creates a matplotlib.pylab plot of an image array with the specified cuts in intensity - applied. This routine is used by L{saveBitmap} and L{saveContourOverlayBitmap}, which both + applied. This routine is used by saveBitmap and saveContourOverlayBitmap, which both produce output as .png, .jpg, etc. images. - - @type imageData: np array - @param imageData: image data array - @type cutLevels: list - @param cutLevels: sets the image scaling - available options: - - pixel values: cutLevels=[low value, high value]. - - histogram equalisation: cutLevels=["histEq", number of bins ( e.g. 1024)] - - relative: cutLevels=["relative", cut per cent level (e.g. 99.5)] - - smart: cutLevels=["smart", cut per cent level (e.g. 99.5)] - ["smart", 99.5] seems to provide good scaling over a range of different images. - @rtype: dictionary - @return: image section (np.array), matplotlib image normalisation (matplotlib.colors.Normalize), in the format {'image', 'norm'}. - - @note: If cutLevels[0] == "histEq", then only {'image'} is returned. - + + Args: + imageData (numpy.ndarray): image data array + cutLevels (list): sets the image scaling - available options: + + - pixel values: ``[low value, high value]`` + - histogram equalisation: ``["histEq", number of bins (e.g. 1024)]`` + - relative: ``["relative", cut per cent level (e.g. 99.5)]`` + - smart: ``["smart", cut per cent level (e.g. 99.5)]`` + + ``["smart", 99.5]`` seems to provide good scaling over a range of different images. + + Returns: + dict: image section (numpy.ndarray) and matplotlib image normalisation + (matplotlib.colors.Normalize) in the format ``{'image', 'norm'}``. If + ``cutLevels[0] == "histEq"``, only ``{'image'}`` is returned. + """ - + oImWidth=imageData.shape[1] oImHeight=imageData.shape[0] - + # Optional histogram equalisation if cutLevels[0]=="histEq": - + imageData=histEq(imageData, cutLevels[1]) anorm=pylab.Normalize(imageData.min(), imageData.max()) - + elif cutLevels[0]=="relative": - + # this turns image data into 1D array then sorts - sorted=np.sort(np.ravel(imageData)) + sorted=np.sort(np.ravel(imageData)) maxValue=sorted.max() minValue=sorted.min() - + # want to discard the top and bottom specified topCutIndex=len(sorted-1) \ -int(math.floor(float((100.0-cutLevels[1])/100.0)*len(sorted-1))) @@ -524,17 +511,17 @@ def intensityCutImage(imageData, cutLevels): topCut=sorted[topCutIndex] bottomCut=sorted[bottomCutIndex] anorm=pylab.Normalize(bottomCut, topCut) - + elif cutLevels[0]=="smart": - + # this turns image data into 1Darray then sorts - sorted=np.sort(np.ravel(imageData)) + sorted=np.sort(np.ravel(imageData)) maxValue=sorted.max() minValue=sorted.min() numBins=10000 # 0.01 per cent accuracy binWidth=(maxValue-minValue)/float(numBins) histogram=ndimage.histogram(sorted, minValue, maxValue, numBins) - + # Find the bin with the most pixels in it, set that as our minimum # Then search through the bins until we get to a bin with more/or the same number of # pixels in it than the previous one. @@ -544,27 +531,27 @@ def intensityCutImage(imageData, cutLevels): backgroundValue=histogram.max() foundBackgroundBin=False foundTopBin=False - lastBin=-10000 + lastBin=-10000 for i in range(len(histogram)): - + if histogram[i]>=lastBin and foundBackgroundBin==True: - - # Added a fudge here to stop us picking for top bin a bin within + + # Added a fudge here to stop us picking for top bin a bin within # 10 percent of the background pixel value if (minValue+(binWidth*i))>bottomBinValue*1.1: topBinValue=minValue+(binWidth*i) foundTopBin=True break - + if histogram[i]==backgroundValue and foundBackgroundBin==False: bottomBinValue=minValue+(binWidth*i) foundBackgroundBin=True lastBin=histogram[i] - + if foundTopBin==False: topBinValue=maxValue - + #Now we apply relative scaling to this smartClipped=np.clip(sorted, bottomBinValue, topBinValue) topCutIndex=len(smartClipped-1) \ @@ -574,10 +561,10 @@ def intensityCutImage(imageData, cutLevels): bottomCut=smartClipped[bottomCutIndex] anorm=pylab.Normalize(bottomCut, topCut) else: - + # Normalise using given cut levels anorm=pylab.Normalize(cutLevels[0], cutLevels[1]) - + if cutLevels[0]=="histEq": return {'image': imageData.copy()} else: @@ -587,18 +574,17 @@ def intensityCutImage(imageData, cutLevels): def resampleToTanProjection(imageData, imageWCS, outputPixDimensions=[600, 600]): """Resamples an image and WCS to a tangent plane projection. Purely for plotting purposes (e.g., ensuring RA, dec. coordinate axes perpendicular). - - @type imageData: np array - @param imageData: image data array - @type imageWCS: astWCS.WCS - @param imageWCS: astWCS.WCS object - @type outputPixDimensions: list - @param outputPixDimensions: [width, height] of output image in pixels - @rtype: dictionary - @return: image data (np array), updated astWCS WCS object for image, in format {'data', 'wcs'}. - + + Args: + imageData (numpy.ndarray): image data array + imageWCS (astWCS.WCS): WCS object for the image + outputPixDimensions (list): [width, height] of output image in pixels + + Returns: + dict: image data (numpy.ndarray) and updated WCS object in format ``{'data', 'wcs'}`` + """ - + RADeg, decDeg=imageWCS.getCentreWCSCoords() xPixelScale=imageWCS.getXPixelSizeDeg() yPixelScale=imageWCS.getYPixelSizeDeg() @@ -626,46 +612,43 @@ def resampleToTanProjection(imageData, imageWCS, outputPixDimensions=[600, 600]) newWCS=astWCS.WCS(newHead, mode='pyfits') newImage=np.zeros([ySizePix, xSizePix]) - tanImage=resampleToWCS(newImage, newWCS, imageData, imageWCS, highAccuracy=True, + tanImage=resampleToWCS(newImage, newWCS, imageData, imageWCS, highAccuracy=True, onlyOverlapping=False) - - return tanImage - + + return tanImage + #--------------------------------------------------------------------------------------------------- def resampleToWCS(im1Data, im1WCS, im2Data, im2WCS, highAccuracy = False, onlyOverlapping = True): - """Resamples data corresponding to second image (with data im2Data, WCS im2WCS) onto the WCS - of the first image (im1Data, im1WCS). The output, resampled image is of the pixel same - dimensions of the first image. This routine is for assisting in plotting - performing - photometry on the output is not recommended. - + """Resamples data corresponding to second image (with data im2Data, WCS im2WCS) onto the WCS + of the first image (im1Data, im1WCS). The output, resampled image is of the pixel same + dimensions of the first image. This routine is for assisting in plotting - performing + photometry on the output is not recommended. + Set highAccuracy == True to sample every corresponding pixel in each image; otherwise only every nth pixel (where n is the ratio of the image scales) will be sampled, with values in between being set using a linear interpolation (much faster). - + Set onlyOverlapping == True to speed up resampling by only resampling the overlapping area defined by both image WCSs. - - @type im1Data: np array - @param im1Data: image data array for first image - @type im1WCS: astWCS.WCS - @param im1WCS: astWCS.WCS object corresponding to im1Data - @type im2Data: np array - @param im2Data: image data array for second image (to be resampled to match first image) - @type im2WCS: astWCS.WCS - @param im2WCS: astWCS.WCS object corresponding to im2Data - @type highAccuracy: bool - @param highAccuracy: if True, sample every corresponding pixel in each image; otherwise, sample - every nth pixel, where n = the ratio of the image scales. - @type onlyOverlapping: bool - @param onlyOverlapping: if True, only consider the overlapping area defined by both image WCSs - (speeds things up) - @rtype: dictionary - @return: np image data array and associated WCS in format {'data', 'wcs'} - + + Args: + im1Data (numpy.ndarray): image data array for first image + im1WCS (astWCS.WCS): WCS object corresponding to im1Data + im2Data (numpy.ndarray): image data array for second image (to be resampled to match + first image) + im2WCS (astWCS.WCS): WCS object corresponding to im2Data + highAccuracy (bool): if True, sample every corresponding pixel in each image; otherwise, + sample every nth pixel, where n = the ratio of the image scales + onlyOverlapping (bool): if True, only consider the overlapping area defined by both image + WCSs (speeds things up) + + Returns: + dict: numpy image data array and associated WCS in format ``{'data', 'wcs'}`` + """ - + resampledData=np.zeros(im1Data.shape) - + # Find overlap - speed things up # But have a border so as not to require the overlap to be perfect # There's also no point in oversampling image 1 if it's much higher res than image 2 @@ -685,7 +668,7 @@ def resampleToWCS(im1Data, im1WCS, im2Data, im2WCS, highAccuracy = False, onlyOv else: xPixStep=1 yPixStep=1 - + if onlyOverlapping == True: overlap=astWCS.findWCSOverlap(im1WCS, im2WCS) xOverlap=[overlap['wcs1Pix'][0], overlap['wcs1Pix'][1]] @@ -716,7 +699,7 @@ def resampleToWCS(im1Data, im1WCS, im2Data, im2WCS, highAccuracy = False, onlyOv xMax=im1Data.shape[1] yMin=0 yMax=im1Data.shape[0] - + for x in range(xMin, xMax, xPixStep): for y in range(yMin, yMax, yPixStep): RA, dec=im1WCS.pix2wcs(x, y) @@ -738,44 +721,44 @@ def resampleToWCS(im1Data, im1WCS, im2Data, im2WCS, highAccuracy = False, onlyOv index2data=interpolate.interp1d(np.arange(0, vals.shape[0], 1), vals) interpedVals=index2data(np.arange(0, vals.shape[0]-1, 1.0/yPixStep)) resampledData[yMin:yMin+interpedVals.shape[0], col]=interpedVals - + # Note: should really just copy im1WCS keywords into im2WCS and return that # Only a problem if we're using this for anything other than plotting return {'data': resampledData, 'wcs': im1WCS.copy()} - + #--------------------------------------------------------------------------------------------------- def generateContourOverlay(backgroundImageData, backgroundImageWCS, contourImageData, contourImageWCS, \ contourLevels, contourSmoothFactor = 0, highAccuracy = False): - """Rescales an image array to be used as a contour overlay to have the same dimensions as the - background image, and generates a set of contour levels. The image array from which the contours - are to be generated will be resampled to the same dimensions as the background image data, and - can be optionally smoothed using a Gaussian filter. The sigma of the Gaussian filter + """Rescales an image array to be used as a contour overlay to have the same dimensions as the + background image, and generates a set of contour levels. The image array from which the contours + are to be generated will be resampled to the same dimensions as the background image data, and + can be optionally smoothed using a Gaussian filter. The sigma of the Gaussian filter (contourSmoothFactor) is specified in arcsec. - - @type backgroundImageData: np array - @param backgroundImageData: background image data array - @type backgroundImageWCS: astWCS.WCS - @param backgroundImageWCS: astWCS.WCS object of the background image data array - @type contourImageData: np array - @param contourImageData: image data array from which contours are to be generated - @type contourImageWCS: astWCS.WCS - @param contourImageWCS: astWCS.WCS object corresponding to contourImageData - @type contourLevels: list - @param contourLevels: sets the contour levels - available options: - - values: contourLevels=[list of values specifying each level] - - linear spacing: contourLevels=['linear', min level value, max level value, number - of levels] - can use "min", "max" to automatically set min, max levels from image data - - log spacing: contourLevels=['log', min level value, max level value, number of - levels] - can use "min", "max" to automatically set min, max levels from image data - @type contourSmoothFactor: float - @param contourSmoothFactor: standard deviation (in arcsec) of Gaussian filter for - pre-smoothing of contour image data (set to 0 for no smoothing) - @type highAccuracy: bool - @param highAccuracy: if True, sample every corresponding pixel in each image; otherwise, sample - every nth pixel, where n = the ratio of the image scales. - - """ - + + Args: + backgroundImageData (numpy.ndarray): background image data array + backgroundImageWCS (astWCS.WCS): WCS object of the background image + contourImageData (numpy.ndarray): image data array from which contours are to be generated + contourImageWCS (astWCS.WCS): WCS object corresponding to contourImageData + contourLevels (list): sets the contour levels - available options: + + - values: list of values specifying each level + - linear spacing: ``['linear', min level value, max level value, number of levels]`` + - can use ``"min"``, ``"max"`` to automatically set min, max levels from image data + - log spacing: ``['log', min level value, max level value, number of levels]`` + - can use ``"min"``, ``"max"`` to automatically set min, max levels from image data + + contourSmoothFactor (float): standard deviation (in arcsec) of Gaussian filter for + pre-smoothing of contour image data (set to 0 for no smoothing) + highAccuracy (bool): if True, sample every corresponding pixel in each image; otherwise, + sample every nth pixel, where n = the ratio of the image scales + + Returns: + dict: resampled contour image array and contour levels in the format + ``{'scaledImage', 'contourLevels'}`` + + """ + # For compromise between speed and accuracy, scale a copy of the background # image down to a scale that is one pixel = 1/5 of a pixel in the contour image # But only do this if it has CDij keywords as we know how to scale those @@ -783,13 +766,13 @@ def generateContourOverlay(backgroundImageData, backgroundImageWCS, contourImage xScaleFactor=backgroundImageWCS.getXPixelSizeDeg()/(contourImageWCS.getXPixelSizeDeg()/5.0) yScaleFactor=backgroundImageWCS.getYPixelSizeDeg()/(contourImageWCS.getYPixelSizeDeg()/5.0) scaledBackground=scaleImage(backgroundImageData, backgroundImageWCS, (xScaleFactor, yScaleFactor)) - scaled=resampleToWCS(scaledBackground['data'], scaledBackground['wcs'], + scaled=resampleToWCS(scaledBackground['data'], scaledBackground['wcs'], contourImageData, contourImageWCS, highAccuracy = highAccuracy) scaledContourData=scaled['data'] scaledContourWCS=scaled['wcs'] scaledBackground=True else: - scaled=resampleToWCS(backgroundImageData, backgroundImageWCS, + scaled=resampleToWCS(backgroundImageData, backgroundImageWCS, contourImageData, contourImageWCS, highAccuracy = highAccuracy) scaledContourData=scaled['data'] scaledContourWCS=scaled['wcs'] @@ -798,7 +781,7 @@ def generateContourOverlay(backgroundImageData, backgroundImageWCS, contourImage if contourSmoothFactor != None and contourSmoothFactor > 0: sigmaPix=(contourSmoothFactor/3600.0)/scaledContourWCS.getPixelSizeDeg() scaledContourData=ndimage.gaussian_filter(scaledContourData, sigmaPix) - + # Various ways of setting the contour levels # If just a list is passed in, use those instead if contourLevels[0] == "linear": @@ -809,14 +792,14 @@ def generateContourOverlay(backgroundImageData, backgroundImageWCS, contourImage if contourLevels[2] == "max": xMax=contourImageData.flatten().max() else: - xMax=float(contourLevels[2]) + xMax=float(contourLevels[2]) nLevels=contourLevels[3] xStep=(xMax-xMin)/(nLevels-1) cLevels=[] for j in range(nLevels+1): level=xMin+j*xStep cLevels.append(level) - + elif contourLevels[0] == "log": if contourLevels[1] == "min": xMin=contourImageData.flatten().min() @@ -825,7 +808,7 @@ def generateContourOverlay(backgroundImageData, backgroundImageWCS, contourImage if contourLevels[2] == "max": xMax=contourImageData.flatten().max() else: - xMax=float(contourLevels[2]) + xMax=float(contourLevels[2]) if xMin <= 0.0: raise Exception("minimum contour level set to <= 0 and log scaling chosen.") xLogMin=math.log10(xMin) @@ -836,67 +819,65 @@ def generateContourOverlay(backgroundImageData, backgroundImageWCS, contourImage prevLevel=0 for j in range(nLevels+1): level=math.pow(10, xLogMin+j*xLogStep) - cLevels.append(level) - + cLevels.append(level) + else: cLevels=contourLevels - - # Now blow the contour image data back up to the size of the original image + + # Now blow the contour image data back up to the size of the original image if scaledBackground == True: scaledBack=scaleImage(scaledContourData, scaledContourWCS, (1.0/xScaleFactor, 1.0/yScaleFactor))['data'] else: scaledBack=scaledContourData - + return {'scaledImage': scaledBack, 'contourLevels': cLevels} - + #--------------------------------------------------------------------------------------------------- def saveBitmap(outputFileName, imageData, cutLevels, size, colorMapName): """Makes a bitmap image from an image array; the image format is specified by the filename extension. (e.g. ".jpg" =JPEG, ".png"=PNG). - - @type outputFileName: string - @param outputFileName: filename of output bitmap image - @type imageData: np array - @param imageData: image data array - @type cutLevels: list - @param cutLevels: sets the image scaling - available options: - - pixel values: cutLevels=[low value, high value]. - - histogram equalisation: cutLevels=["histEq", number of bins ( e.g. 1024)] - - relative: cutLevels=["relative", cut per cent level (e.g. 99.5)] - - smart: cutLevels=["smart", cut per cent level (e.g. 99.5)] - ["smart", 99.5] seems to provide good scaling over a range of different images. - @type size: int - @param size: size of output image in pixels - @type colorMapName: string - @param colorMapName: name of a standard matplotlib colormap, e.g. "hot", "cool", "gray" - etc. (do "help(pylab.colormaps)" in the Python interpreter to see available options) - - """ - + + Args: + outputFileName (str): filename of output bitmap image + imageData (numpy.ndarray): image data array + cutLevels (list): sets the image scaling - available options: + + - pixel values: ``[low value, high value]`` + - histogram equalisation: ``["histEq", number of bins (e.g. 1024)]`` + - relative: ``["relative", cut per cent level (e.g. 99.5)]`` + - smart: ``["smart", cut per cent level (e.g. 99.5)]`` + + ``["smart", 99.5]`` seems to provide good scaling over a range of different images. + + size (int): size of output image in pixels + colorMapName (str): name of a standard matplotlib colormap, e.g. "hot", "cool", "gray" + + """ + cut=intensityCutImage(imageData, cutLevels) - + # Make plot aspectR=float(cut['image'].shape[0])/float(cut['image'].shape[1]) pylab.figure(figsize=(10,10*aspectR)) pylab.axes([0,0,1,1]) - + try: colorMap=pylab.cm.get_cmap(colorMapName) except AssertionError: raise Exception(colorMapName+" is not a defined matplotlib colormap.") - + if cutLevels[0]=="histEq": pylab.imshow(cut['image'], interpolation="bilinear", origin='lower', cmap=colorMap) - + else: pylab.imshow(cut['image'], interpolation="bilinear", norm=cut['norm'], origin='lower', cmap=colorMap) pylab.axis("off") - - pylab.savefig("out_astImages.png") + + pylab.savefig("out_astImages.png") pylab.close("all") - + try: from PIL import Image except: @@ -904,7 +885,7 @@ def saveBitmap(outputFileName, imageData, cutLevels, size, colorMapName): im=Image.open("out_astImages.png") im.thumbnail((int(size),int(size))) im.save(outputFileName) - + os.remove("out_astImages.png") #--------------------------------------------------------------------------------------------------- @@ -914,64 +895,56 @@ def saveContourOverlayBitmap(outputFileName, backgroundImageData, backgroundImag """Makes a bitmap image from an image array, with a set of contours generated from a second image array overlaid. The image format is specified by the file extension (e.g. ".jpg"=JPEG, ".png"=PNG). The image array from which the contours are to be generated - can optionally be pre-smoothed using a Gaussian filter. - - @type outputFileName: string - @param outputFileName: filename of output bitmap image - @type backgroundImageData: np array - @param backgroundImageData: background image data array - @type backgroundImageWCS: astWCS.WCS - @param backgroundImageWCS: astWCS.WCS object of the background image data array - @type cutLevels: list - @param cutLevels: sets the image scaling - available options: - - pixel values: cutLevels=[low value, high value]. - - histogram equalisation: cutLevels=["histEq", number of bins ( e.g. 1024)] - - relative: cutLevels=["relative", cut per cent level (e.g. 99.5)] - - smart: cutLevels=["smart", cut per cent level (e.g. 99.5)] - ["smart", 99.5] seems to provide good scaling over a range of different images. - @type size: int - @param size: size of output image in pixels - @type colorMapName: string - @param colorMapName: name of a standard matplotlib colormap, e.g. "hot", "cool", "gray" - etc. (do "help(pylab.colormaps)" in the Python interpreter to see available options) - @type contourImageData: np array - @param contourImageData: image data array from which contours are to be generated - @type contourImageWCS: astWCS.WCS - @param contourImageWCS: astWCS.WCS object corresponding to contourImageData - @type contourSmoothFactor: float - @param contourSmoothFactor: standard deviation (in pixels) of Gaussian filter for - pre-smoothing of contour image data (set to 0 for no smoothing) - @type contourLevels: list - @param contourLevels: sets the contour levels - available options: - - values: contourLevels=[list of values specifying each level] - - linear spacing: contourLevels=['linear', min level value, max level value, number - of levels] - can use "min", "max" to automatically set min, max levels from image data - - log spacing: contourLevels=['log', min level value, max level value, number of - levels] - can use "min", "max" to automatically set min, max levels from image data - @type contourColor: string - @param contourColor: color of the overlaid contours, specified by the name of a standard - matplotlib color, e.g., "black", "white", "cyan" - etc. (do "help(pylab.colors)" in the Python interpreter to see available options) - @type contourWidth: int - @param contourWidth: width of the overlaid contours - - """ - + can optionally be pre-smoothed using a Gaussian filter. + + Args: + outputFileName (str): filename of output bitmap image + backgroundImageData (numpy.ndarray): background image data array + backgroundImageWCS (astWCS.WCS): WCS object of the background image + cutLevels (list): sets the image scaling - available options: + + - pixel values: ``[low value, high value]`` + - histogram equalisation: ``["histEq", number of bins (e.g. 1024)]`` + - relative: ``["relative", cut per cent level (e.g. 99.5)]`` + - smart: ``["smart", cut per cent level (e.g. 99.5)]`` + + ``["smart", 99.5]`` seems to provide good scaling over a range of different images. + + size (int): size of output image in pixels + colorMapName (str): name of a standard matplotlib colormap, e.g. "hot", "cool", "gray" + contourImageData (numpy.ndarray): image data array from which contours are to be generated + contourImageWCS (astWCS.WCS): WCS object corresponding to contourImageData + contourSmoothFactor (float): standard deviation (in pixels) of Gaussian filter for + pre-smoothing of contour image data (set to 0 for no smoothing) + contourLevels (list): sets the contour levels - available options: + + - values: list of values specifying each level + - linear spacing: ``['linear', min level value, max level value, number of levels]`` + - can use ``"min"``, ``"max"`` to automatically set min, max levels from image data + - log spacing: ``['log', min level value, max level value, number of levels]`` + - can use ``"min"``, ``"max"`` to automatically set min, max levels from image data + + contourColor (str): color of the overlaid contours, specified by the name of a standard + matplotlib color, e.g., "black", "white", "cyan" + contourWidth (int): width of the overlaid contours + + """ + cut=intensityCutImage(backgroundImageData, cutLevels) - + # Make plot of just the background image aspectR=float(cut['image'].shape[0])/float(cut['image'].shape[1]) pylab.figure(figsize=(10,10*aspectR)) pylab.axes([0,0,1,1]) - + try: colorMap=pylab.cm.get_cmap(colorMapName) except AssertionError: raise Exception(colorMapName+" is not a defined matplotlib colormap.") - + if cutLevels[0]=="histEq": pylab.imshow(cut['image'], interpolation="bilinear", origin='lower', cmap=colorMap) - + else: pylab.imshow(cut['image'], interpolation="bilinear", norm=cut['norm'], origin='lower', cmap=colorMap) @@ -981,78 +954,74 @@ def saveContourOverlayBitmap(outputFileName, backgroundImageData, backgroundImag # Add the contours contourData=generateContourOverlay(backgroundImageData, backgroundImageWCS, contourImageData, \ contourImageWCS, contourLevels, contourSmoothFactor) - + pylab.contour(contourData['scaledImage'], contourData['contourLevels'], colors=contourColor, - linewidths=contourWidth) - - pylab.savefig("out_astImages.png") + linewidths=contourWidth) + + pylab.savefig("out_astImages.png") pylab.close("all") - + try: from PIL import Image except: raise Exception("astImages.saveContourOverlayBitmap requires the Python Imaging Library to be installed") - + im=Image.open("out_astImages.png") im.thumbnail((int(size),int(size))) im.save(outputFileName) - + os.remove("out_astImages.png") - + #--------------------------------------------------------------------------------------------------- def saveFITS(outputFileName, imageData, imageWCS = None): """Writes an image array to a new .fits file. - - @type outputFileName: string - @param outputFileName: filename of output FITS image - @type imageData: np array - @param imageData: image data array - @type imageWCS: astWCS.WCS object - @param imageWCS: image WCS object - - @note: If imageWCS=None, the FITS image will be written with a rudimentary header containing - no meta data. - + + Args: + outputFileName (str): filename of output FITS image + imageData (numpy.ndarray): image data array + imageWCS (astWCS.WCS or None): image WCS object; if None, the FITS image will be + written with a rudimentary header containing no meta data. + """ - + if os.path.exists(outputFileName): os.remove(outputFileName) - + # There a fudge here for handling both pyfits and astropy.io.fits headers # Removed from version 0.10.0+ (supporting astropy only) if imageWCS != None: hdu=pyfits.PrimaryHDU(None, imageWCS.header) else: hdu=pyfits.PrimaryHDU(None, None) - + newImg=pyfits.HDUList() hdu.data=imageData newImg.append(hdu) newImg.writeto(outputFileName) newImg.close() - + #--------------------------------------------------------------------------------------------------- def histEq(inputArray, numBins): - """Performs histogram equalisation of the input np array. - - @type inputArray: np array - @param inputArray: image data array - @type numBins: int - @param numBins: number of bins in which to perform the operation (e.g. 1024) - @rtype: np array - @return: image data array - + """Performs histogram equalisation of the input numpy array. + + Args: + inputArray (numpy.ndarray): image data array + numBins (int): number of bins in which to perform the operation (e.g. 1024) + + Returns: + numpy.ndarray: equalised image data array + """ - + imageData=inputArray - + # histogram equalisation: we want an equal number of pixels in each intensity range - sortedDataIntensities=np.sort(np.ravel(imageData)) + sortedDataIntensities=np.sort(np.ravel(imageData)) median=np.median(sortedDataIntensities) - + # Make cumulative histogram of data values, simple min-max used to set bin sizes and range dataCumHist=np.zeros(numBins) - minIntensity=sortedDataIntensities.min() + minIntensity=sortedDataIntensities.min() maxIntensity=sortedDataIntensities.max() histRange=maxIntensity-minIntensity binWidth=histRange/float(numBins-1) @@ -1063,54 +1032,53 @@ def histEq(inputArray, numBins): onesRange=list(range(binNumber, numBins)) np.put(addArray, onesRange, onesArray) dataCumHist=dataCumHist+addArray - + # Make ideal cumulative histogram idealValue=dataCumHist.max()/float(numBins) idealCumHist=np.arange(idealValue, dataCumHist.max()+idealValue, idealValue) - + # Map the data to the ideal for y in range(imageData.shape[0]): for x in range(imageData.shape[1]): # Get index corresponding to dataIntensity intensityBin=int(math.ceil((imageData[y][x]-minIntensity)/binWidth)) - + # Guard against rounding errors (happens rarely I think) if intensityBin<0: intensityBin=0 if intensityBin>len(dataCumHist)-1: intensityBin=len(dataCumHist)-1 - + # Get the cumulative frequency corresponding intensity level in the data dataCumFreq=dataCumHist[intensityBin] - + # Get the index of the corresponding ideal cumulative frequency idealBin=np.searchsorted(idealCumHist, dataCumFreq) idealIntensity=(idealBin*binWidth)+minIntensity - imageData[y][x]=idealIntensity - + imageData[y][x]=idealIntensity + return imageData #--------------------------------------------------------------------------------------------------- def normalise(inputArray, clipMinMax): """Clips the inputArray in intensity and normalises the array such that minimum and maximum - values are 0, 1. Clip in intensity is specified by clipMinMax, a list in the format - [clipMin, clipMax] - + values are 0, 1. Clip in intensity is specified by clipMinMax, a list in the format + [clipMin, clipMax]. + Used for normalising image arrays so that they can be turned into RGB arrays that matplotlib - can plot (see L{astPlots.ImagePlot}). - - @type inputArray: np array - @param inputArray: image data array - @type clipMinMax: list - @param clipMinMax: [minimum value of clipped array, maximum value of clipped array] - @rtype: np array - @return: normalised array with minimum value 0, maximum value 1 + can plot (see astPlots.ImagePlot). + + Args: + inputArray (numpy.ndarray): image data array + clipMinMax (list): [minimum value of clipped array, maximum value of clipped array] + + Returns: + numpy.ndarray: normalised array with minimum value 0, maximum value 1 """ clipped=inputArray.clip(clipMinMax[0], clipMinMax[1]) slope=1.0/(clipMinMax[1]-clipMinMax[0]) intercept=-clipMinMax[0]*slope clipped=clipped*slope+intercept - + return clipped - diff --git a/astLib/astPlots.py b/astLib/astPlots.py index 5505880..b789921 100755 --- a/astLib/astPlots.py +++ b/astLib/astPlots.py @@ -2,8 +2,8 @@ (c) 2007-2024 Matt Hilton -This module provides the matplotlib powered ImagePlot class, which is designed to be flexible. -ImagePlots can have RA, Dec. coordinate axes, contour overlays, and have objects marked in them, +This module provides the matplotlib powered ImagePlot class, which is designed to be flexible. +ImagePlots can have RA, Dec. coordinate axes, contour overlays, and have objects marked in them, using WCS coordinates. RGB plots are supported too. """ @@ -27,17 +27,17 @@ def u(x): else: def u(x): return x - -DEC_TICK_STEPS=[{'deg': 1.0/60.0/60.0, 'unit': "s"}, + +DEC_TICK_STEPS=[{'deg': 1.0/60.0/60.0, 'unit': "s"}, {'deg': 2.0/60.0/60.0, 'unit': "s"}, - {'deg': 5.0/60.0/60.0, 'unit': "s"}, + {'deg': 5.0/60.0/60.0, 'unit': "s"}, {'deg': 10.0/60.0/60.0, 'unit': "s"}, {'deg': 30.0/60.0/60.0, 'unit': "s"}, {'deg': 1.0/60.0, 'unit': "m"}, {'deg': 2.0/60.0, 'unit': "m"}, {'deg': 5.0/60.0, 'unit': "m"}, {'deg': 15.0/60.0, 'unit': "m"}, - {'deg': 30.0/60.0, 'unit': "m"}, + {'deg': 30.0/60.0, 'unit': "m"}, {'deg': 1.0, 'unit': "d"}, {'deg': 2.0, 'unit': "d"}, {'deg': 4.0, 'unit': "d"}, @@ -50,9 +50,9 @@ def u(x): RA_TICK_STEPS=[ {'deg': (0.5/60.0/60.0/24.0)*360.0, 'unit': "s"}, {'deg': (1.0/60.0/60.0/24.0)*360.0, 'unit': "s"}, - {'deg': (2.0/60.0/60.0/24.0)*360.0, 'unit': "s"}, - {'deg': (4.0/60.0/60.0/24.0)*360.0, 'unit': "s"}, - {'deg': (5.0/60.0/60.0/24.0)*360.0, 'unit': "s"}, + {'deg': (2.0/60.0/60.0/24.0)*360.0, 'unit': "s"}, + {'deg': (4.0/60.0/60.0/24.0)*360.0, 'unit': "s"}, + {'deg': (5.0/60.0/60.0/24.0)*360.0, 'unit': "s"}, {'deg': (10.0/60.0/60.0/24.0)*360.0, 'unit': "s"}, {'deg': (20.0/60.0/60.0/24.0)*360.0, 'unit': "s"}, {'deg': (30.0/60.0/60.0/24.0)*360.0, 'unit': "s"}, @@ -61,7 +61,7 @@ def u(x): {'deg': (5.0/60.0/24.0)*360.0, 'unit': "m"}, {'deg': (10.0/60.0/24.0)*360.0, 'unit': "m"}, {'deg': (20.0/60.0/24.0)*360.0, 'unit': "m"}, - {'deg': (30.0/60.0/24.0)*360.0, 'unit': "m"}, + {'deg': (30.0/60.0/24.0)*360.0, 'unit': "m"}, {'deg': (1.0/24.0)*360.0, 'unit': "h"}, {'deg': (3.0/24.0)*360.0, 'unit': "h"}, {'deg': (6.0/24.0)*360.0, 'unit': "h"}, @@ -82,21 +82,21 @@ def u(x): class ImagePlot: """This class describes a matplotlib image plot containing an astronomical image with an associated WCS. - - Objects within the image boundaries can be marked by passing their WCS coordinates to - L{ImagePlot.addPlotObjects}. - - Other images can be overlaid using L{ImagePlot.addContourOverlay}. - + + Objects within the image boundaries can be marked by passing their WCS coordinates to + addPlotObjects. + + Other images can be overlaid using addContourOverlay. + For images rotated with North at the top, East at the left (as can be done using - L{astImages.clipRotatedImageSectionWCS} or L{astImages.resampleToTanProjection}, WCS coordinate - axes can be plotted, with tick marks set appropriately for the image size. Otherwise, a compass + astImages.clipRotatedImageSectionWCS or astImages.resampleToTanProjection), WCS coordinate + axes can be plotted, with tick marks set appropriately for the image size. Otherwise, a compass can be plotted showing the directions of North and East in the image. RGB images are also supported. - + The plot can of course be tweaked further after creation using matplotlib/pylab commands. - + """ def __init__(self, imageData, imageWCS, axes = [0.1,0.1,0.8,0.8], \ cutLevels = ["smart", 99.5], colorMapName = "gray", title = None, axesLabels = "sexagesimal", \ @@ -104,61 +104,52 @@ def __init__(self, imageData, imageWCS, axes = [0.1,0.1,0.8,0.8], \ colorBar = False, interpolation = "bilinear"): """Makes an ImagePlot from the given image array and astWCS. For coordinate axes to work, the image and WCS should have been rotated such that East is at the left, North is at the top - (see e.g. L{astImages.clipRotatedImageSectionWCS}, or L{astImages.resampleToTanProjection}). - + (see e.g. astImages.clipRotatedImageSectionWCS, or astImages.resampleToTanProjection). + If imageData is given as a list in the format [r, g, b], a color RGB plot will be made. However, in this case the cutLevels must be specified manually for each component as a list - i.e. cutLevels = [[r min, r max], [g min, g max], [b min, b max]]. In this case of course, the colorMap will be ignored. All r, g, b image arrays must have the same dimensions. - + Set axesLabels = None to make a plot without coordinate axes plotted. - + The axes can be marked in either sexagesimal or decimal celestial coordinates. If RATickSteps - or decTickSteps are set to "auto", the appropriate axis scales will be determined automatically - from the size of the image array and associated WCS. The tick step sizes can be overidden. - If the coordinate axes are in sexagesimal format a dictionary in the format {'deg', 'unit'} is - needed (see L{RA_TICK_STEPS} and L{DEC_TICK_STEPS} for examples). If the coordinate axes are in + or decTickSteps are set to "auto", the appropriate axis scales will be determined automatically + from the size of the image array and associated WCS. The tick step sizes can be overidden. + If the coordinate axes are in sexagesimal format a dictionary in the format {'deg', 'unit'} is + needed (see RA_TICK_STEPS and DEC_TICK_STEPS for examples). If the coordinate axes are in decimal format, the tick step size is specified simply in RA, dec decimal degrees. - - @type imageData: np array or list - @param imageData: image data array or list of np arrays [r, g, b] - @type imageWCS: astWCS.WCS - @param imageWCS: astWCS.WCS object - @type axes: list - @param axes: specifies where in the current figure to draw the finder chart (see pylab.axes) - @type cutLevels: list - @param cutLevels: sets the image scaling - available options: - - pixel values: cutLevels=[low value, high value]. - - histogram equalisation: cutLevels=["histEq", number of bins ( e.g. 1024)] - - relative: cutLevels=["relative", cut per cent level (e.g. 99.5)] - - smart: cutLevels=["smart", cut per cent level (e.g. 99.5)] - ["smart", 99.5] seems to provide good scaling over a range of different images. - Note that for RGB images, cut levels must be specified manually i.e. as a list: - [[r min, rmax], [g min, g max], [b min, b max]] - @type colorMapName: string - @param colorMapName: name of a standard matplotlib colormap, e.g. "hot", "cool", "gray" - etc. (do "help(pylab.colormaps)" in the Python interpreter to see available options) - @type title: string - @param title: optional title for the plot - @type axesLabels: string - @param axesLabels: either "sexagesimal" (for H:M:S, D:M:S), "decimal" (for decimal degrees) - or None (for no coordinate axes labels) - @type axesFontFamily: string - @param axesFontFamily: matplotlib fontfamily, e.g. 'serif', 'sans-serif' etc. - @type axesFontSize: float - @param axesFontSize: font size of axes labels and titles (in points) - @type colorBar: bool - @param colorBar: if True, plot a vertical color bar at the side of the image indicating the intensity - scale. - @type interpolation: string - @param interpolation: interpolation to apply to the image plot (see the documentation for - the matplotlib.pylab.imshow command) - + + Args: + imageData (numpy.ndarray or list): image data array or list of np arrays [r, g, b] + imageWCS (astWCS.WCS): astWCS.WCS object + axes (list): specifies where in the current figure to draw the finder chart (see pylab.axes) + cutLevels (list): sets the image scaling - available options: + + - pixel values: ``cutLevels=[low value, high value]`` + - histogram equalisation: ``cutLevels=["histEq", number of bins (e.g. 1024)]`` + - relative: ``cutLevels=["relative", cut per cent level (e.g. 99.5)]`` + - smart: ``cutLevels=["smart", cut per cent level (e.g. 99.5)]`` + + ``["smart", 99.5]`` seems to provide good scaling over a range of different images. + For RGB images, cut levels must be specified manually as a list: + ``[[r min, r max], [g min, g max], [b min, b max]]`` + colorMapName (str): name of a standard matplotlib colormap, e.g. "hot", "cool", "gray" + title (str): optional title for the plot + axesLabels (str): either "sexagesimal" (for H:M:S, D:M:S), "decimal" (for decimal degrees) + or None (for no coordinate axes labels) + axesFontFamily (str): matplotlib fontfamily, e.g. 'serif', 'sans-serif' etc. + axesFontSize (float): font size of axes labels and titles (in points) + colorBar (bool): if True, plot a vertical color bar at the side of the image indicating + the intensity scale + interpolation (str): interpolation to apply to the image plot (see the documentation for + the matplotlib.pylab.imshow command) + """ - + self.RADeg, self.decDeg=imageWCS.getCentreWCSCoords() self.wcs=imageWCS - + # Handle case where imageData is [r, g, b] if type(imageData) == list: if len(imageData) == 3: @@ -178,7 +169,7 @@ def __init__(self, imageData, imageWCS, axes = [0.1,0.1,0.8,0.8], \ else: self.data=imageData self.rgbImage=False - + self.axes=pylab.axes(axes) self.cutLevels=cutLevels self.colorMapName=colorMapName @@ -187,14 +178,14 @@ def __init__(self, imageData, imageWCS, axes = [0.1,0.1,0.8,0.8], \ self.colorBar=colorBar self.axesFontSize=axesFontSize self.axesFontFamily=axesFontFamily - + self.flipXAxis=False self.flipYAxis=False - + self.interpolation=interpolation - + if self.axesLabels != None: - + # Allow user to override the automatic coord tick spacing if self.axesLabels == "sexagesimal": if RATickSteps != "auto": @@ -214,33 +205,33 @@ def __init__(self, imageData, imageWCS, axes = [0.1,0.1,0.8,0.8], \ raise Exception("decTickSteps needs to be a float (if not 'auto') for decimal axes labels") self.RATickSteps=RATickSteps self.decTickSteps=decTickSteps - + self.calcWCSAxisLabels(axesLabels = self.axesLabels) - + # this list stores objects to overplot, add to it using addPlotObjects() - self.plotObjects=[] - + self.plotObjects=[] + # this list stores image data to overlay as contours, add to it using addContourOverlay() self.contourOverlays=[] - + self.draw() def draw(self): """Redraws the ImagePlot. - + """ - + pylab.axes(self.axes) pylab.cla() - + if self.title != None: pylab.title(self.title) try: colorMap=pylab.cm.get_cmap(self.colorMapName) except AssertionError: raise Exception(self.colorMapName+"is not a defined matplotlib colormap.") - + if self.rgbImage == False: self.cutImage=astImages.intensityCutImage(self.data, self.cutLevels) if self.cutLevels[0]=="histEq": @@ -250,14 +241,14 @@ def draw(self): origin='lower', cmap=colorMap) else: pylab.imshow(self.data, interpolation="bilinear", origin='lower') - + if self.colorBar == True: pylab.colorbar(shrink=0.8) - + for c in self.contourOverlays: - pylab.contour(c['contourData']['scaledImage'], c['contourData']['contourLevels'], + pylab.contour(c['contourData']['scaledImage'], c['contourData']['contourLevels'], colors=c['color'], linewidths=c['width']) - + for p in self.plotObjects: for x, y, l in zip(p['x'], p['y'], p['objLabels']): if p['symbol'] == "circle": @@ -265,43 +256,43 @@ def draw(self): linewidth=p['width']) self.axes.add_patch(c) elif p['symbol'] == "box": - c=patches.Rectangle((x-p['sizePix']/2, y-p['sizePix']/2), p['sizePix'], p['sizePix'], + c=patches.Rectangle((x-p['sizePix']/2, y-p['sizePix']/2), p['sizePix'], p['sizePix'], fill=p['fill'], color=p['color'], linewidth=p['width']) self.axes.add_patch(c) elif p['symbol'] == "cross": - pylab.plot([x-p['sizePix']/2, x+p['sizePix']/2], [y, y], linestyle='-', + pylab.plot([x-p['sizePix']/2, x+p['sizePix']/2], [y, y], linestyle='-', linewidth=p['width'], color= p['color']) - pylab.plot([x, x], [y-p['sizePix']/2, y+p['sizePix']/2], linestyle='-', + pylab.plot([x, x], [y-p['sizePix']/2, y+p['sizePix']/2], linestyle='-', linewidth=p['width'], color= p['color']) elif p['symbol'] == "diamond": - c=patches.RegularPolygon([x, y], 4, radius=p['sizePix']/2, orientation=0, + c=patches.RegularPolygon([x, y], 4, radius=p['sizePix']/2, orientation=0, color=p['color'], fill=p['fill'], linewidth=p['width']) self.axes.add_patch(c) if l != None: pylab.text(x, y+p['sizePix']/1.5, l, horizontalalignment='center', \ fontsize=p['objLabelSize'], color=p['color']) - + if p['symbol'] == "compass": x=p['x'][0] y=p['y'][0] ra=p['RA'][0] dec=p['dec'][0] - + westPoint,eastPoint,southPoint,northPoint=astCoords.calcRADecSearchBox(ra, dec, p['sizeArcSec']/3600.0/2.0) northPix=self.wcs.wcs2pix(ra, northPoint) eastPix=self.wcs.wcs2pix(eastPoint, dec) - + edx=eastPix[0]-x edy=eastPix[1]-y ndx=northPix[0]-x ndy=northPix[1]-y nArrow=patches.Arrow(x, y, ndx, ndy, edgecolor=p['color'], facecolor=p['color'], width=p['width']) - eArrow=patches.Arrow(x, y, edx, edy, edgecolor=p['color'], facecolor=p['color'], width=p['width']) + eArrow=patches.Arrow(x, y, edx, edy, edgecolor=p['color'], facecolor=p['color'], width=p['width']) self.axes.add_patch(nArrow) self.axes.add_patch(eArrow) - pylab.text(x+ndx+ndx*0.2, y+ndy+ndy*0.2, "N", horizontalalignment='center', + pylab.text(x+ndx+ndx*0.2, y+ndy+ndy*0.2, "N", horizontalalignment='center', verticalalignment='center', fontsize=p['objLabelSize'], color=p['color']) - pylab.text(x+edx+edx*0.2, y+edy+edy*0.2, "E", horizontalalignment='center', + pylab.text(x+edx+edx*0.2, y+edy+edy*0.2, "E", horizontalalignment='center', verticalalignment='center', fontsize=p['objLabelSize'], color=p['color']) if p['symbol'] == "scaleBar": @@ -309,7 +300,7 @@ def draw(self): y=p['y'][0] ra=p['RA'][0] dec=p['dec'][0] - + westPoint,eastPoint,southPoint,northPoint=astCoords.calcRADecSearchBox(ra, dec, p['sizeArcSec']/3600.0/2.0) northPix=self.wcs.wcs2pix(ra, northPoint) eastPix=self.wcs.wcs2pix(eastPoint, dec) @@ -317,16 +308,16 @@ def draw(self): edy=eastPix[1]-y ndx=northPix[0]-x ndy=northPix[1]-y - + if p['style'] == "arrows": - eArrow=patches.Arrow(x, y, edx, edy, edgecolor=p['color'], facecolor=p['color'], width=p['width']) - wArrow=patches.Arrow(x, y, -edx, edy, edgecolor=p['color'], facecolor=p['color'], width=p['width']) + eArrow=patches.Arrow(x, y, edx, edy, edgecolor=p['color'], facecolor=p['color'], width=p['width']) + wArrow=patches.Arrow(x, y, -edx, edy, edgecolor=p['color'], facecolor=p['color'], width=p['width']) self.axes.add_patch(eArrow) self.axes.add_patch(wArrow) elif p['style'] == "whiskers": - ewArrow=patches.FancyArrowPatch(posA = (x+edx, y), posB = (x-edx,y+edy), edgecolor=p['color'], facecolor=p['color'], linewidth = p['width'], arrowstyle = '|-|') + ewArrow=patches.FancyArrowPatch(posA = (x+edx, y), posB = (x-edx,y+edy), edgecolor=p['color'], facecolor=p['color'], linewidth = p['width'], arrowstyle = '|-|') self.axes.add_patch(ewArrow) - + # Work out label if p['scaleBarLabel'] == None: scaleLabel=None @@ -338,9 +329,9 @@ def draw(self): scaleLabel="%.0f %s" % (p['sizeArcSec']/3600.0, DEG) else: scaleLabel=p['scaleBarLabel'] - pylab.text(x, y+0.025*self.data.shape[1], scaleLabel, horizontalalignment='center', + pylab.text(x, y+0.025*self.data.shape[1], scaleLabel, horizontalalignment='center', verticalalignment='center', fontsize=p['objLabelSize'], color=p['color']) - + if self.axesLabels != None: pylab.xticks(self.ticsRA[0], self.ticsRA[1], weight='normal', family=self.axesFontFamily, \ fontsize=self.axesFontSize) @@ -353,7 +344,7 @@ def draw(self): pylab.yticks([], []) pylab.xlabel("") pylab.ylabel("") - + if self.flipXAxis == False: pylab.xlim(0, self.data.shape[1]-1) else: @@ -366,44 +357,38 @@ def draw(self): def addContourOverlay(self, contourImageData, contourWCS, tag, levels = ["linear", "min", "max", 5], width = 1, color = "white", smooth = 0, highAccuracy = False): - """Adds image data to the ImagePlot as a contour overlay. The contours can be removed using - L{removeContourOverlay}. If a contour overlay already exists with this tag, it will be replaced. - - @type contourImageData: np array - @param contourImageData: image data array from which contours are to be generated - @type contourWCS: astWCS.WCS - @param contourWCS: astWCS.WCS object for the image to be contoured - @type tag: string - @param tag: identifying tag for this set of contours - @type levels: list - @param levels: sets the contour levels - available options: - - values: contourLevels=[list of values specifying each level] - - linear spacing: contourLevels=['linear', min level value, max level value, number - of levels] - can use "min", "max" to automatically set min, max levels from image data - - log spacing: contourLevels=['log', min level value, max level value, number of - levels] - can use "min", "max" to automatically set min, max levels from image data - @type width: int - @param width: width of the overlaid contours - @type color: string - @param color: color of the overlaid contours, specified by the name of a standard - matplotlib color, e.g., "black", "white", "cyan" - etc. (do "help(pylab.colors)" in the Python interpreter to see available options) - @type smooth: float - @param smooth: standard deviation (in arcsec) of Gaussian filter for - pre-smoothing of contour image data (set to 0 for no smoothing) - @type highAccuracy: bool - @param highAccuracy: if True, sample every corresponding pixel in each image; otherwise, sample - every nth pixel, where n = the ratio of the image scales. - + """Adds image data to the ImagePlot as a contour overlay. The contours can be removed using + removeContourOverlay. If a contour overlay already exists with this tag, it will be replaced. + + Args: + contourImageData (numpy.ndarray): image data array from which contours are to be generated + contourWCS (astWCS.WCS): astWCS.WCS object for the image to be contoured + tag (str): identifying tag for this set of contours + levels (list): sets the contour levels - available options: + + - values: ``contourLevels=[list of values specifying each level]`` + - linear spacing: ``contourLevels=['linear', min level value, max level value, number of levels]`` + (use "min", "max" to automatically set min, max levels from image data) + - log spacing: ``contourLevels=['log', min level value, max level value, number of levels]`` + (use "min", "max" to automatically set min, max levels from image data) + + width (int): width of the overlaid contours + color (str): color of the overlaid contours, specified by the name of a standard matplotlib + color, e.g., "black", "white", "cyan" + smooth (float): standard deviation (in arcsec) of Gaussian filter for pre-smoothing of + contour image data (set to 0 for no smoothing) + highAccuracy (bool): if True, sample every corresponding pixel in each image; otherwise, + sample every nth pixel, where n = the ratio of the image scales + """ - + if self.rgbImage == True: backgroundData=self.data[:,:,0] else: backgroundData=self.data contourData=astImages.generateContourOverlay(backgroundData, self.wcs, contourImageData, \ contourWCS, levels, smooth, highAccuracy = highAccuracy) - + alreadyGot=False for c in self.contourOverlays: if c['tag'] == tag: @@ -412,83 +397,74 @@ def addContourOverlay(self, contourImageData, contourWCS, tag, levels = ["linear c['color']=color c['width']=width alreadyGot=True - + if alreadyGot == False: self.contourOverlays.append({'contourData': contourData, 'tag': tag, 'color': color, \ 'width': width}) self.draw() - + def removeContourOverlay(self, tag): """Removes the contourOverlay from the ImagePlot corresponding to the tag. - - @type tag: string - @param tag: tag for contour overlay in ImagePlot.contourOverlays to be removed - + + Args: + tag (str): tag for contour overlay in ImagePlot.contourOverlays to be removed + """ - + index=0 for p in self.contourOverlays: if p['tag'] == tag: self.plotObjects.remove(self.plotObjects[index]) index=index+1 self.draw() - - + + def addPlotObjects(self, objRAs, objDecs, tag, symbol="circle", size=4.0, width=1.0, color="yellow",\ fill = False, objLabels = None, objLabelSize = 12.0): - """Add objects with RA, dec coords objRAs, objDecs to the ImagePlot. Only objects that fall within + """Add objects with RA, dec coords objRAs, objDecs to the ImagePlot. Only objects that fall within the image boundaries will be plotted. - + symbol specifies the type of symbol with which to mark the object in the image. The following values are allowed: - "circle" - "box" - "cross" - "diamond" - + size specifies the diameter in arcsec of the symbol (if plotSymbol == "circle"), or the width of the box in arcsec (if plotSymbol == "box") - + width specifies the thickness of the symbol lines in pixels - + color can be any valid matplotlib color (e.g. "red", "green", etc.) - + The objects can be removed from the plot by using removePlotObjects(), and then calling - draw(). If the ImagePlot already has a set of plotObjects with the same tag, they will be + draw(). If the ImagePlot already has a set of plotObjects with the same tag, they will be replaced. - - @type objRAs: np array or list - @param objRAs: object RA coords in decimal degrees - @type objDecs: np array or list - @param objDecs: corresponding object Dec. coords in decimal degrees - @type tag: string - @param tag: identifying tag for this set of objects - @type symbol: string - @param symbol: either "circle", "box", "cross", or "diamond" - @type size: float - @param size: size of symbols to plot (radius in arcsec, or width of box) - @type width: float - @param width: width of symbols in pixels - @type color: string - @param color: any valid matplotlib color string, e.g. "red", "green" etc. - @type fill: bool - @param color: if True, fill symbols - @type objLabels: list - @param objLabels: text labels to plot next to objects in figure - @type objLabelSize: float - @param objLabelSize: size of font used for object labels (in points) - + + Args: + objRAs (numpy.ndarray or list): object RA coords in decimal degrees + objDecs (numpy.ndarray or list): corresponding object Dec. coords in decimal degrees + tag (str): identifying tag for this set of objects + symbol (str): either "circle", "box", "cross", or "diamond" + size (float): size of symbols to plot (radius in arcsec, or width of box) + width (float): width of symbols in pixels + color (str): any valid matplotlib color string, e.g. "red", "green" etc. + fill (bool): if True, fill symbols + objLabels (list): text labels to plot next to objects in figure + objLabelSize (float): size of font used for object labels (in points) + """ - + pixCoords=self.wcs.wcs2pix(objRAs, objDecs) - + xMax=self.data.shape[1] yMax=self.data.shape[0] - + if objLabels is None: objLabels=[None]*len(objRAs) - + xInPlot=[] yInPlot=[] RAInPlot=[] @@ -501,15 +477,15 @@ def addPlotObjects(self, objRAs, objDecs, tag, symbol="circle", size=4.0, width= RAInPlot.append(r) decInPlot.append(d) labelInPlot.append(l) - + xInPlot=np.array(xInPlot) yInPlot=np.array(yInPlot) RAInPlot=np.array(RAInPlot) decInPlot=np.array(decInPlot) - + # Size of symbols in pixels in plot - converted from arcsec sizePix=(size/3600.0)/self.wcs.getPixelSizeDeg() - + alreadyGot=False for p in self.plotObjects: if p['tag'] == tag: @@ -527,54 +503,52 @@ def addPlotObjects(self, objRAs, objDecs, tag, symbol="circle", size=4.0, width= p['objLabelSize']=objLabelSize p['fill']=fill alreadyGot=True - + if alreadyGot == False: self.plotObjects.append({'x': xInPlot, 'y': yInPlot, 'RA': RAInPlot, 'dec': decInPlot, - 'tag': tag, 'objLabels': labelInPlot, 'symbol': symbol, + 'tag': tag, 'objLabels': labelInPlot, 'symbol': symbol, 'sizePix': sizePix, 'width': width, 'color': color, 'objLabelSize': objLabelSize, 'sizeArcSec': size, 'fill': fill}) self.draw() - - + + def removePlotObjects(self, tag): """Removes the plotObjects from the ImagePlot corresponding to the tag. The plot must be redrawn for the change to take effect. - - @type tag: string - @param tag: tag for set of objects in ImagePlot.plotObjects to be removed - + + Args: + tag (str): tag for set of objects in ImagePlot.plotObjects to be removed + """ - + index=0 for p in self.plotObjects: if p['tag'] == tag: self.plotObjects.remove(self.plotObjects[index]) index=index+1 self.draw() - - + + def addCompass(self, location, sizeArcSec, color = "white", fontSize = 12, \ width = 20.0): - """Adds a compass to the ImagePlot at the given location ('N', 'NE', 'E', 'SE', 'S', - 'SW', 'W', or 'NW'). Note these aren't directions on the WCS coordinate grid, they are - relative positions on the plot - so N is top centre, NE is top right, SW is bottom right etc.. + """Adds a compass to the ImagePlot at the given location ('N', 'NE', 'E', 'SE', 'S', + 'SW', 'W', or 'NW'). Note these aren't directions on the WCS coordinate grid, they are + relative positions on the plot - so N is top centre, NE is top right, SW is bottom right etc.. Alternatively, pixel coordinates (x, y) in the image can be given. - - @type location: string or tuple - @param location: location in the plot where the compass is drawn: - - string: N, NE, E, SE, S, SW, W or NW - - tuple: (x, y) - @type sizeArcSec: float - @param sizeArcSec: length of the compass arrows on the plot in arc seconds - @type color: string - @param color: any valid matplotlib color string - @type fontSize: float - @param fontSize: size of font used to label N and E, in points - @type width: float - @param width: width of arrows used to mark compass - + + Args: + location (str or tuple): location in the plot where the compass is drawn: + + - string: N, NE, E, SE, S, SW, W or NW + - tuple: (x, y) + + sizeArcSec (float): length of the compass arrows on the plot in arc seconds + color (str): any valid matplotlib color string + fontSize (float): size of font used to label N and E, in points + width (float): width of arrows used to mark compass + """ - + if type(location) == str: cRADeg, cDecDeg=self.wcs.getCentreWCSCoords() RAMin, RAMax, decMin, decMax=self.wcs.getImageMinMaxWCSCoords() @@ -591,7 +565,7 @@ def addCompass(self, location, sizeArcSec, color = "white", fontSize = 12, \ foundLocation=False x=cy y=cx - if self.wcs.isFlipped() == False: + if self.wcs.isFlipped() == False: if location.find("N") != -1: y=Y-2*yBufferPix foundLocation=True @@ -625,7 +599,7 @@ def addCompass(self, location, sizeArcSec, color = "white", fontSize = 12, \ RADeg, decDeg=self.wcs.pix2wcs(x, y) else: raise Exception("didn't understand location for scale bar - should be string or tuple.") - + alreadyGot=False for p in self.plotObjects: if p['tag'] == "compass": @@ -641,10 +615,10 @@ def addCompass(self, location, sizeArcSec, color = "white", fontSize = 12, \ p['color']=color p['objLabelSize']=fontSize alreadyGot=True - + if alreadyGot == False: self.plotObjects.append({'x': [x], 'y': [y], 'RA': [RADeg], 'dec': [decDeg], - 'tag': "compass", 'objLabels': [None], 'symbol': "compass", + 'tag': "compass", 'objLabels': [None], 'symbol': "compass", 'width': width, 'color': color, 'objLabelSize': fontSize, 'sizeArcSec': sizeArcSec}) self.draw() @@ -652,30 +626,26 @@ def addCompass(self, location, sizeArcSec, color = "white", fontSize = 12, \ def addScaleBar(self, location, sizeArcSec, color = "white", fontSize = 12, \ width = 20.0, label = None, style = "whiskers"): - """Adds a scale bar to the ImagePlot at the given location ('N', 'NE', 'E', 'SE', 'S', - 'SW', 'W', or 'NW'). Note these aren't directions on the WCS coordinate grid, they are - relative positions on the plot - so N is top centre, NE is top right, SW is bottom right etc.. + """Adds a scale bar to the ImagePlot at the given location ('N', 'NE', 'E', 'SE', 'S', + 'SW', 'W', or 'NW'). Note these aren't directions on the WCS coordinate grid, they are + relative positions on the plot - so N is top centre, NE is top right, SW is bottom right etc.. Alternatively, pixel coordinates (x, y) in the image can be given. - - @type location: string or tuple - @param location: location in the plot where the compass is drawn: - - string: N, NE, E, SE, S, SW, W or NW - - tuple: (x, y) - @type sizeArcSec: float - @param sizeArcSec: scale length to indicate on the plot in arc seconds - @type color: string - @param color: any valid matplotlib color string - @type fontSize: float - @param fontSize: size of font used to label N and E, in points - @type width: float - @param width: width of arrow used to mark scale - @type label: string - @param label: overrides the displayed label if not None (if None, label is the angular size) - @type style: string - @param style: either "whiskers" or "arrows" - + + Args: + location (str or tuple): location in the plot where the scale bar is drawn: + + - string: N, NE, E, SE, S, SW, W or NW + - tuple: (x, y) + + sizeArcSec (float): scale length to indicate on the plot in arc seconds + color (str): any valid matplotlib color string + fontSize (float): size of font used to label the scale bar, in points + width (float): width of arrow used to mark scale + label (str): overrides the displayed label if not None (if None, label is the angular size) + style (str): either "whiskers" or "arrows" + """ - + # Work out where the scale bar is going in WCS coords from the relative location given if type(location) == str: cRADeg, cDecDeg=self.wcs.getCentreWCSCoords() @@ -727,7 +697,7 @@ def addScaleBar(self, location, sizeArcSec, color = "white", fontSize = 12, \ RADeg, decDeg=self.wcs.pix2wcs(x, y) else: raise Exception("didn't understand location for scale bar - should be string or tuple.") - + alreadyGot=False for p in self.plotObjects: if p['tag'] == "scaleBar": @@ -745,46 +715,46 @@ def addScaleBar(self, location, sizeArcSec, color = "white", fontSize = 12, \ p['scaleBarLabel']=label p['style']=style alreadyGot=True - + if alreadyGot == False: self.plotObjects.append({'x': [x], 'y': [y], 'RA': [RADeg], 'dec': [decDeg], - 'tag': "scaleBar", 'objLabels': [None], 'symbol': "scaleBar", + 'tag': "scaleBar", 'objLabels': [None], 'symbol': "scaleBar", 'width': width, 'color': color, 'objLabelSize': fontSize, 'sizeArcSec': sizeArcSec, 'scaleBarLabel': label, 'style': style}) self.draw() - + def calcWCSAxisLabels(self, axesLabels = "decimal"): - """This function calculates the positions of coordinate labels for the RA and Dec axes of the + """This function calculates the positions of coordinate labels for the RA and Dec axes of the ImagePlot. The tick steps are calculated automatically unless self.RATickSteps, - self.decTickSteps are set to values other than "auto" (see L{ImagePlot.__init__}). - + self.decTickSteps are set to values other than "auto" (see __init__). + The ImagePlot must be redrawn for changes to be applied. - - @type axesLabels: string - @param axesLabels: either "sexagesimal" (for H:M:S, D:M:S), "decimal" (for decimal degrees), - or None for no coordinate axes labels - + + Args: + axesLabels (str): either "sexagesimal" (for H:M:S, D:M:S), "decimal" (for decimal degrees), + or None for no coordinate axes labels + """ - + # Label equinox on axes equinox=self.wcs.getEquinox() if equinox<1984: equinoxLabel="B"+str(int(equinox)) else: equinoxLabel="J"+str(int(equinox)) - + self.axesLabels=axesLabels - + ticsDict=self.getTickSteps() - + # Manual override - note: no minor tick marks anymore, but may want to bring them back if self.RATickSteps != "auto": ticsDict['major']['RA']=self.RATickSteps if self.decTickSteps != "auto": ticsDict['major']['dec']=self.decTickSteps - + RALocs=[] decLocs=[] RALabels=[] @@ -803,7 +773,7 @@ def calcWCSAxisLabels(self, axesLabels = "decimal"): decDegStep=ticsDict[key]['dec'] else: raise Exception("axesLabels must be either 'sexagesimal' or 'decimal'") - + # xArray=np.arange(0, self.data.shape[1], 1) xArray=np.linspace(0, self.data.shape[1], self.data.shape[1]*2) yArray=np.arange(0, self.data.shape[0], 1) @@ -821,14 +791,14 @@ def calcWCSAxisLabels(self, axesLabels = "decimal"): RAMax=ras.max() decMin=decs.min() decMax=decs.max() - + # Work out if wrapped around midRAPix, midDecPix=self.wcs.wcs2pix((RAEdges[1]+RAEdges[0])/2.0, (decMax+decMin)/2.0) if midRAPix < 0 or midRAPix > self.wcs.header['NAXIS1']: wrappedRA=True else: wrappedRA=False - + # Note RA, dec work in opposite sense below because E at left if ras[1] < ras[0]: self.flipXAxis=False @@ -842,7 +812,7 @@ def calcWCSAxisLabels(self, axesLabels = "decimal"): else: self.flipYAxis=False dec2y=interpolate.interp1d(decs, yArray, kind='linear') - + if wrappedRA == False: RAPlotMin=RADegStep*math.modf(RAMin/RADegStep)[1] RAPlotMax=RADegStep*math.modf(RAMax/RADegStep)[1] @@ -878,7 +848,7 @@ def calcWCSAxisLabels(self, axesLabels = "decimal"): if decPlotMax >= decMax: decPlotMax=decPlotMax-decDegStep decDegs=np.arange(decPlotMin, decPlotMax+0.0001, decDegStep) - + if key == "major": if axesLabels == "sexagesimal": for r in RADegs: @@ -941,10 +911,10 @@ def calcWCSAxisLabels(self, axesLabels = "decimal"): elif ticsDict[key]['dec']['unit'] == 'm': dString=dString+DEG+mString+PRIME else: - dString=dString+DEG+mString+PRIME+sString+DOUBLE_PRIME + dString=dString+DEG+mString+PRIME+sString+DOUBLE_PRIME decLabels.append(dString) elif axesLabels == "decimal": - + if wrappedRA == False: RALabels=RALabels+RADegs.tolist() else: @@ -955,7 +925,7 @@ def calcWCSAxisLabels(self, axesLabels = "decimal"): nonNegativeLabels.append(r) RALabels=RALabels+nonNegativeLabels decLabels=decLabels+decDegs.tolist() - + # Format RALabels, decLabels to same number of d.p. dpNumRA=len(str(ticsDict['major']['RA']).split(".")[-1]) dpNumDec=len(str(ticsDict['major']['dec']).split(".")[-1]) @@ -964,8 +934,8 @@ def calcWCSAxisLabels(self, axesLabels = "decimal"): RALabels[i]=fString % (RALabels[i]) for i in range(len(decLabels)): fString="%."+str(dpNumDec)+"f" - decLabels[i]=fString % (decLabels[i]) - + decLabels[i]=fString % (decLabels[i]) + if key == 'minor': RALabels=RALabels+RADegs.shape[0]*[''] decLabels=decLabels+decDegs.shape[0]*[''] @@ -986,32 +956,33 @@ def calcWCSAxisLabels(self, axesLabels = "decimal"): self.ticsRA=[RALocs, RALabels] self.ticsDec=[decLocs, decLabels] - + def save(self, fileName): - """Saves the ImagePlot in any format that matplotlib can understand, as determined from the + """Saves the ImagePlot in any format that matplotlib can understand, as determined from the fileName extension. - - @type fileName: string - @param fileName: path where plot will be written - + + Args: + fileName (str): path where plot will be written + """ - + pylab.draw() pylab.savefig(fileName) - - + + def getTickSteps(self): """Chooses the appropriate WCS coordinate tick steps for the plot based on its size. Whether the ticks are decimal or sexagesimal is set by self.axesLabels. - - Note: minor ticks not used at the moment. - - @rtype: dictionary - @return: tick step sizes for major, minor plot ticks, in format {'major', 'minor'} - + + Note: + Minor ticks not used at the moment. + + Returns: + dict: tick step sizes for major, minor plot ticks, in format ``{'major', 'minor'}`` + """ - + # Aim for 5 major tick marks on a plot xArray=np.arange(0, self.data.shape[1], 1) yArray=np.arange(0, self.data.shape[0], 1) @@ -1026,7 +997,7 @@ def getTickSteps(self): RAMax=RAEdges.max() decMin=decs.min() decMax=decs.max() - + # Work out if wrapped around midRAPix, midDecPix=self.wcs.wcs2pix((RAEdges[1]+RAEdges[0])/2.0, (decMax+decMin)/2.0) if midRAPix < 0 or midRAPix > self.wcs.header['NAXIS1']: @@ -1043,12 +1014,12 @@ def getTickSteps(self): ticsDict['major']={} ticsDict['minor']={} if self.axesLabels == "sexagesimal": - + matchIndex = 0 for i in range(len(RA_TICK_STEPS)): if RAWidthDeg/2.5 > RA_TICK_STEPS[i]['deg']: matchIndex = i - + ticsDict['major']['RA']=RA_TICK_STEPS[matchIndex] ticsDict['minor']['RA']=RA_TICK_STEPS[matchIndex-1] @@ -1056,32 +1027,31 @@ def getTickSteps(self): for i in range(len(DEC_TICK_STEPS)): if decHeightDeg/2.5 > DEC_TICK_STEPS[i]['deg']: matchIndex = i - + ticsDict['major']['dec']=DEC_TICK_STEPS[matchIndex] ticsDict['minor']['dec']=DEC_TICK_STEPS[matchIndex-1] - + return ticsDict - + elif self.axesLabels == "decimal": - + matchIndex = 0 for i in range(len(DECIMAL_TICK_STEPS)): if RAWidthDeg/2.5 > DECIMAL_TICK_STEPS[i]: matchIndex = i - + ticsDict['major']['RA']=DECIMAL_TICK_STEPS[matchIndex] ticsDict['minor']['RA']=DECIMAL_TICK_STEPS[matchIndex-1] - + matchIndex = 0 for i in range(len(DECIMAL_TICK_STEPS)): if decHeightDeg/2.5 > DECIMAL_TICK_STEPS[i]: matchIndex = i - + ticsDict['major']['dec']=DECIMAL_TICK_STEPS[matchIndex] ticsDict['minor']['dec']=DECIMAL_TICK_STEPS[matchIndex-1] - + return ticsDict - + else: raise Exception("axesLabels must be either 'sexagesimal' or 'decimal'") - diff --git a/astLib/astSED.py b/astLib/astSED.py index 2889aac..f5afd29 100644 --- a/astLib/astSED.py +++ b/astLib/astSED.py @@ -1,10 +1,10 @@ """Module for performing calculations on Spectral Energy Distributions (SEDs). -(c) 2007-2013 Matt Hilton +(c) 2007-2013 Matt Hilton This module provides classes for manipulating SEDs, in particular the Bruzual & Charlot 2003, Maraston -2005, and Percival et al 2009 stellar population synthesis models are currently supported. Functions are -provided for calculating the evolution of colours and magnitudes in these models with redshift etc., and +2005, and Percival et al 2009 stellar population synthesis models are currently supported. Functions are +provided for calculating the evolution of colours and magnitudes in these models with redshift etc., and for fitting broadband photometry using these models. """ @@ -35,35 +35,35 @@ class Passband: """This class describes a filter transmission curve. Passband objects are created by loading data from from text files containing wavelength in angstroms in the first column, relative transmission efficiency - in the second column (whitespace delimited). For example, to create a Passband object for the 2MASS J + in the second column (whitespace delimited). For example, to create a Passband object for the 2MASS J filter: - + passband=astSED.Passband("J_2MASS.res") - + where "J_2MASS.res" is a file in the current working directory that describes the filter. - + Wavelength units can be specified as 'angstroms', 'nanometres' or 'microns'; if either of the latter, they will be converted to angstroms. - + """ def __init__(self, fileName, normalise = True, inputUnits = 'angstroms', wavelengthColumn = 0, transmissionColumn = 1): - + inFile=open(fileName, "r") lines=inFile.readlines() - + wavelength=[] transmission=[] for line in lines: - + if line[0] != "#" and len(line) > 3: - + bits=line.split() transmission.append(float(bits[transmissionColumn])) wavelength.append(float(bits[wavelengthColumn])) - + self.wavelength=np.array(wavelength) self.transmission=np.array(transmission) - + if inputUnits == 'angstroms': pass elif inputUnits == 'nanometres': @@ -77,68 +77,66 @@ def __init__(self, fileName, normalise = True, inputUnits = 'angstroms', wavelen self.wavelength=self.wavelength*1e10 else: raise Exception("didn't understand passband input units") - + # Sort into ascending order of wavelength otherwise normalisation will be wrong merged=np.array([self.wavelength, self.transmission]).transpose() sortedMerged=np.array(sorted(merged, key=operator.itemgetter(0))) self.wavelength=sortedMerged[:, 0] self.transmission=sortedMerged[:, 1] - + if normalise == True: self.transmission=self.transmission/np.trapezoid(self.transmission, self.wavelength) - + # Store a ready-to-go interpolation object to speed calculation of fluxes up self.interpolator=interpolate.interp1d(self.wavelength, self.transmission, kind='linear') def asList(self): """Returns a two dimensional list of [wavelength, transmission], suitable for plotting by gnuplot. - - @rtype: list - @return: list in format [wavelength, transmission] - + + Returns: + list: list in format [wavelength, transmission] + """ - + listData=[] for l, f in zip(self.wavelength, self.transmission): listData.append([l, f]) - + return listData - + def rescale(self, maxTransmission): """Rescales the passband so that maximum value of the transmission is equal to maxTransmission. Useful for plotting. - - @type maxTransmission: float - @param maxTransmission: maximum value of rescaled transmission curve - + + Args: + maxTransmission (float): maximum value of rescaled transmission curve + """ - + self.transmission=self.transmission*(maxTransmission/self.transmission.max()) def plot(self, xmin = 'min', xmax = 'max', maxTransmission = None): - """Plots the passband, rescaling the maximum of the tranmission curve to maxTransmission if + """Plots the passband, rescaling the maximum of the transmission curve to maxTransmission if required. - - @type xmin: float or 'min' - @param xmin: minimum of the wavelength range of the plot - @type xmax: float or 'max' - @param xmax: maximum of the wavelength range of the plot - @type maxTransmission: float - @param maxTransmission: maximum value of rescaled transmission curve - + + Args: + xmin (float or 'min'): minimum of the wavelength range of the plot + xmax (float or 'max'): maximum of the wavelength range of the plot + maxTransmission (float): maximum value of rescaled transmission curve + """ - + if maxTransmission != None: self.rescale(maxTransmission) - + pylab.matplotlib.interactive(True) pylab.plot(self.wavelength, self.transmission) - + if xmin == 'min': xmin=self.wavelength.min() if xmax == 'max': xmax=self.wavelength.max() - + pylab.xlim(xmin, xmax) pylab.xlabel("Wavelength") pylab.ylabel("Relative Flux") @@ -146,63 +144,61 @@ def plot(self, xmin = 'min', xmax = 'max', maxTransmission = None): def effectiveWavelength(self): """Calculates effective wavelength for the passband. This is the same as equation (3) of Carter et al. 2009. - - @rtype: float - @return: effective wavelength of the passband, in Angstroms - + + Returns: + float: effective wavelength of the passband, in Angstroms + """ - + a=np.trapezoid(self.transmission*self.wavelength, self.wavelength) b=np.trapezoid(self.transmission/self.wavelength, self.wavelength) effWavelength=np.sqrt(a/b) - + return effWavelength #------------------------------------------------------------------------------------------------------------ class TopHatPassband(Passband): """This class generates a passband with a top hat response between the given wavelengths. - + """ - + def __init__(self, wavelengthMin, wavelengthMax, normalise = True): """Generates a passband object with top hat response between wavelengthMin, wavelengthMax. Units are assumed to be Angstroms. - - @type wavelengthMin: float - @param wavelengthMin: minimum of the wavelength range of the passband - @type wavelengthMax: float - @param wavelengthMax: maximum of the wavelength range of the passband - @type normalise: bool - @param normalise: if True, scale such that total area under the passband over the wavelength - range is 1. - + + Args: + wavelengthMin (float): minimum of the wavelength range of the passband + wavelengthMax (float): maximum of the wavelength range of the passband + normalise (bool): if True, scale such that total area under the passband over the wavelength + range is 1 + """ - + self.wavelength=np.arange(wavelengthMin, wavelengthMax+10, 10, dtype = float) self.transmission=np.ones(self.wavelength.shape, dtype = float) - + if normalise == True: self.transmission=self.transmission/np.trapezoid(self.transmission, self.wavelength) - + # Store a ready-to-go interpolation object to speed calculation of fluxes up self.interpolator=interpolate.interp1d(self.wavelength, self.transmission, kind='linear') - - + + #------------------------------------------------------------------------------------------------------------ class SED: """This class describes a Spectral Energy Distribution (SED). - + To create a SED object, lists (or np arrays) of wavelength and relative flux must be provided. The SED - can optionally be redshifted. The wavelength units of SEDs are assumed to be Angstroms - flux + can optionally be redshifted. The wavelength units of SEDs are assumed to be Angstroms - flux calculations using Passband and SED objects specified with different wavelength units will be incorrect. - - The L{StellarPopulation} class (and derivatives) can be used to extract SEDs for specified ages from e.g. + + The StellarPopulation class (and derivatives) can be used to extract SEDs for specified ages from e.g. the Bruzual & Charlot 2003 or Maraston 2005 models. - + """ - + def __init__(self, wavelength = [], flux = [], z = 0.0, ageGyr = None, normalise = False, label = None): - + # We keep a copy of the wavelength, flux at z = 0, as it's more robust to copy that # to self.wavelength, flux and redshift it, rather than repeatedly redshifting the same # arrays back and forth @@ -212,14 +208,14 @@ def __init__(self, wavelength = [], flux = [], z = 0.0, ageGyr = None, normalise self.flux=np.array(flux) self.z=z self.label=label # plain text label, handy for using in photo-z codes - + # Store the intrinsic (i.e. unextincted) flux in case we change extinction self.EBMinusV=0.0 self.intrinsic_z0flux=np.array(flux) - + if normalise == True: self.normalise() - + if z != 0.0: self.redshift(z) @@ -227,28 +223,28 @@ def __init__(self, wavelength = [], flux = [], z = 0.0, ageGyr = None, normalise def copy(self): - """Copies the SED, returning a new SED object - - @rtype: L{SED} object - @return: SED - + """Copies the SED, returning a new SED object. + + Returns: + SED: copy of the SED + """ - - newSED=SED(wavelength = self.z0wavelength, flux = self.z0flux, z = self.z, ageGyr = self.ageGyr, + + newSED=SED(wavelength = self.z0wavelength, flux = self.z0flux, z = self.z, ageGyr = self.ageGyr, normalise = False, label = self.label) - + return newSED - - + + def loadFromFile(self, fileName): """Loads SED from a white space delimited file in the format wavelength, flux. Lines beginning with # are ignored. - - @type fileName: string - @param fileName: path to file containing wavelength, flux data - + + Args: + fileName (str): path to file containing wavelength, flux data + """ - + inFile=open(fileName, "r") lines=inFile.readlines() inFile.close() @@ -260,12 +256,12 @@ def loadFromFile(self, fileName): bits=line.split() wavelength.append(float(bits[0])) flux.append(float(bits[1])) - + # Sort SED so wavelength is in ascending order if wavelength[0] > wavelength[-1]: wavelength.reverse() flux.reverse() - + self.z0wavelength=np.array(wavelength) self.z0flux=np.array(flux) self.wavelength=np.array(wavelength) @@ -275,51 +271,50 @@ def loadFromFile(self, fileName): def writeToFile(self, fileName): """Writes SED to a white space delimited file in the format wavelength, flux. - - @type fileName: string - @param fileName: path to file - + + Args: + fileName (str): path to file + """ - + outFile=open(fileName, "w") for l, f in zip(self.wavelength, self.flux): outFile.write(str(l)+" "+str(f)+"\n") outFile.close() - + def asList(self): """Returns a two dimensional list of [wavelength, flux], suitable for plotting by gnuplot. - - @rtype: list - @return: list in format [wavelength, flux] - + + Returns: + list: list in format [wavelength, flux] + """ - + listData=[] for l, f in zip(self.wavelength, self.flux): listData.append([l, f]) - + return listData - + def plot(self, xmin = 'min', xmax = 'max'): """Produces a simple (wavelength, flux) plot of the SED. - - @type xmin: float or 'min' - @param xmin: minimum of the wavelength range of the plot - @type xmax: float or 'max' - @param xmax: maximum of the wavelength range of the plot - + + Args: + xmin (float or 'min'): minimum of the wavelength range of the plot + xmax (float or 'max'): maximum of the wavelength range of the plot + """ - + pylab.matplotlib.interactive(True) pylab.plot(self.wavelength, self.flux) - + if xmin == 'min': xmin=self.wavelength.min() if xmax == 'max': xmax=self.wavelength.max() - + # Sensible y scale plotMask=np.logical_and(np.greater(self.wavelength, xmin), np.less(self.wavelength, xmax)) plotMax=self.flux[plotMask].max() @@ -327,51 +322,51 @@ def plot(self, xmin = 'min', xmax = 'max'): pylab.xlim(xmin, xmax) pylab.xlabel("Wavelength") pylab.ylabel("Relative Flux") - + def integrate(self, wavelengthMin = 'min', wavelengthMax = 'max'): """Calculates flux in SED within given wavelength range. - - @type wavelengthMin: float or 'min' - @param wavelengthMin: minimum of the wavelength range - @type wavelengthMax: float or 'max' - @param wavelengthMax: maximum of the wavelength range - @rtype: float - @return: relative flux - + + Args: + wavelengthMin (float or 'min'): minimum of the wavelength range + wavelengthMax (float or 'max'): maximum of the wavelength range + + Returns: + float: relative flux + """ if wavelengthMin == 'min': wavelengthMin=self.wavelength.min() if wavelengthMax == 'max': wavelengthMax=self.wavelength.max() - + mask=np.logical_and(np.greater(self.wavelength, wavelengthMin), \ np.less(self.wavelength, wavelengthMax)) flux=np.trapezoid(self.flux[mask], self.wavelength[mask]) - + return flux - + def smooth(self, smoothPix): """Smooths SED.flux with a uniform (boxcar) filter of width smoothPix. Cannot be undone. - - @type smoothPix: int - @param smoothPix: size of uniform filter applied to SED, in pixels - + + Args: + smoothPix (int): size of uniform filter applied to SED, in pixels + """ smoothed=ndimage.uniform_filter1d(self.flux, smoothPix) self.flux=smoothed - + def redshift(self, z): """Redshifts the SED to redshift z. - - @type z: float - @param z: redshift - + + Args: + z (float): redshift + """ - + # We have to conserve energy so the area under the redshifted SED has to be equal to # the area under the unredshifted SED, otherwise magnitude calculations will be wrong # when comparing SEDs at different zs @@ -379,184 +374,179 @@ def redshift(self, z): self.flux=np.zeros(self.z0flux.shape[0]) self.wavelength=self.wavelength+self.z0wavelength self.flux=self.flux+self.z0flux - + z0TotalFlux=np.trapezoid(self.z0wavelength, self.z0flux) self.wavelength=self.wavelength*(1.0+z) zTotalFlux=np.trapezoid(self.wavelength, self.flux) self.flux=self.flux*(z0TotalFlux/zTotalFlux) self.z=z - + def normalise(self, minWavelength = 'min', maxWavelength = 'max'): """Normalises the SED such that the area under the specified wavelength range is equal to 1. - - @type minWavelength: float or 'min' - @param minWavelength: minimum wavelength of range over which to normalise SED - @type maxWavelength: float or 'max' - @param maxWavelength: maximum wavelength of range over which to normalise SED - + + Args: + minWavelength (float or 'min'): minimum wavelength of range over which to normalise SED + maxWavelength (float or 'max'): maximum wavelength of range over which to normalise SED + """ if minWavelength == 'min': minWavelength=self.wavelength.min() if maxWavelength == 'max': maxWavelength=self.wavelength.max() - + lowCut=np.greater(self.wavelength, minWavelength) highCut=np.less(self.wavelength, maxWavelength) totalCut=np.logical_and(lowCut, highCut) sedFluxSlice=self.flux[totalCut] sedWavelengthSlice=self.wavelength[totalCut] - + self.flux=self.flux/np.trapezoid(abs(sedFluxSlice), sedWavelengthSlice)#self.wavelength) def normaliseToMag(self, ABMag, passband): """Normalises the SED to match the flux equivalent to the given AB magnitude in the given passband. - - @type ABMag: float - @param ABMag: AB magnitude to which the SED is to be normalised at the given passband - @type passband: an L{Passband} object - @param passband: passband at which normalisation to AB magnitude is calculated - + + Args: + ABMag (float): AB magnitude to which the SED is to be normalised at the given passband + passband (Passband): passband at which normalisation to AB magnitude is calculated + """ - + magFlux=mag2Flux(ABMag, 0.0, passband) sedFlux=self.calcFlux(passband) norm=magFlux[0]/sedFlux self.flux=self.flux*norm self.z0flux=self.z0flux*norm - + def matchFlux(self, matchSED, minWavelength, maxWavelength): """Matches the flux in the wavelength range given by minWavelength, maxWavelength to the flux in the same region in matchSED. Useful for plotting purposes. - - @type matchSED: L{SED} object - @param matchSED: SED to match flux to - @type minWavelength: float - @param minWavelength: minimum of range in which to match flux of current SED to matchSED - @type maxWavelength: float - @param maxWavelength: maximum of range in which to match flux of current SED to matchSED - + + Args: + matchSED (SED): SED to match flux to + minWavelength (float): minimum of range in which to match flux of current SED to matchSED + maxWavelength (float): maximum of range in which to match flux of current SED to matchSED + """ - + interpMatch=interpolate.interp1d(matchSED.wavelength, matchSED.flux, kind='linear') interpSelf=interpolate.interp1d(self.wavelength, self.flux, kind='linear') - + wavelengthRange=np.arange(minWavelength, maxWavelength, 5.0) - + matchFlux=np.trapezoid(interpMatch(wavelengthRange), wavelengthRange) selfFlux=np.trapezoid(interpSelf(wavelengthRange), wavelengthRange) - + self.flux=self.flux*(matchFlux/selfFlux) - + def calcFlux(self, passband): """Calculates flux in the given passband. - - @type passband: L{Passband} object - @param passband: filter passband through which to calculate the flux from the SED - @rtype: float - @return: flux - + + Args: + passband (Passband): filter passband through which to calculate the flux from the SED + + Returns: + float: flux + """ lowCut=np.greater(self.wavelength, passband.wavelength.min()) highCut=np.less(self.wavelength, passband.wavelength.max()) totalCut=np.logical_and(lowCut, highCut) sedFluxSlice=self.flux[totalCut] sedWavelengthSlice=self.wavelength[totalCut] - - # Use linear interpolation to rebin the passband to the same dimensions as the + + # Use linear interpolation to rebin the passband to the same dimensions as the # part of the SED we're interested in - sedInBand=passband.interpolator(sedWavelengthSlice)*sedFluxSlice + sedInBand=passband.interpolator(sedWavelengthSlice)*sedFluxSlice totalFlux=np.trapezoid(sedInBand*sedWavelengthSlice, sedWavelengthSlice) totalFlux=totalFlux/np.trapezoid(passband.interpolator(sedWavelengthSlice)\ *sedWavelengthSlice, sedWavelengthSlice) - - return totalFlux - + + return totalFlux + def calcMag(self, passband, addDistanceModulus = True, magType = "Vega"): """Calculates magnitude in the given passband. If addDistanceModulus == True, then the distance modulus (5.0*log10*(dl*1e5), where dl is the luminosity distance - in Mpc at the redshift of the L{SED}) is added. - - @type passband: L{Passband} object - @param passband: filter passband through which to calculate the magnitude from the SED - @type addDistanceModulus: bool - @param addDistanceModulus: if True, adds 5.0*log10*(dl*1e5) to the mag returned, where - dl is the luminosity distance (Mpc) corresponding to the SED z - @type magType: string - @param magType: either "Vega" or "AB" - @rtype: float - @return: magnitude through the given passband on the specified magnitude system - + in Mpc at the redshift of the SED) is added. + + Args: + passband (Passband): filter passband through which to calculate the magnitude from the SED + addDistanceModulus (bool): if True, adds 5.0*log10*(dl*1e5) to the mag returned, where + dl is the luminosity distance (Mpc) corresponding to the SED z + magType (str): either "Vega" or "AB" + + Returns: + float: magnitude through the given passband on the specified magnitude system + """ f1=self.calcFlux(passband) if magType == "Vega": f2=VEGA.calcFlux(passband) elif magType == "AB": f2=AB.calcFlux(passband) - + mag=-2.5*math.log10(f1/f2) if magType == "Vega": mag=mag+0.026 # Add 0.026 because Vega has V=0.026 (e.g. Bohlin & Gilliland 2004) - + if self.z > 0.0 and addDistanceModulus == True: appMag=5.0*math.log10(astCalc.dl(self.z)*1e5)+mag else: appMag=mag - + return appMag - + def calcColour(self, passband1, passband2, magType = "Vega"): """Calculates the colour passband1-passband2. - - @type passband1: L{Passband} object - @param passband1: filter passband through which to calculate the first magnitude - @type passband2: L{Passband} object - @param passband1: filter passband through which to calculate the second magnitude - @type magType: string - @param magType: either "Vega" or "AB" - @rtype: float - @return: colour defined by passband1 - passband2 on the specified magnitude system + + Args: + passband1 (Passband): filter passband through which to calculate the first magnitude + passband2 (Passband): filter passband through which to calculate the second magnitude + magType (str): either "Vega" or "AB" + + Returns: + float: colour defined by passband1 - passband2 on the specified magnitude system """ mag1=self.calcMag(passband1, magType = magType, addDistanceModulus = True) mag2=self.calcMag(passband2, magType = magType, addDistanceModulus = True) - + colour=mag1-mag2 return colour - + def getSEDDict(self, passbands): """This is a convenience function for pulling out fluxes from a SED for a given set of passbands - in the same format as made by L{mags2SEDDict} - designed to make fitting code simpler. - - @type passbands: list of L{Passband} objects - @param passbands: list of passbands through which fluxes will be calculated - + in the same format as made by mags2SEDDict - designed to make fitting code simpler. + + Args: + passbands (list): list of Passband objects through which fluxes will be calculated + """ - + flux=[] wavelength=[] for p in passbands: flux.append(self.calcFlux(p)) wavelength.append(p.effectiveWavelength()) - + SEDDict={} SEDDict['flux']=np.array(flux) SEDDict['wavelength']=np.array(wavelength) - + return SEDDict def addEmissionLines(self, scaleFactor = 1.0, modifyIntrinsicFlux = True): """Add emission lines to the SED (mostly) in the style of Ilbert+2009. - @type scaleFactor: float - @param scaleFactor: line flux will be scaled by this factor. + Args: + scaleFactor (float): line flux will be scaled by this factor """ @@ -581,21 +571,21 @@ def addEmissionLines(self, scaleFactor = 1.0, modifyIntrinsicFlux = True): self.intrinsic_z0flux=self.z0flux.copy() self.redshift(self.z) - + def extinctionCalzetti(self, EBMinusV): """Applies the Calzetti et al. 2000 (ApJ, 533, 682) extinction law to the SED with the given E(B-V) amount of extinction. R_v' = 4.05 is assumed (see equation (5) of Calzetti et al.). - - @type EBMinusV: float - @param EBMinusV: extinction E(B-V), in magnitudes - + + Args: + EBMinusV (float): extinction E(B-V), in magnitudes + """ - + self.EBMinusV=EBMinusV - + # All done in rest frame self.z0flux=self.intrinsic_z0flux.copy() - + # Allow us to set EBMinusV == 0 to turn extinction off if EBMinusV > 0: # Note that EBMinusV is assumed to be Es as in equations (2) - (5) @@ -626,7 +616,7 @@ def extinctionCalzetti(self, EBMinusV): intercept=interpolator(1200.0)-(slope*1200.0) mask=np.less(self.z0wavelength, 1200.0) kPrime[mask]=slope*self.z0wavelength[mask]+intercept - + # Long wavelengths interpolator=interpolate.interp1d(self.z0wavelength, kPrimeLong, kind='linear') slope=(interpolator(21900.0)-interpolator(22000.0))/(21900.0-22000.0) @@ -635,26 +625,26 @@ def extinctionCalzetti(self, EBMinusV): kPrime[mask]=slope*self.z0wavelength[mask]+intercept except: raise Exception("This SED has a wavelength range that doesn't cover ~1200-22000 Angstroms") - + # Never let go negative kPrime[np.less_equal(kPrime, 0.0)]=1e-6 - + reddening=np.power(10, 0.4*EBMinusV*kPrime) self.z0flux=self.z0flux/reddening self.redshift(self.z) - + #------------------------------------------------------------------------------------------------------------ class VegaSED(SED): """This class stores the SED of Vega, used for calculation of magnitudes on the Vega system. - + The Vega SED used is taken from Bohlin 2007 (http://adsabs.harvard.edu/abs/2007ASPC..364..315B), and is available from the STScI CALSPEC library (http://www.stsci.edu/hst/observatory/cdbs/calspec.html). - + """ - + def __init__(self, normalise = False): - + VEGA_SED_PATH=astLib.__path__[0]+os.path.sep+"data"+os.path.sep+"bohlin2006_Vega.sed" # from HST CALSPEC with open(VEGA_SED_PATH, "r") as inFile: @@ -663,41 +653,41 @@ def __init__(self, normalise = False): wavelength=[] flux=[] for line in lines: - + if line[0] != "#" and len(line) > 3: - + bits=line.split() flux.append(float(bits[1])) wavelength.append(float(bits[0])) - + self.wavelength=np.array(wavelength) self.flux=np.array(flux, dtype=np.float64) - + # We may want to redshift reference SEDs to calculate rest-frame colors from SEDs at different zs self.z0wavelength=np.array(wavelength) self.z0flux=np.array(flux, dtype=np.float64) self.z=0.0 - + #if normalise == True: #self.flux=self.flux/np.trapezoid(self.flux, self.wavelength) #self.z0flux=self.z0flux/np.trapezoid(self.z0flux, self.z0wavelength) - + #------------------------------------------------------------------------------------------------------------ class StellarPopulation: """This class describes a stellar population model, either a Simple Stellar Population (SSP) or a Composite Stellar Population (CSP), such as the models of Bruzual & Charlot 2003 or Maraston 2005. - + The constructor for this class can be used for generic SSPs or CSPs stored in white space delimited text files, containing columns for age, wavelength, and flux. Columns are counted from 0 ... n. Lines starting with # are ignored. - - The classes L{M05Model} (for Maraston 2005 models), L{BC03Model} (for Bruzual & Charlot 2003 models), and - L{P09Model} (for Percival et al. 2009 models) are derived from this class. The only difference between + + The classes M05Model (for Maraston 2005 models), BC03Model (for Bruzual & Charlot 2003 models), and + P09Model (for Percival et al. 2009 models) are derived from this class. The only difference between them is the code used to load in the model data. - + """ def __init__(self, fileName, ageColumn = 0, wavelengthColumn = 1, fluxColumn = 2): - + with open(fileName, "r") as inFile: lines=inFile.readlines() @@ -715,7 +705,7 @@ def __init__(self, fileName, ageColumn = 0, wavelengthColumn = 1, fluxColumn = 2 self.ages.append(age) if wavelength not in self.wavelengths: self.wavelengths.append(wavelength) - + # Construct a grid of flux - rows correspond to each wavelength, columns to age self.fluxGrid=np.zeros([len(self.ages), len(self.wavelengths)]) for line in lines: @@ -724,34 +714,33 @@ def __init__(self, fileName, ageColumn = 0, wavelengthColumn = 1, fluxColumn = 2 sedAge=float(bits[ageColumn]) sedWavelength=float(bits[wavelengthColumn]) sedFlux=float(bits[fluxColumn]) - + row=self.ages.index(sedAge) column=self.wavelengths.index(sedWavelength) - + self.fluxGrid[row][column]=sedFlux def getSED(self, ageGyr, z = 0.0, normalise = False, label = None): """Extract a SED for given age. Do linear interpolation between models if necessary. - - @type ageGyr: float - @param ageGyr: age of the SED in Gyr - @type z: float - @param z: redshift the SED from z = 0 to z = z - @type normalise: bool - @param normalise: normalise the SED to have area 1 - @rtype: L{SED} object - @return: SED - + + Args: + ageGyr (float): age of the SED in Gyr + z (float): redshift the SED from z = 0 to z = z + normalise (bool): normalise the SED to have area 1 + + Returns: + SED: SED object + """ - + if ageGyr in self.ages: - + flux=self.fluxGrid[self.ages.index(ageGyr)] sed=SED(self.wavelengths, flux, z = z, normalise = normalise, label = label) return sed - + else: - + # Use interpolation, iterating over each wavelength column flux=[] for i in range(len(self.wavelengths)): @@ -764,22 +753,19 @@ def getSED(self, ageGyr, z = 0.0, normalise = False, label = None): def getColourEvolution(self, passband1, passband2, zFormation, zStepSize = 0.05, magType = "Vega"): """Calculates the evolution of the colour observed through passband1 - passband2 for the StellarPopulation with redshift, from z = 0 to z = zFormation. - - @type passband1: L{Passband} object - @param passband1: filter passband through which to calculate the first magnitude - @type passband2: L{Passband} object - @param passband2: filter passband through which to calculate the second magnitude - @type zFormation: float - @param zFormation: formation redshift of the StellarPopulation - @type zStepSize: float - @param zStepSize: size of interval in z at which to calculate model colours - @type magType: string - @param magType: either "Vega" or "AB" - @rtype: dictionary - @return: dictionary of np.arrays in format {'z', 'colour'} - + + Args: + passband1 (Passband): filter passband through which to calculate the first magnitude + passband2 (Passband): filter passband through which to calculate the second magnitude + zFormation (float): formation redshift of the StellarPopulation + zStepSize (float): size of interval in z at which to calculate model colours + magType (str): either "Vega" or "AB" + + Returns: + dict: dictionary of np.arrays in format ``{'z', 'colour'}`` + """ - + zSteps=int(math.ceil(zFormation/zStepSize)) zData=[] colourData=[] @@ -793,34 +779,29 @@ def getColourEvolution(self, passband1, passband2, zFormation, zStepSize = 0.05, zData=np.array(zData) colourData=np.array(colourData) - + return {'z': zData, 'colour': colourData} - - def getMagEvolution(self, passband, magNormalisation, zNormalisation, zFormation, zStepSize = 0.05, + + def getMagEvolution(self, passband, magNormalisation, zNormalisation, zFormation, zStepSize = 0.05, onePlusZSteps = False, magType = "Vega"): """Calculates the evolution with redshift (from z = 0 to z = zFormation) of apparent magnitude - in the observed frame through the passband for the StellarPopulation, normalised to magNormalisation + in the observed frame through the passband for the StellarPopulation, normalised to magNormalisation (apparent) at z = zNormalisation. - - @type passband: L{Passband} object - @param passband: filter passband through which to calculate the magnitude - @type magNormalisation: float - @param magNormalisation: sets the apparent magnitude of the SED at zNormalisation - @type zNormalisation: float - @param zNormalisation: the redshift at which the magnitude normalisation is carried out - @type zFormation: float - @param zFormation: formation redshift of the StellarPopulation - @type zStepSize: float - @param zStepSize: size of interval in z at which to calculate model magnitudes - @type onePlusZSteps: bool - @param onePlusZSteps: if True, zSteps are (1+z)*zStepSize, otherwise zSteps are linear - @type magType: string - @param magType: either "Vega" or "AB" - @rtype: dictionary - @return: dictionary of np.arrays in format {'z', 'mag'} - + + Args: + passband (Passband): filter passband through which to calculate the magnitude + magNormalisation (float): sets the apparent magnitude of the SED at zNormalisation + zNormalisation (float): the redshift at which the magnitude normalisation is carried out + zFormation (float): formation redshift of the StellarPopulation + zStepSize (float): size of interval in z at which to calculate model magnitudes + onePlusZSteps (bool): if True, zSteps are (1+z)*zStepSize, otherwise zSteps are linear + magType (str): either "Vega" or "AB" + + Returns: + dict: dictionary of np.arrays in format ``{'z', 'mag'}`` + """ - + # Count upwards in z steps as interpolation doesn't work if array ordered z decreasing zSteps=int(math.ceil(zFormation/zStepSize)) zData=[] @@ -844,13 +825,13 @@ def getMagEvolution(self, passband, magNormalisation, zNormalisation, zFormation zData=np.array(zData) magData=np.array(magData) - + # Do the normalisation interpolator=interpolate.interp1d(zData, magData, kind='linear') modelNormMag=interpolator(zNormalisation) normConstant=magNormalisation-modelNormMag magData=magData+normConstant - + return {'z': zData, 'mag': magData} def calcEvolutionCorrection(self, zFrom, zTo, zFormation, passband, magType = "Vega"): @@ -858,58 +839,55 @@ def calcEvolutionCorrection(self, zFrom, zTo, zFormation, passband, magType = "V from redshift zFrom to redshift zTo, where the stellarPopulation is assumed to be formed at redshift zFormation. - @type zFrom: float - @param zFormation: redshift to evolution correct from - @type zTo: float - @param zTo: redshift to evolution correct to - @type zFormation: float - @param zFormation: formation redshift of the StellarPopulation - @type passband: L{Passband} object - @param passband: filter passband through which to calculate magnitude - @type magType: string - @param magType: either "Vega" or "AB" - @rtype: float - @return: evolution correction in magnitudes in the rest frame - + Args: + zFrom (float): redshift to evolution correct from + zTo (float): redshift to evolution correct to + zFormation (float): formation redshift of the StellarPopulation + passband (Passband): filter passband through which to calculate magnitude + magType (str): either "Vega" or "AB" + + Returns: + float: evolution correction in magnitudes in the rest frame + """ - + ageFrom=astCalc.tl(zFormation)-astCalc.tl(zFrom) ageTo=astCalc.tl(zFormation)-astCalc.tl(zTo) - + fromSED=self.getSED(ageFrom) toSED=self.getSED(ageTo) - + fromMag=fromSED.calcMag(passband, magType = magType, addDistanceModulus = False) toMag=toSED.calcMag(passband, magType = magType, addDistanceModulus = False) - + return fromMag-toMag - + #------------------------------------------------------------------------------------------------------------ class M05Model(StellarPopulation): """This class describes a Maraston 2005 stellar population model. To load a composite stellar population model (CSP) for a tau = 0.1 Gyr burst of star formation, solar metallicity, Salpeter IMF: - + m05csp = astSED.M05Model(M05_DIR+"/csp_e_0.10_z02_salp.sed_agb") - + where M05_DIR is set to point to the directory where the Maraston 2005 models are stored on your system. - + The file format of the Maraston 2005 simple stellar poulation (SSP) models is different to the file format used for the CSPs, and this needs to be specified using the fileType parameter. To load a SSP with solar metallicity, red horizontal branch morphology: - + m05ssp = astSED.M05Model(M05_DIR+"/sed.ssz002.rhb", fileType = "ssp") - + The wavelength units of SEDs from M05 models are Angstroms, with flux in units of erg/s/Angstrom. - + """ def __init__(self, fileName, fileType = "csp"): - + self.modelFamily="M05" inFile=open(fileName, "r") lines=inFile.readlines() inFile.close() - + self.fileName=fileName if fileType == "csp": @@ -922,7 +900,7 @@ def __init__(self, fileName, fileType = "csp"): fluxColumn=3 else: raise Exception("fileType must be 'ssp' or 'csp'") - + # Extract a list of model ages and valid wavelengths from the file self.ages=[] self.wavelengths=[] @@ -935,7 +913,7 @@ def __init__(self, fileName, fileType = "csp"): self.ages.append(age) if wavelength not in self.wavelengths: self.wavelengths.append(wavelength) - + # Construct a grid of flux - rows correspond to each wavelength, columns to age self.fluxGrid=np.zeros([len(self.ages), len(self.wavelengths)]) for line in lines: @@ -944,39 +922,39 @@ def __init__(self, fileName, fileType = "csp"): sedAge=float(bits[ageColumn]) sedWavelength=float(bits[wavelengthColumn]) sedFlux=float(bits[fluxColumn]) - + row=self.ages.index(sedAge) column=self.wavelengths.index(sedWavelength) - + self.fluxGrid[row][column]=sedFlux - + #------------------------------------------------------------------------------------------------------------ class BC03Model(StellarPopulation): """This class describes a Bruzual & Charlot 2003 stellar population model, extracted from a GALAXEV .ised file using the galaxevpl program that is included in GALAXEV. The file format is white space delimited, with wavelength in the first column. Subsequent columns contain the model fluxes for SEDs of different - ages, as specified when running galaxevpl. The age corresponding to each flux column is taken from the + ages, as specified when running galaxevpl. The age corresponding to each flux column is taken from the comment line beginning "# Age (yr)", and is converted to Gyr. For example, to load a tau = 0.1 Gyr burst of star formation, solar metallicity, Salpeter IMF model stored in a file (created by galaxevpl) called "csp_lr_solar_0p1Gyr.136": - + bc03model = BC03Model("csp_lr_solar_0p1Gyr.136") - The wavelength units of SEDs from BC03 models are Angstroms. Flux is converted into units of + The wavelength units of SEDs from BC03 models are Angstroms. Flux is converted into units of erg/s/Angstrom (the units in the files output by galaxevpl are LSun/Angstrom). """ - + def __init__(self, fileName): - + self.modelFamily="BC03" self.fileName=fileName inFile=open(fileName, "r") lines=inFile.readlines() inFile.close() - + # Extract a list of model ages - BC03 ages are in years, so convert to Gyr self.ages=[] for line in lines: @@ -984,7 +962,7 @@ def __init__(self, fileName): rawAges=line[line.find("# Age (yr)")+10:].split() for age in rawAges: self.ages.append(float(age)/1e9) - + # Extract a list of valid wavelengths from the file # If we have many ages in the file, this is more complicated... lambdaLinesCount=0 @@ -999,15 +977,15 @@ def __init__(self, fileName): for i in range(startFluxDataLine, len(lines), lambdaLinesCount): line=lines[i] bits=line.split() - self.wavelengths.append(float(bits[0])) - + self.wavelengths.append(float(bits[0])) + # Construct a grid of flux - rows correspond to each wavelength, columns to age self.fluxGrid=np.zeros([len(self.ages), len(self.wavelengths)]) for i in range(startFluxDataLine, len(lines), lambdaLinesCount): line=lines[i] bits=[] for k in range(i, i+lambdaLinesCount): - bits=bits+lines[k].split() + bits=bits+lines[k].split() ageFluxes=bits[1:] sedWavelength=float(bits[0]) column=self.wavelengths.index(sedWavelength) @@ -1017,24 +995,24 @@ def __init__(self, fileName): # Convert flux into erg/s/Angstrom - native units of galaxevpl files are LSun/Angstrom self.fluxGrid=self.fluxGrid*3.826e33 - + #------------------------------------------------------------------------------------------------------------ class P09Model(StellarPopulation): - """This class describes a Percival et al 2009 (BaSTI; http://albione.oa-teramo.inaf.it) stellar - population model. We assume that the synthetic spectra for each model are unpacked under the directory + """This class describes a Percival et al 2009 (BaSTI; http://albione.oa-teramo.inaf.it) stellar + population model. We assume that the synthetic spectra for each model are unpacked under the directory pointed to by fileName. - - The wavelength units of SEDs from P09 models are converted to Angstroms. Flux is converted into units of + + The wavelength units of SEDs from P09 models are converted to Angstroms. Flux is converted into units of erg/s/Angstrom (the units in the BaSTI low-res spectra are 4.3607e-33 erg/s/m). - + """ - + def __init__(self, fileName): - + self.modelFamily="P09" files=glob.glob(fileName+os.path.sep+"*.t??????") - + self.fileName=fileName # Map end of filenames to ages in Gyr @@ -1046,7 +1024,7 @@ def __init__(self, fileName): self.ages.append(ageGyr) extensionAgeMap[ext]=ageGyr self.ages.sort() - + # Construct a grid of flux - rows correspond to each wavelength, columns to age self.wavelengths=None self.fluxGrid=None @@ -1067,43 +1045,41 @@ def __init__(self, fileName): self.wavelengths=wavelength if self.fluxGrid == None: self.fluxGrid=np.zeros([len(self.ages), len(self.wavelengths)]) - self.fluxGrid[i]=flux + self.fluxGrid[i]=flux # Convert flux into erg/s/Angstrom - native units in BaSTI files are 4.3607e-33 erg/s/m self.fluxGrid=self.fluxGrid/4.3607e-33/1e10 - + #------------------------------------------------------------------------------------------------------------ def makeModelSEDDictList(modelList, z, passbandsList, labelsList = [], EBMinusVList = [0.0], forceYoungerThanUniverse = True): - """This routine makes a list of SEDDict dictionaries (see L{mags2SEDDict}) for fitting using - L{fitSEDDict}. This speeds up the fitting as this allows us to calculate model SED magnitudes only once, + """This routine makes a list of SEDDict dictionaries (see mags2SEDDict) for fitting using + fitSEDDict. This speeds up the fitting as this allows us to calculate model SED magnitudes only once, if all objects to be fitted are at the same redshift. We add some meta data to the modelSEDDicts (e.g. the model file names). - - The effect of extinction by dust (assuming the Calzetti et al. 2000 law) can be included by giving a list + + The effect of extinction by dust (assuming the Calzetti et al. 2000 law) can be included by giving a list of E(B-V) values. - + If forceYoungerThanUniverse == True, ages which are older than the universe at the given z will not be included. - - @type modelList: list of L{StellarPopulation} model objects - @param modelList: list of StellarPopulation models to include - @type z: float - @param z: redshift to apply to all stellar population models in modelList - @type EBMinusVList: list - @param EBMinusVList: list of E(B-V) extinction values to apply to all models, in magnitudes - @type labelsList: list - @param labelsList: optional list used for labelling passbands in output SEDDicts - @type forceYoungerThanUniverse: bool - @param forceYoungerThanUniverse: if True, do not allow models that exceed the age of the universe at z - @rtype: list - @return: list of dictionaries containing model fluxes, to be used as input to L{fitSEDDict}. - + + Args: + modelList (list): list of StellarPopulation models to include + z (float): redshift to apply to all stellar population models in modelList + passbandsList (list): list of Passband objects + labelsList (list): optional list used for labelling passbands in output SEDDicts + EBMinusVList (list): list of E(B-V) extinction values to apply to all models, in magnitudes + forceYoungerThanUniverse (bool): if True, do not allow models that exceed the age of the universe at z + + Returns: + list: list of dictionaries containing model fluxes, to be used as input to fitSEDDict + """ - + # Otherwise if this is the case we won't actually make any model SEDDicts ... if EBMinusVList == []: EBMinusVList=[0.0] - + modelSEDDictList=[] for m in range(len(modelList)): testAges=np.array(modelList[m].ages) @@ -1121,37 +1097,39 @@ def makeModelSEDDictList(modelList, z, passbandsList, labelsList = [], EBMinusVL modelSEDDict['E(B-V)']=EBMinusV modelSEDDict['ageGyr']=t modelSEDDict['z']=z - modelSEDDict['fileName']=modelList[m].fileName + modelSEDDict['fileName']=modelList[m].fileName modelSEDDict['modelListIndex']=m modelSEDDictList.append(modelSEDDict) - + return modelSEDDictList - + #------------------------------------------------------------------------------------------------------------ def fitSEDDict(SEDDict, modelSEDDictList): - """Fits the given SED dictionary (made using L{mags2SEDDict}) with the given list of model SED - dictionaries. The latter should be made using L{makeModelSEDDictList}, and entries for fluxes should + """Fits the given SED dictionary (made using mags2SEDDict) with the given list of model SED + dictionaries. The latter should be made using makeModelSEDDictList, and entries for fluxes should correspond directly between the model and SEDDict. - + Returns a dictionary with best fit values. - - @type SEDDict: dictionary, in format of L{mags2SEDDict} - @param SEDDict: dictionary of observed fluxes and uncertainties, in format of L{mags2SEDDict} - @type modelSEDDictList: list of dictionaries, in format of L{makeModelSEDDictList} - @param modelSEDDictList: list of dictionaries containing fluxes of models to be fitted to the observed fluxes listed in the SEDDict. This should be made using L{makeModelSEDDictList}. - @rtype: dictionary - @return: results of the fitting - keys: - - 'minChiSq': minimum chi squared value of best fit - - 'chiSqContrib': corresponding contribution at each passband to the minimum chi squared value - - 'ageGyr': the age in Gyr of the best fitting model - - 'modelFileName': the file name of the stellar population model corresponding to the best fit - - 'modelListIndex': the index of the best fitting model in the input modelSEDDictList - - 'norm': the normalisation that the best fit model should be multiplied by to match the SEDDict - - 'z': the redshift of the best fit model - - 'E(B-V)': the extinction, E(B-V), in magnitudes, of the best fit model - + + Args: + SEDDict (dict): dictionary of observed fluxes and uncertainties, in format of mags2SEDDict + modelSEDDictList (list): list of dictionaries containing fluxes of models to be fitted to the + observed fluxes listed in the SEDDict, made using makeModelSEDDictList + + Returns: + dict: results of the fitting - keys: + + - ``'minChiSq'``: minimum chi squared value of best fit + - ``'chiSqContrib'``: corresponding contribution at each passband to the minimum chi squared value + - ``'ageGyr'``: the age in Gyr of the best fitting model + - ``'modelFileName'``: the file name of the stellar population model corresponding to the best fit + - ``'modelListIndex'``: the index of the best fitting model in the input modelSEDDictList + - ``'norm'``: the normalisation that the best fit model should be multiplied by to match the SEDDict + - ``'z'``: the redshift of the best fit model + - ``'E(B-V)'``: the extinction, E(B-V), in magnitudes, of the best fit model + """ - + modelFlux=[] for modelSEDDict in modelSEDDictList: modelFlux.append(modelSEDDict['flux']) @@ -1170,37 +1148,38 @@ def fitSEDDict(SEDDict, modelSEDDictList): bestChiSq=minChiSq bestChiSqContrib=((sedFlux[bestMatchIndex]-norms[bestMatchIndex]*modelFlux[bestMatchIndex])**2)\ /sedFluxErr[bestMatchIndex]**2 - - resultsDict={'minChiSq': bestChiSq, + + resultsDict={'minChiSq': bestChiSq, 'chiSqContrib': bestChiSqContrib, 'allChiSqValues': chiSq, - 'ageGyr': modelSEDDictList[bestMatchIndex]['ageGyr'], + 'ageGyr': modelSEDDictList[bestMatchIndex]['ageGyr'], 'modelFileName': modelSEDDictList[bestMatchIndex]['fileName'], 'modelListIndex': modelSEDDictList[bestMatchIndex]['modelListIndex'], - 'norm': bestNorm, - 'z': modelSEDDictList[bestMatchIndex]['z'], + 'norm': bestNorm, + 'z': modelSEDDictList[bestMatchIndex]['z'], 'E(B-V)': modelSEDDictList[bestMatchIndex]['E(B-V)']} - + return resultsDict - + #------------------------------------------------------------------------------------------------------------ def mags2SEDDict(ABMags, ABMagErrs, passbands): """Takes a set of corresponding AB magnitudes, uncertainties, and passbands, and - returns a dictionary with keys 'flux', 'fluxErr', 'wavelength'. Fluxes are in units of + returns a dictionary with keys 'flux', 'fluxErr', 'wavelength'. Fluxes are in units of erg/s/cm^2/Angstrom, wavelength in Angstroms. These dictionaries are the staple diet of the - L{fitSEDDict} routine. - - @type ABMags: list or np array - @param ABMags: AB magnitudes, specified in corresponding order to passbands and ABMagErrs - @type ABMagErrs: list or np array - @param ABMagErrs: AB magnitude errors, specified in corresponding order to passbands and ABMags - @type passbands: list of L{Passband} objects - @param passbands: passband objects, specified in corresponding order to ABMags and ABMagErrs - @rtype: dictionary - @return: dictionary with keys {'flux', 'fluxErr', 'wavelength'}, suitable for input to L{fitSEDDict} + fitSEDDict routine. + + Args: + ABMags (list or numpy.ndarray): AB magnitudes, specified in corresponding order to passbands + and ABMagErrs + ABMagErrs (list or numpy.ndarray): AB magnitude errors, specified in corresponding order to + passbands and ABMags + passbands (list): list of Passband objects, specified in corresponding order to ABMags and ABMagErrs + + Returns: + dict: dictionary with keys ``{'flux', 'fluxErr', 'wavelength'}``, suitable for input to fitSEDDict """ - + flux=[] fluxErr=[] wavelength=[] @@ -1209,34 +1188,33 @@ def mags2SEDDict(ABMags, ABMagErrs, passbands): flux.append(f) fluxErr.append(err) wavelength.append(p.effectiveWavelength()) - + SEDDict={} SEDDict['flux']=np.array(flux) SEDDict['fluxErr']=np.array(fluxErr) SEDDict['wavelength']=np.array(wavelength) - + return SEDDict - + #------------------------------------------------------------------------------------------------------------ def mag2Flux(ABMag, ABMagErr, passband): """Converts given AB magnitude and uncertainty into flux, in erg/s/cm^2/Angstrom. - - @type ABMag: float - @param ABMag: magnitude on AB system in passband - @type ABMagErr: float - @param ABMagErr: uncertainty in AB magnitude in passband - @type passband: L{Passband} object - @param passband: L{Passband} object at which ABMag was measured - @rtype: list - @return: [flux, fluxError], in units of erg/s/cm^2/Angstrom - + + Args: + ABMag (float): magnitude on AB system in passband + ABMagErr (float): uncertainty in AB magnitude in passband + passband (Passband): Passband object at which ABMag was measured + + Returns: + list: [flux, fluxError], in units of erg/s/cm^2/Angstrom + """ - + fluxJy=(10**23.0)*10**(-(ABMag+48.6)/2.5) # AB mag aLambda=3e-13 # for conversion to erg s-1 cm-2 angstrom-1 with lambda in microns effLMicron=passband.effectiveWavelength()*(1e-10/1e-6) fluxWLUnits=aLambda*fluxJy/effLMicron**2 - + fluxJyErr=(10**23.0)*10**(-(ABMag-ABMagErr+48.6)/2.5) # AB mag fluxWLUnitsErr=aLambda*fluxJyErr/effLMicron**2 fluxWLUnitsErr=fluxWLUnitsErr-fluxWLUnits @@ -1246,16 +1224,15 @@ def mag2Flux(ABMag, ABMagErr, passband): #------------------------------------------------------------------------------------------------------------ def flux2Mag(flux, fluxErr, passband): """Converts given flux and uncertainty in erg/s/cm^2/Angstrom into AB magnitudes. - - @type flux: float - @param flux: flux in erg/s/cm^2/Angstrom in passband - @type fluxErr: float - @param fluxErr: uncertainty in flux in passband, in erg/s/cm^2/Angstrom - @type passband: L{Passband} object - @param passband: L{Passband} object at which ABMag was measured - @rtype: list - @return: [ABMag, ABMagError], in AB magnitudes - + + Args: + flux (float): flux in erg/s/cm^2/Angstrom in passband + fluxErr (float): uncertainty in flux in passband, in erg/s/cm^2/Angstrom + passband (Passband): Passband object at which ABMag was measured + + Returns: + list: [ABMag, ABMagError], in AB magnitudes + """ # aLambda = 3x10-5 for effective wavelength in angstroms @@ -1264,56 +1241,57 @@ def flux2Mag(flux, fluxErr, passband): fluxJy=(flux*effLMicron**2)/aLambda mag=-2.5*np.log10(fluxJy/10**23)-48.6 - + fluxErrJy=(fluxErr*effLMicron**2)/aLambda magErr=mag-(-2.5*np.log10((fluxJy+fluxErrJy)/10**23)-48.6) - + return [mag, magErr] #------------------------------------------------------------------------------------------------------------ def mag2Jy(ABMag): - """Converts an AB magnitude into flux density in Jy - - @type ABMag: float - @param ABMag: AB magnitude - @rtype: float - @return: flux density in Jy - + """Converts an AB magnitude into flux density in Jy. + + Args: + ABMag (float): AB magnitude + + Returns: + float: flux density in Jy + """ - + fluxJy=((10**23)*10**(-(float(ABMag)+48.6)/2.5)) - + return fluxJy #------------------------------------------------------------------------------------------------------------ def Jy2Mag(fluxJy): - """Converts flux density in Jy into AB magnitude - - @type fluxJy: float - @param fluxJy: flux density in Jy - @rtype: float - @return: AB magnitude - + """Converts flux density in Jy into AB magnitude. + + Args: + fluxJy (float): flux density in Jy + + Returns: + float: AB magnitude + """ - + ABMag=-2.5*(np.log10(fluxJy)-23.0)-48.6 - + return ABMag - + #------------------------------------------------------------------------------------------------------------ # Data VEGA=VegaSED() -"""The SED of Vega, used for calculation of magnitudes on the Vega system (L{SED} object).""" +"""The SED of Vega, used for calculation of magnitudes on the Vega system (SED object).""" # AB SED has constant flux density 3631 Jy AB=SED(wavelength = np.logspace(1, 8, int(1e5)), flux = np.ones(1000000)) -"""Flat spectrum SED, used for calculation of magnitudes on the AB system (L{SED} object).""" +"""Flat spectrum SED, used for calculation of magnitudes on the AB system (SED object).""" AB.flux=(3e-5*3631)/(AB.wavelength**2) AB.z0flux=AB.flux[:] # Solar SED from HST CALSPEC (http://www.stsci.edu/hst/observatory/cdbs/calspec.html) SOL=SED() -"""The SED of the Sun (L{SED} object).""" +"""The SED of the Sun (SED object).""" SOL.loadFromFile(astLib.__path__[0]+os.path.sep+"data"+os.path.sep+"sun_reference_stis_001.ascii") - diff --git a/astLib/astStats.py b/astLib/astStats.py index 02ccc73..f99a8fe 100755 --- a/astLib/astStats.py +++ b/astLib/astStats.py @@ -1,6 +1,6 @@ """Module for performing statistical calculations. -(c) 2007-2012 Matt Hilton +(c) 2007-2012 Matt Hilton (c) 2013-2014 Matt Hilton & Steven Boada @@ -11,10 +11,10 @@ Some routines may fail if they are passed lists with few items and encounter a `divide by zero' error. Where this occurs, the function will return None. An error message will be printed to the console when this happens if astStats.REPORT_ERRORS=True (the default). Testing if an -astStats function returns None can be used to handle errors in scripts. +astStats function returns None can be used to handle errors in scripts. -For extensive statistics modules, the Python bindings for GNU R (U{http://rpy.sourceforge.net}), or -SciPy (U{http://www.scipy.org}) are suggested. +For extensive statistics modules, the Python bindings for GNU R (http://rpy.sourceforge.net), or +SciPy (http://www.scipy.org) are suggested. """ @@ -27,25 +27,27 @@ #--------------------------------------------------------------------------------------------------- def mean(dataList): """Calculates the mean average of a list of numbers. - - @type dataList: list or numpy array - @param dataList: input data, must be a one dimensional list - @rtype: float - @return: mean average - + + Args: + dataList (list or numpy.ndarray): input data, must be a one dimensional list + + Returns: + float: mean average + """ return numpy.mean(dataList) - + #--------------------------------------------------------------------------------------------------- def weightedMean(dataList): """Calculates the weighted mean average of a two dimensional list (value, weight) of numbers. - - @type dataList: list - @param dataList: input data, must be a two dimensional list in format [value, weight] - @rtype: float - @return: weighted mean average - + + Args: + dataList (list): input data, must be a two dimensional list in format [value, weight] + + Returns: + float: weighted mean average + """ sum=0 weightSum=0 @@ -61,24 +63,26 @@ def weightedMean(dataList): #--------------------------------------------------------------------------------------------------- def stdev(dataList): """Calculates the (sample) standard deviation of a list of numbers. - - @type dataList: list or numpy array - @param dataList: input data, must be a one dimensional list - @rtype: float - @return: standard deviation - + + Args: + dataList (list or numpy.ndarray): input data, must be a one dimensional list + + Returns: + float: standard deviation + """ return numpy.std(dataList) - + #--------------------------------------------------------------------------------------------------- def rms(dataList): """Calculates the root mean square of a list of numbers. - - @type dataList: list - @param dataList: input data, must be a one dimensional list - @rtype: float - @return: root mean square - + + Args: + dataList (list): input data, must be a one dimensional list + + Returns: + float: root mean square + """ dataListSq=[] for item in dataList: @@ -87,18 +91,17 @@ def rms(dataList): rms=math.sqrt(listMeanSq) return rms - + #--------------------------------------------------------------------------------------------------- def weightedStdev(dataList): - """Calculates the weighted (sample) standard deviation of a list of numbers. - - @type dataList: list - @param dataList: input data, must be a two dimensional list in format [value, weight] - @rtype: float - @return: weighted standard deviation - - @note: Returns None if an error occurs. - + """Calculates the weighted (sample) standard deviation of a list of numbers. + + Args: + dataList (list): input data, must be a two dimensional list in format [value, weight] + + Returns: + float or None: weighted standard deviation, or None if an error occurs. + """ listMean=weightedMean(dataList) sum=0 @@ -108,7 +111,7 @@ def weightedStdev(dataList): if item[1]>0.0: sum=sum+float((item[0]-listMean)/item[1])*float((item[0]-listMean)/item[1]) wSum=wSum+float(1.0/item[1])*float(1.0/item[1]) - + if len(dataList)>1: nFactor=float(len(dataList))/float(len(dataList)-1) stdev=math.sqrt(nFactor*(sum/wSum)) @@ -117,28 +120,30 @@ def weightedStdev(dataList): print("""ERROR: astStats.weightedStdev() : dataList contains < 2 items.""") stdev=None return stdev - + #--------------------------------------------------------------------------------------------------- def median(dataList): """Calculates the median of a list of numbers. - - @type dataList: list or numpy array - @param dataList: input data, must be a one dimensional list - @rtype: float - @return: median average - - """ + + Args: + dataList (list or numpy.ndarray): input data, must be a one dimensional list + + Returns: + float: median average + + """ return numpy.median(dataList) - + #--------------------------------------------------------------------------------------------------- def modeEstimate(dataList): """Returns an estimate of the mode of a set of values by mode=(3*median)-(2*mean). - - @type dataList: list - @param dataList: input data, must be a one dimensional list - @rtype: float - @return: estimate of mode average - + + Args: + dataList (list): input data, must be a one dimensional list + + Returns: + float: estimate of mode average + """ mode=(3*median(dataList))-(2*mean(dataList)) @@ -147,39 +152,38 @@ def modeEstimate(dataList): #--------------------------------------------------------------------------------------------------- def MAD(dataList): """Calculates the Median Absolute Deviation of a list of numbers. - - @type dataList: list - @param dataList: input data, must be a one dimensional list - @rtype: float - @return: median absolute deviation - + + Args: + dataList (list): input data, must be a one dimensional list + + Returns: + float: median absolute deviation + """ listMedian=median(dataList) - + # Calculate |x-M| values diffModuli=[] for item in dataList: diffModuli.append(math.fabs(item-listMedian)) - + MAD=median(diffModuli) - + return MAD #--------------------------------------------------------------------------------------------------- def biweightLocation(dataList, tuningConstant): """Calculates the biweight location estimator (like a robust average) of a list of numbers. - - @type dataList: list - @param dataList: input data, must be a one dimensional list - @type tuningConstant: float - @param tuningConstant: 6.0 is recommended. - @rtype: float - @return: biweight location - - @note: Returns None if an error occurs. - - """ + + Args: + dataList (list): input data, must be a one dimensional list + tuningConstant (float): 6.0 is recommended. + + Returns: + float or None: biweight location, or None if an error occurs. + + """ C=tuningConstant listMedian=median(dataList) listMAD=MAD(dataList) @@ -187,7 +191,7 @@ def biweightLocation(dataList, tuningConstant): uValues=[] for item in dataList: uValues.append((item-listMedian)/(C*listMAD)) - + top=0 # numerator equation (5) Beers et al if you like bottom=0 # denominator for i in range(len(uValues)): @@ -195,36 +199,34 @@ def biweightLocation(dataList, tuningConstant): top=top+((dataList[i]-listMedian) \ *(1.0-(uValues[i]*uValues[i])) \ *(1.0-(uValues[i]*uValues[i]))) - + bottom=bottom+((1.0-(uValues[i]*uValues[i])) \ *(1.0-(uValues[i]*uValues[i]))) - + CBI=listMedian+(top/bottom) - + else: if REPORT_ERRORS==True: print("""ERROR: astStats: biweightLocation() : MAD() returned 0.""") return None - + return CBI #--------------------------------------------------------------------------------------------------- def biweightScale(dataList, tuningConstant): """Calculates the biweight scale estimator (like a robust standard deviation) of a list - of numbers. - - @type dataList: list - @param dataList: input data, must be a one dimensional list - @type tuningConstant: float - @param tuningConstant: 9.0 is recommended. - @rtype: float - @return: biweight scale - - @note: Returns None if an error occurs. - - """ + of numbers. + + Args: + dataList (list): input data, must be a one dimensional list + tuningConstant (float): 9.0 is recommended. + + Returns: + float or None: biweight scale, or None if an error occurs. + + """ C=tuningConstant - + # Calculate |x-M| values and u values listMedian=median(dataList) listMAD=MAD(dataList) @@ -239,11 +241,11 @@ def biweightScale(dataList, tuningConstant): if REPORT_ERRORS==True: print("""ERROR: astStats.biweightScale() : divide by zero error.""") return None - + top=0 # numerator equation (9) Beers et al bottom=0 valCount=0 # Count values where u<1 only - + for i in range(len(uValues)): # Skip u values >1 if math.fabs(uValues[i])<=1.0: @@ -252,7 +254,7 @@ def biweightScale(dataList, tuningConstant): top=top+((diffModuli[i]*diffModuli[i])*u4Term) bottom=bottom+(u2Term*(1.0-(5.0*(uValues[i]*uValues[i])))) valCount=valCount+1 - + top=math.sqrt(top) bottom=math.fabs(bottom) @@ -262,24 +264,22 @@ def biweightScale(dataList, tuningConstant): #--------------------------------------------------------------------------------------------------- def biweightClipped(dataList, tuningConstant, sigmaCut): """Iteratively calculates biweight location and scale, using sigma clipping, for a list - of values. The calculation is performed on the first column of a multi-dimensional + of values. The calculation is performed on the first column of a multi-dimensional list; other columns are ignored. - - @type dataList: list - @param dataList: input data, must contain more than 5 values - @type tuningConstant: float - @param tuningConstant: 6.0 is recommended for location estimates, 9.0 is recommended for - scale estimates - @type sigmaCut: float - @param sigmaCut: sigma clipping to apply - @rtype: dictionary - @return: estimate of biweight location, scale, and list of non-clipped data, in the format - {'biweightLocation', 'biweightScale', 'dataList'} - - @note: Returns None if an error occurs. - - """ - + + Args: + dataList (list): input data, must contain more than 5 values + tuningConstant (float): 6.0 is recommended for location estimates, 9.0 is recommended + for scale estimates + sigmaCut (float): sigma clipping to apply + + Returns: + dict or None: estimate of biweight location, scale, and list of non-clipped data in + the format ``{'biweightLocation', 'biweightScale', 'dataList'}``, or None if an + error occurs. + + """ + iterations=0 clippedValues=[] for row in dataList: @@ -287,28 +287,28 @@ def biweightClipped(dataList, tuningConstant, sigmaCut): clippedValues.append(row[0]) else: clippedValues.append(row) - + cbi=None sbi=None clippedData=None while iterations<11 and len(clippedValues)>5: - - cbi=biweightLocation(clippedValues, tuningConstant) + + cbi=biweightLocation(clippedValues, tuningConstant) sbi=biweightScale(clippedValues, tuningConstant) - + # check for either biweight routine falling over # happens when feed in lots of similar numbers - # e.g. when bootstrapping with a small sample + # e.g. when bootstrapping with a small sample if cbi==None or sbi==None: - + if REPORT_ERRORS==True: print("""ERROR: astStats : biweightClipped() : divide by zero error.""") - + return None - + else: - + clippedValues=[] clippedData=[] for row in dataList: @@ -322,27 +322,27 @@ def biweightClipped(dataList, tuningConstant, sigmaCut): and row 4: - + m=listCopy.mean() s=listCopy.std() - + listCopy=listCopy[numpy.less(abs(listCopy), abs(m+sigmaCut*s))] - + iterations=iterations+1 - + return {'clippedMean': m, 'clippedStdev': s, 'numPoints': listCopy.shape[0]} - + #--------------------------------------------------------------------------------------------------- def clippedWeightedLSFit(dataList, sigmaCut): """Performs a weighted least squares fit on a list of numbers with sigma clipping. Minimum number of data points is 5. - - @type dataList: list - @param dataList: input data, must be a three dimensional list in format [x, y, y weight] - @rtype: dictionary - @return: slope and intercept on y-axis, with associated errors, in the format - {'slope', 'intercept', 'slopeError', 'interceptError'} - - @note: Returns None if an error occurs. - + + Args: + dataList (list): input data, must be a three dimensional list in format [x, y, y weight] + sigmaCut (float): sigma clipping to apply + + Returns: + dict or None: slope and intercept on y-axis with associated errors in the format + ``{'slope', 'intercept', 'slopeError', 'interceptError'}``, or None if an + error occurs. + """ - + iterations=0 clippedValues=[] for row in dataList: clippedValues.append(row) - + while iterations<11 and len(clippedValues)>4: - + fitResults=weightedLSFit(clippedValues, "errors") - + if fitResults['slope'] == None: - + if REPORT_ERRORS==True: print("""ERROR: astStats : clippedWeightedLSFit() : divide by zero error.""") - + return None - + else: - + clippedValues=[] for row in dataList: - + # Trim points more than sigmaCut*sigma away from the fitted line fit=fitResults['slope']*row[0]+fitResults['intercept'] res=row[1]-fit if abs(res)/row[2] < sigmaCut: clippedValues.append(row) - + iterations=iterations+1 - + # store the number of values that made it through the clipping process fitResults['numDataPoints']=len(clippedValues) - + return fitResults - + #--------------------------------------------------------------------------------------------------- def weightedLSFit(dataList, weightType): """Performs a weighted least squares fit on a three dimensional list of numbers [x, y, y error]. - - @type dataList: list - @param dataList: input data, must be a three dimensional list in format [x, y, y error] - @type weightType: string - @param weightType: if "errors", weights are calculated assuming the input data is in the - format [x, y, error on y]; if "weights", the weights are assumed to be already calculated and - stored in a fourth column [x, y, error on y, weight] (as used by e.g. L{astStats.biweightLSFit}) - @rtype: dictionary - @return: slope and intercept on y-axis, with associated errors, in the format - {'slope', 'intercept', 'slopeError', 'interceptError'} - - @note: Returns None if an error occurs. - + + Args: + dataList (list): input data, must be a three dimensional list in format [x, y, y error] + weightType (str): if ``"errors"``, weights are calculated assuming the input data is in + the format [x, y, error on y]; if ``"weights"``, the weights are assumed to be already + calculated and stored in a fourth column [x, y, error on y, weight] (as used by e.g. + biweightLSFit) + + Returns: + dict or None: slope and intercept on y-axis with associated errors in the format + ``{'slope', 'intercept', 'slopeError', 'interceptError'}``, or None if an + error occurs. + """ if weightType == "weights": sumW=0 @@ -527,7 +528,7 @@ def weightedLSFit(dataList, weightType): sumWXX=sumWXX+(W*item[0]*item[0]) sumW=sumW+W #print sumW, sumWXX, sumWX - + try: m=((sumW*sumWXY)-(sumWX*sumWY)) \ /((sumW*sumWXX)-(sumWX*sumWX)) @@ -535,7 +536,7 @@ def weightedLSFit(dataList, weightType): if REPORT_ERRORS == True: print("ERROR: astStats.weightedLSFit() : divide by zero error.") return None - + try: c=((sumWXX*sumWY)-(sumWX*sumWXY)) \ /((sumW*sumWXX)-(sumWX*sumWX)) @@ -543,38 +544,38 @@ def weightedLSFit(dataList, weightType): if REPORT_ERRORS == True: print("ERROR: astStats.weightedLSFit() : divide by zero error.") return None - + sumRes=0 for item in dataList: - + sumRes=sumRes+((item[1]-(m*item[0])-c) \ *(item[1]-(m*item[0])-c)) - + sigma=math.sqrt((1.0/(n-2))*sumRes) - + # Can get div0 errors here so check # When biweight fitting converges this shouldn't happen - if (n*sumWXX)-(sumWX*sumWX)>0.0: - + if (n*sumWXX)-(sumWX*sumWX)>0.0: + mSigma=(sigma*math.sqrt(n)) \ /math.sqrt((n*sumWXX)-(sumWX*sumWX)) - + cSigma=(sigma*math.sqrt(sumWXX)) \ /math.sqrt((n*sumWXX)-(sumWX*sumWX)) - + else: - + if REPORT_ERRORS==True: print("""ERROR: astStats.weightedLSFit() : divide by zero error.""") return None - + else: if REPORT_ERRORS==True: print("""ERROR: astStats.weightedLSFit() : dataList contains < 5 items.""") return None - + elif weightType == "errors": sumX=0 sumY=0 @@ -588,46 +589,44 @@ def weightedLSFit(dataList, weightType): sumXY=sumXY+((item[0]*item[1])/(item[2]*item[2])) sumXX=sumXX+((item[0]*item[0])/(item[2]*item[2])) sumSigma=sumSigma+(1.0/(item[2]*item[2])) - delta=(sumSigma*sumXX)-(sumX*sumX) + delta=(sumSigma*sumXX)-(sumX*sumX) m=((sumSigma*sumXY)-(sumX*sumY))/delta c=((sumXX*sumY)-(sumX*sumXY))/delta mSigma=math.sqrt(sumSigma/delta) cSigma=math.sqrt(sumXX/delta) - + return {'slope':m, 'intercept':c, 'slopeError':mSigma, 'interceptError':cSigma} - + #--------------------------------------------------------------------------------------------------- def biweightLSFit(dataList, tuningConstant, sigmaCut = None): """Performs a weighted least squares fit, where the weights used are the biweight transforms of the residuals to the previous best fit .i.e. the procedure is iterative, and converges very quickly (iterations is set to 10 by default). Minimum number of data points is 10. - + This seems to give slightly different results to the equivalent R routine, so use at your own risk! - - @type dataList: list - @param dataList: input data, must be a three dimensional list in format [x, y, y weight] - @type tuningConstant: float - @param tuningConstant: 6.0 is recommended for location estimates, 9.0 is recommended for - scale estimates - @type sigmaCut: float - @param sigmaCut: sigma clipping to apply (set to None if not required) - @rtype: dictionary - @return: slope and intercept on y-axis, with associated errors, in the format - {'slope', 'intercept', 'slopeError', 'interceptError'} - - @note: Returns None if an error occurs. - + + Args: + dataList (list): input data, must be a three dimensional list in format [x, y, y weight] + tuningConstant (float): 6.0 is recommended for location estimates, 9.0 is recommended + for scale estimates + sigmaCut (float or None): sigma clipping to apply; set to None if not required + + Returns: + dict or None: slope and intercept on y-axis with associated errors in the format + ``{'slope', 'intercept', 'slopeError', 'interceptError'}``, or None if an + error occurs. + """ dataCopy=[] for row in dataList: dataCopy.append(row) - + # First perform unweighted fit, then calculate residuals results=OLSFit(dataCopy) origLen=len(dataCopy) @@ -637,9 +636,9 @@ def biweightLSFit(dataList, tuningConstant, sigmaCut = None): res=[] for item in dataCopy: res.append((m*item[0]+c)-item[1]) - + if len(res)>5: - # For clipping, trim away things >3 sigma + # For clipping, trim away things >3 sigma # away from median if sigmaCut != None: absRes=[] @@ -652,41 +651,40 @@ def biweightLSFit(dataList, tuningConstant, sigmaCut = None): and len(dataCopy)>2: del dataCopy[count] del res[count] - + # Index of datalist gets out of # sync with absRes as we delete # items - count=count-1 - + count=count-1 + count=count+1 - + # Biweight transform residuals weights=biweightTransform(res, tuningConstant) - - # Perform weighted fit, using biweight transforms + + # Perform weighted fit, using biweight transforms # of residuals as weight wData=[] for i in range(len(dataCopy)): wData.append([dataCopy[i][0], dataCopy[i][1], dataCopy[i][2], weights[i][1]]) - + results=weightedLSFit(wData, "weights") return results - + #--------------------------------------------------------------------------------------------------- def cumulativeBinner(data, binMin, binMax, binTotal): """Bins the input data cumulatively. - - @param data: input data, must be a one dimensional list - @type binMin: float - @param binMin: minimum value from which to bin data - @type binMax: float - @param binMax: maximum value from which to bin data - @type binTotal: int - @param binTotal: number of bins - @rtype: list - @return: binned data, in format [bin centre, frequency] - + + Args: + data (list): input data, must be a one dimensional list + binMin (float): minimum value from which to bin data + binMax (float): maximum value from which to bin data + binTotal (int): number of bins + + Returns: + list: binned data in format [bin centre, frequency] + """ #Bin data binStep=float(binMax-binMin)/binTotal @@ -697,28 +695,27 @@ def cumulativeBinner(data, binMin, binMax, binTotal): for item in data: if item>(binMin+(i*binStep)): bins[i]=bins[i]+1.0/totalItems - + # Gnuplot requires points at bin midpoints coords=[] for i in range(binTotal): coords.append([binMin+(float(i+0.5)*binStep), bins[i]]) - + return coords #--------------------------------------------------------------------------------------------------- def binner(data, binMin, binMax, binTotal): - """Bins the input data.. - - @param data: input data, must be a one dimensional list - @type binMin: float - @param binMin: minimum value from which to bin data - @type binMax: float - @param binMax: maximum value from which to bin data - @type binTotal: int - @param binTotal: number of bins - @rtype: list - @return: binned data, in format [bin centre, frequency] - + """Bins the input data. + + Args: + data (list): input data, must be a one dimensional list + binMin (float): minimum value from which to bin data + binMax (float): maximum value from which to bin data + binTotal (int): number of bins + + Returns: + list: binned data in format [bin centre, frequency] + """ #Bin data binStep=float(binMax-binMin)/binTotal @@ -729,28 +726,28 @@ def binner(data, binMin, binMax, binTotal): if item>(binMin+(i*binStep)) \ and item<=(binMin+((i+1)*binStep)): bins[i]=bins[i]+1 - + # Gnuplot requires points at bin midpoints coords=[] for i in range(binTotal): coords.append([binMin+(float(i+0.5)*binStep), bins[i]]) - + return coords #--------------------------------------------------------------------------------------------------- def weightedBinner(data, weights, binMin, binMax, binTotal): """Bins the input data, recorded frequency is sum of weights in bin. - - @param data: input data, must be a one dimensional list - @type binMin: float - @param binMin: minimum value from which to bin data - @type binMax: float - @param binMax: maximum value from which to bin data - @type binTotal: int - @param binTotal: number of bins - @rtype: list - @return: binned data, in format [bin centre, frequency] - + + Args: + data (list): input data, must be a one dimensional list + weights (list): weights corresponding to each data value + binMin (float): minimum value from which to bin data + binMax (float): maximum value from which to bin data + binTotal (int): number of bins + + Returns: + list: binned data in format [bin centre, frequency] + """ #Bin data binStep=float(binMax-binMin)/binTotal @@ -761,12 +758,12 @@ def weightedBinner(data, weights, binMin, binMax, binTotal): if item>(binMin+(i*binStep)) \ and item<=(binMin+((i+1)*binStep)): bins[i]=bins[i]+weight - + # Gnuplot requires points at bin midpoints coords=[] for i in range(binTotal): coords.append([binMin+(float(i+0.5)*binStep), bins[i]]) - + return coords - + #--------------------------------------------------------------------------------------------------- diff --git a/astLib/astWCS.py b/astLib/astWCS.py index 1384f8f..b944ef4 100644 --- a/astLib/astWCS.py +++ b/astLib/astWCS.py @@ -7,7 +7,7 @@ This is a higher level interface to some of the routines in PyWCSTools (distributed with astLib). PyWCSTools is a simple SWIG wrapping of WCSTools by Jessica Mink -(U{http://tdc-www.harvard.edu/software/wcstools/}). It is intended is to make +(http://tdc-www.harvard.edu/software/wcstools/). It is intended is to make this interface complete enough such that direct use of PyWCSTools is unnecessary. @@ -25,7 +25,7 @@ # FITS convention. NUMPY_MODE = True """If True (default), pixel coordinates accepted/returned by routines such as -L{astWCS.WCS.pix2wcs}, L{astWCS.WCS.wcs2pix} have (0, 0) as the origin. Set +astWCS.WCS.pix2wcs, astWCS.WCS.wcs2pix have (0, 0) as the origin. Set to False to make these routines accept/return pixel coords with (1, 1) as the origin (i.e. to match the FITS convention, default behaviour prior to astLib version 0.3.0).""" @@ -45,15 +45,15 @@ class WCS: Coordinate System (WCS) contained in the header of a FITS image. Conversions between pixel and WCS coordinates can also be performed. - To create a WCS object from a FITS file called "test.fits", simply: + To create a WCS object from a FITS file called "test.fits", simply:: - WCS=astWCS.WCS("test.fits") + WCS=astWCS.WCS("test.fits") - Likewise, to create a WCS object from the pyfits.header of "test.fits": + Likewise, to create a WCS object from the pyfits.header of "test.fits":: - img=pyfits.open("test.fits") - header=img[0].header - WCS=astWCS.WCS(header, mode = "pyfits") + img=pyfits.open("test.fits") + header=img[0].header + WCS=astWCS.WCS(header, mode = "pyfits") """ @@ -63,29 +63,28 @@ def __init__(self, headerSource, extensionName = 0, mode = "image", zapKeywords header of the specified .fits image, or from a pyfits.header object. Set mode = "pyfits" if the headerSource is a pyfits.header. - For some images from some archives, particular header keywords such as + For some images from some archives, particular header keywords such as COMMENT or HISTORY may contain unprintable strings. If you encounter this, try setting zapKeywords = ['COMMENT', 'HISTORY'] (for example). - - @type headerSource: string or pyfits.header - @param headerSource: filename of input .fits image, or a pyfits.header - object - @type extensionName: int or string - @param extensionName: name or number of .fits extension in which image - data is stored - @type mode: string - @param mode: set to "image" if headerSource is a .fits file name, or - set to "astropy" if headerSource is an astropy.io.fits.Header object - (setting to this to "pyfits" is equivalent, for backwards - compatibility with existing code) - @type zapKeywords: list - @param: zapKeywords: keywords to remove from the header before making - astWCS object. - @type: useAstropyWCS: bool - @param: useAstropyWCS: if True, use astropy.wcs to perform WCS - coordinate conversions in wcs2pix, pix2wcs (if False, use PyWCSTools) - - @note: The meta data provided by headerSource is stored in WCS.header + + Args: + headerSource (str or astropy.io.fits.Header): filename of input + .fits image, or an astropy.io.fits.Header object + extensionName (int or str): name or number of .fits extension in + which image data is stored + mode (str): set to ``"image"`` if headerSource is a .fits file + name, or set to ``"astropy"`` if headerSource is an + astropy.io.fits.Header object (setting to ``"pyfits"`` is + equivalent, for backwards compatibility with existing code) + zapKeywords (list): keywords to remove from the header before + making the WCS object + useAstropyWCS (bool): if True, use astropy.wcs to perform WCS + coordinate conversions in wcs2pix, pix2wcs (if False, use + PyWCSTools) + naxis (int): number of WCS axes to use + + Note: + The meta data provided by headerSource is stored in WCS.header as a pyfits.header object. """ @@ -112,14 +111,14 @@ def __init__(self, headerSource, extensionName = 0, mode = "image", zapKeywords for count in range(self.headerSource.count(z)): self.headerSource.remove(z) self.header=headerSource - + # Scan for CUNIT values that upset astropy.wcs for i in (1, 2): if 'CUNIT%d' % (i) in self.header.keys(): if self.header['CUNIT%d' % (i)] == '' or self.header['CUNIT%d' % (i)] == 'degree'\ or self.header['CUNIT%d' % (i)] == 'degrees' or self.header['CUNIT%d' % (i)] == 'DEG': self.header['CUNIT%d' % (i)]='deg' - + # This enables a shim to allow code written for astLib to use astropy.wcs underneath self.useAstropyWCS=useAstropyWCS if NUMPY_MODE == True: @@ -133,13 +132,13 @@ def __init__(self, headerSource, extensionName = 0, mode = "image", zapKeywords def copy(self): """Copies the WCS object to a new object. - @rtype: astWCS.WCS object - @return: WCS object + Returns: + WCS: a new WCS object copied from this one """ # This only sets up a new WCS object, doesn't do a deep copy - ret = WCS(self.headerSource, self.extensionName, self.mode, + ret = WCS(self.headerSource, self.extensionName, self.mode, useAstropyWCS = self.useAstropyWCS) # This fixes copy bug @@ -164,15 +163,15 @@ def updateFromHeader(self): newHead.append((i[0], i[1])) else: newHead.append(('HIERARCH '+i[0], i[1])) - + # Workaround for ZPN bug when PV2_3 == 0 (as in, e.g., ESO WFI images) if "PV2_3" in list(newHead.keys()) and newHead['PV2_3'] == 0 and newHead['CTYPE1'] == 'RA---ZPN': newHead["PV2_3"]=1e-15 - + cardstring = "" for card in newHead.cards: cardstring = cardstring+str(card) - + if self.useAstropyWCS == True: self.AWCS = apywcs.WCS(self.header, naxis = self.naxis) # For astropy.wcs shim self.WCSStructure = wcs.wcsinit(cardstring) @@ -182,8 +181,8 @@ def getCentreWCSCoords(self): """Returns the RA and dec coordinates (in decimal degrees) at the centre of the WCS. - @rtype: list - @return: coordinates in decimal degrees in format [RADeg, decDeg] + Returns: + list: coordinates in decimal degrees in format [RADeg, decDeg] """ full = wcs.wcsfull(self.WCSStructure) @@ -198,9 +197,9 @@ def getFullSizeSkyDeg(self): decimal degrees on the sky (i.e., with the projection taken into account). - @rtype: list - @return: width and height of image in decimal degrees on the sky in - format [width, height] + Returns: + list: width and height of image in decimal degrees on the sky in + format [width, height] """ full = wcs.wcsfull(self.WCSStructure) @@ -214,9 +213,9 @@ def getHalfSizeDeg(self): """Returns the half-width, half-height of the image according to the WCS in RA and dec degrees. - @rtype: list - @return: half-width and half-height of image in R.A., dec. decimal - degrees in format [half-width, half-height] + Returns: + list: half-width and half-height of image in R.A., dec. decimal + degrees in format [half-width, half-height] """ half = wcs.wcssize(self.WCSStructure) @@ -230,8 +229,8 @@ def getImageMinMaxWCSCoords(self): """Returns the minimum, maximum WCS coords defined by the size of the parent image (as defined by the NAXIS keywords in the image header). - @rtype: list - @return: [minimum R.A., maximum R.A., minimum Dec., maximum Dec.] + Returns: + list: [minimum R.A., maximum R.A., minimum Dec., maximum Dec.] """ @@ -263,8 +262,12 @@ def wcs2pix(self, RADeg, decDeg): coordinates (given in decimal degrees). RADeg, decDeg can be single floats, or lists or np arrays. - @rtype: list - @return: pixel coordinates in format [x, y] + Args: + RADeg (float or list or numpy.ndarray): R.A. in decimal degrees + decDeg (float or list or numpy.ndarray): dec. in decimal degrees + + Returns: + list: pixel coordinates in format [x, y] """ if self.useAstropyWCS == False: @@ -296,7 +299,7 @@ def wcs2pix(self, RADeg, decDeg): pixCoords[0] = pixCoords[0]-1 pixCoords[1] = pixCoords[1]-1 pixCoords = [pixCoords[0], pixCoords[1]] - + else: # astropy.wcs shim if self.naxis == 2: @@ -314,8 +317,12 @@ def pix2wcs(self, x, y): """Returns the WCS coordinates corresponding to the input pixel coordinates. - @rtype: list - @return: WCS coordinates in format [RADeg, decDeg] + Args: + x (float or list or numpy.ndarray): pixel x coordinate(s) + y (float or list or numpy.ndarray): pixel y coordinate(s) + + Returns: + list: WCS coordinates in format [RADeg, decDeg] """ if self.useAstropyWCS == False: @@ -342,7 +349,7 @@ def pix2wcs(self, x, y): else: raise Exception("Not handling NAXIS > 3 with astropy.wcs shim") WCSCoords = np.array(WCSCoords)[:2].transpose().tolist() - + return WCSCoords @@ -350,11 +357,15 @@ def coordsAreInImage(self, RADeg, decDeg): """Returns True if the given RA, dec coordinate is within the image boundaries. - @rtype: bool - @return: True if coordinate within image, False if not. + Args: + RADeg (float): R.A. in decimal degrees + decDeg (float): dec. in decimal degrees + + Returns: + bool: True if coordinate is within the image, False otherwise """ - + pixCoords = self.wcs2pix(RADeg, decDeg) if pixCoords[0] >= 0 and pixCoords[0] < self.header['NAXIS1'] and \ pixCoords[1] >= 0 and pixCoords[1] < self.header['NAXIS2']: @@ -367,8 +378,8 @@ def getRotationDeg(self): """Returns the rotation angle in degrees around the axis, North through East. - @rtype: float - @return: rotation angle in degrees + Returns: + float: rotation angle in degrees """ return self.WCSStructure.rot @@ -377,8 +388,8 @@ def getRotationDeg(self): def isFlipped(self): """Returns 1 if image is reflected around axis, otherwise returns 0. - @rtype: int - @return: 1 if image is flipped, 0 otherwise + Returns: + int: 1 if image is flipped, 0 otherwise """ return self.WCSStructure.imflip @@ -388,8 +399,8 @@ def getPixelSizeDeg(self): """Returns the pixel scale of the WCS. This is the average of the x, y pixel scales. - @rtype: float - @return: pixel size in decimal degrees + Returns: + float: pixel size in decimal degrees """ @@ -401,8 +412,8 @@ def getPixelSizeDeg(self): def getXPixelSizeDeg(self): """Returns the pixel scale along the x-axis of the WCS in degrees. - @rtype: float - @return: pixel size in decimal degrees + Returns: + float: pixel size in decimal degrees """ @@ -414,8 +425,8 @@ def getXPixelSizeDeg(self): def getYPixelSizeDeg(self): """Returns the pixel scale along the y-axis of the WCS in degrees. - @rtype: float - @return: pixel size in decimal degrees + Returns: + float: pixel size in decimal degrees """ @@ -427,8 +438,8 @@ def getYPixelSizeDeg(self): def getEquinox(self): """Returns the equinox of the WCS. - @rtype: float - @return: equinox of the WCS + Returns: + float: equinox of the WCS """ return self.WCSStructure.equinox @@ -437,8 +448,8 @@ def getEquinox(self): def getEpoch(self): """Returns the epoch of the WCS. - @rtype: float - @return: epoch of the WCS + Returns: + float: epoch of the WCS """ return self.WCSStructure.epoch @@ -451,10 +462,19 @@ def findWCSOverlap(wcs1, wcs2): wcs2. Returns these coordinates, plus the corresponding pixel coordinates for each wcs. Useful for clipping overlapping region between two images. - @rtype: dictionary - @return: dictionary with keys 'overlapWCS' (min, max RA, dec of overlap - between wcs1, wcs2) 'wcs1Pix', 'wcs2Pix' (pixel coords in each input - WCS that correspond to 'overlapWCS' coords) + Args: + wcs1 (WCS): first WCS object + wcs2 (WCS): second WCS object + + Returns: + dict: dictionary with keys: + + - ``'overlapWCS'``: [minRA, maxRA, minDec, maxDec] of the overlap + between wcs1 and wcs2 + - ``'wcs1Pix'``: pixel coords in wcs1 corresponding to the overlap + WCS coords + - ``'wcs2Pix'``: pixel coords in wcs2 corresponding to the overlap + WCS coords """ From bd394b81b10be5fe0bc7500cd51d25727e47bcae Mon Sep 17 00:00:00 2001 From: Matt Hilton Date: Tue, 26 May 2026 22:04:10 +0200 Subject: [PATCH 2/5] Fix Sphinx autodoc failure when C extension is not compiled Add autodoc_mock_imports for PyWCSTools._wcs and PyWCSTools.wcs so that Sphinx can import astWCS, astImages, and astPlots without requiring the compiled C extension. Also fix an over-indented RST bullet list in the addPlotObjects docstring that caused a docutils unexpected-indentation error. Co-Authored-By: Claude Sonnet 4.6 --- astLib/astPlots.py | 9 +++++---- docs/conf.py | 3 ++- 2 files changed, 7 insertions(+), 5 deletions(-) diff --git a/astLib/astPlots.py b/astLib/astPlots.py index b789921..f4937af 100755 --- a/astLib/astPlots.py +++ b/astLib/astPlots.py @@ -427,10 +427,11 @@ def addPlotObjects(self, objRAs, objDecs, tag, symbol="circle", size=4.0, width= symbol specifies the type of symbol with which to mark the object in the image. The following values are allowed: - - "circle" - - "box" - - "cross" - - "diamond" + + - "circle" + - "box" + - "cross" + - "diamond" size specifies the diameter in arcsec of the symbol (if plotSymbol == "circle"), or the width of the box in arcsec (if plotSymbol == "box") diff --git a/docs/conf.py b/docs/conf.py index 1ee2b64..9e7a233 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -30,12 +30,13 @@ # ones. extensions = [ 'sphinx.ext.autodoc', - 'sphinx_epytext', 'sphinx.ext.mathjax', 'sphinx.ext.napoleon', #'readthedocs_ext.readthedocs', ] +autodoc_mock_imports = ['PyWCSTools._wcs', 'PyWCSTools.wcs'] + # Add any paths that contain templates here, relative to this directory. templates_path = ['_templates'] From ea8cdf57a2e63492283431e9b64603641c5bdeee Mon Sep 17 00:00:00 2001 From: Matt Hilton Date: Wed, 27 May 2026 04:09:26 +0200 Subject: [PATCH 3/5] Add Sphinx cross-reference links throughout docstrings Replace plain function/class name mentions in docstrings with :func:, :meth:, :class:, and :data: roles so Sphinx generates hyperlinks in the HTML output. Covers cross-module references (e.g. astImages functions cited in astPlots) and intra-module references (e.g. fitSEDDict <-> mags2SEDDict <-> makeModelSEDDictList in astSED). Co-Authored-By: Claude Sonnet 4.6 --- astLib/astCoords.py | 2 +- astLib/astImages.py | 14 +++++------ astLib/astPlots.py | 14 +++++------ astLib/astSED.py | 60 ++++++++++++++++++++++----------------------- astLib/astStats.py | 2 +- astLib/astWCS.py | 2 +- 6 files changed, 47 insertions(+), 47 deletions(-) diff --git a/astLib/astCoords.py b/astLib/astCoords.py index d750983..0f252c6 100755 --- a/astLib/astCoords.py +++ b/astLib/astCoords.py @@ -361,7 +361,7 @@ def calcRADecSearchBox(RADeg, decDeg, radiusSkyDeg): """Calculates minimum and maximum RA, dec coords needed to define a box enclosing a circle of radius radiusSkyDeg around the given RADeg, decDeg coordinates. Useful for freeform queries of e.g. SDSS, UKIDSS etc.. Uses - calcAngSepDeg, so has the same limitations. + :func:`~astLib.astCoords.calcAngSepDeg`, so has the same limitations. Args: RADeg (float): RA coordinate of centre of search region diff --git a/astLib/astImages.py b/astLib/astImages.py index d0b60de..9c18643 100644 --- a/astLib/astImages.py +++ b/astLib/astImages.py @@ -37,10 +37,10 @@ def clipImageSectionWCS(imageData, imageWCS, RADeg, decDeg, clipSizeDeg, returnW coordinates in the original image corresponding to the clipped section. Note that the clip size is specified in degrees on the sky. For projections that have varying - real pixel scale across the map (e.g. CEA), use clipUsingRADecCoords instead. + real pixel scale across the map (e.g. CEA), use :func:`~astLib.astImages.clipUsingRADecCoords` instead. Similarly, this routine will not work for a WCS that has polynomial distortion coefficients - in the header (e.g., CTYPE1 = 'RA---TAN-SIP' etc.) - again clipUsingRADecCoords can be used + in the header (e.g., CTYPE1 = 'RA---TAN-SIP' etc.) - again :func:`~astLib.astImages.clipUsingRADecCoords` can be used in such cases. Args: @@ -177,13 +177,13 @@ def clipRotatedImageSectionWCS(imageData, imageWCS, RADeg, decDeg, clipSizeDeg, The resulting clip is rotated and/or flipped such that North is at the top, and East appears at the left. An updated WCS for the clipped section is also returned. Note that the alignment of the rotated WCS is currently not perfect - however, it is probably good enough in most - cases for use with ImagePlot for plotting purposes. + cases for use with :class:`astLib.astPlots.ImagePlot` for plotting purposes. Note that the clip size is specified in degrees on the sky. For projections that have varying - real pixel scale across the map (e.g. CEA), use clipUsingRADecCoords instead. + real pixel scale across the map (e.g. CEA), use :func:`~astLib.astImages.clipUsingRADecCoords` instead. Similarly, this routine will not work for a WCS that has polynomial distortion coefficients - in the header (e.g., CTYPE1 = 'RA---TAN-SIP' etc.) - again clipUsingRADecCoords can be used + in the header (e.g., CTYPE1 = 'RA---TAN-SIP' etc.) - again :func:`~astLib.astImages.clipUsingRADecCoords` can be used in such cases. Args: @@ -467,7 +467,7 @@ def scaleImage(imageData, imageWCS, scaleFactor): #--------------------------------------------------------------------------------------------------- def intensityCutImage(imageData, cutLevels): """Creates a matplotlib.pylab plot of an image array with the specified cuts in intensity - applied. This routine is used by saveBitmap and saveContourOverlayBitmap, which both + applied. This routine is used by :func:`~astLib.astImages.saveBitmap` and :func:`~astLib.astImages.saveContourOverlayBitmap`, which both produce output as .png, .jpg, etc. images. Args: @@ -1066,7 +1066,7 @@ def normalise(inputArray, clipMinMax): [clipMin, clipMax]. Used for normalising image arrays so that they can be turned into RGB arrays that matplotlib - can plot (see astPlots.ImagePlot). + can plot (see :class:`astLib.astPlots.ImagePlot`). Args: inputArray (numpy.ndarray): image data array diff --git a/astLib/astPlots.py b/astLib/astPlots.py index f4937af..7067072 100755 --- a/astLib/astPlots.py +++ b/astLib/astPlots.py @@ -84,12 +84,12 @@ class ImagePlot: associated WCS. Objects within the image boundaries can be marked by passing their WCS coordinates to - addPlotObjects. + :meth:`~astLib.astPlots.ImagePlot.addPlotObjects`. - Other images can be overlaid using addContourOverlay. + Other images can be overlaid using :meth:`~astLib.astPlots.ImagePlot.addContourOverlay`. For images rotated with North at the top, East at the left (as can be done using - astImages.clipRotatedImageSectionWCS or astImages.resampleToTanProjection), WCS coordinate + :func:`~astLib.astImages.clipRotatedImageSectionWCS` or :func:`~astLib.astImages.resampleToTanProjection`), WCS coordinate axes can be plotted, with tick marks set appropriately for the image size. Otherwise, a compass can be plotted showing the directions of North and East in the image. @@ -104,7 +104,7 @@ def __init__(self, imageData, imageWCS, axes = [0.1,0.1,0.8,0.8], \ colorBar = False, interpolation = "bilinear"): """Makes an ImagePlot from the given image array and astWCS. For coordinate axes to work, the image and WCS should have been rotated such that East is at the left, North is at the top - (see e.g. astImages.clipRotatedImageSectionWCS, or astImages.resampleToTanProjection). + (see e.g. :func:`~astLib.astImages.clipRotatedImageSectionWCS`, or :func:`~astLib.astImages.resampleToTanProjection`). If imageData is given as a list in the format [r, g, b], a color RGB plot will be made. However, in this case the cutLevels must be specified manually for each component as a list - @@ -117,7 +117,7 @@ def __init__(self, imageData, imageWCS, axes = [0.1,0.1,0.8,0.8], \ or decTickSteps are set to "auto", the appropriate axis scales will be determined automatically from the size of the image array and associated WCS. The tick step sizes can be overidden. If the coordinate axes are in sexagesimal format a dictionary in the format {'deg', 'unit'} is - needed (see RA_TICK_STEPS and DEC_TICK_STEPS for examples). If the coordinate axes are in + needed (see :data:`~astLib.astPlots.RA_TICK_STEPS` and :data:`~astLib.astPlots.DEC_TICK_STEPS` for examples). If the coordinate axes are in decimal format, the tick step size is specified simply in RA, dec decimal degrees. Args: @@ -358,7 +358,7 @@ def draw(self): def addContourOverlay(self, contourImageData, contourWCS, tag, levels = ["linear", "min", "max", 5], width = 1, color = "white", smooth = 0, highAccuracy = False): """Adds image data to the ImagePlot as a contour overlay. The contours can be removed using - removeContourOverlay. If a contour overlay already exists with this tag, it will be replaced. + :meth:`~astLib.astPlots.ImagePlot.removeContourOverlay`. If a contour overlay already exists with this tag, it will be replaced. Args: contourImageData (numpy.ndarray): image data array from which contours are to be generated @@ -440,7 +440,7 @@ def addPlotObjects(self, objRAs, objDecs, tag, symbol="circle", size=4.0, width= color can be any valid matplotlib color (e.g. "red", "green", etc.) - The objects can be removed from the plot by using removePlotObjects(), and then calling + The objects can be removed from the plot by using :meth:`~astLib.astPlots.ImagePlot.removePlotObjects`, and then calling draw(). If the ImagePlot already has a set of plotObjects with the same tag, they will be replaced. diff --git a/astLib/astSED.py b/astLib/astSED.py index f5afd29..b94746f 100644 --- a/astLib/astSED.py +++ b/astLib/astSED.py @@ -190,9 +190,9 @@ class SED: To create a SED object, lists (or np arrays) of wavelength and relative flux must be provided. The SED can optionally be redshifted. The wavelength units of SEDs are assumed to be Angstroms - flux - calculations using Passband and SED objects specified with different wavelength units will be incorrect. + calculations using :class:`~astLib.astSED.Passband` and SED objects specified with different wavelength units will be incorrect. - The StellarPopulation class (and derivatives) can be used to extract SEDs for specified ages from e.g. + The :class:`~astLib.astSED.StellarPopulation` class (and derivatives) can be used to extract SEDs for specified ages from e.g. the Bruzual & Charlot 2003 or Maraston 2005 models. """ @@ -409,7 +409,7 @@ def normaliseToMag(self, ABMag, passband): Args: ABMag (float): AB magnitude to which the SED is to be normalised at the given passband - passband (Passband): passband at which normalisation to AB magnitude is calculated + passband (:class:`~astLib.astSED.Passband`): passband at which normalisation to AB magnitude is calculated """ @@ -425,7 +425,7 @@ def matchFlux(self, matchSED, minWavelength, maxWavelength): flux in the same region in matchSED. Useful for plotting purposes. Args: - matchSED (SED): SED to match flux to + matchSED (:class:`~astLib.astSED.SED`): SED to match flux to minWavelength (float): minimum of range in which to match flux of current SED to matchSED maxWavelength (float): maximum of range in which to match flux of current SED to matchSED @@ -446,7 +446,7 @@ def calcFlux(self, passband): """Calculates flux in the given passband. Args: - passband (Passband): filter passband through which to calculate the flux from the SED + passband (:class:`~astLib.astSED.Passband`): filter passband through which to calculate the flux from the SED Returns: float: flux @@ -474,7 +474,7 @@ def calcMag(self, passband, addDistanceModulus = True, magType = "Vega"): in Mpc at the redshift of the SED) is added. Args: - passband (Passband): filter passband through which to calculate the magnitude from the SED + passband (:class:`~astLib.astSED.Passband`): filter passband through which to calculate the magnitude from the SED addDistanceModulus (bool): if True, adds 5.0*log10*(dl*1e5) to the mag returned, where dl is the luminosity distance (Mpc) corresponding to the SED z magType (str): either "Vega" or "AB" @@ -505,8 +505,8 @@ def calcColour(self, passband1, passband2, magType = "Vega"): """Calculates the colour passband1-passband2. Args: - passband1 (Passband): filter passband through which to calculate the first magnitude - passband2 (Passband): filter passband through which to calculate the second magnitude + passband1 (:class:`~astLib.astSED.Passband`): filter passband through which to calculate the first magnitude + passband2 (:class:`~astLib.astSED.Passband`): filter passband through which to calculate the second magnitude magType (str): either "Vega" or "AB" Returns: @@ -522,10 +522,10 @@ def calcColour(self, passband1, passband2, magType = "Vega"): def getSEDDict(self, passbands): """This is a convenience function for pulling out fluxes from a SED for a given set of passbands - in the same format as made by mags2SEDDict - designed to make fitting code simpler. + in the same format as made by :func:`~astLib.astSED.mags2SEDDict` - designed to make fitting code simpler. Args: - passbands (list): list of Passband objects through which fluxes will be calculated + passbands (list): list of :class:`~astLib.astSED.Passband` objects through which fluxes will be calculated """ @@ -681,8 +681,8 @@ class StellarPopulation: files, containing columns for age, wavelength, and flux. Columns are counted from 0 ... n. Lines starting with # are ignored. - The classes M05Model (for Maraston 2005 models), BC03Model (for Bruzual & Charlot 2003 models), and - P09Model (for Percival et al. 2009 models) are derived from this class. The only difference between + The classes :class:`~astLib.astSED.M05Model` (for Maraston 2005 models), :class:`~astLib.astSED.BC03Model` (for Bruzual & Charlot 2003 models), and + :class:`~astLib.astSED.P09Model` (for Percival et al. 2009 models) are derived from this class. The only difference between them is the code used to load in the model data. """ @@ -755,8 +755,8 @@ def getColourEvolution(self, passband1, passband2, zFormation, zStepSize = 0.05, StellarPopulation with redshift, from z = 0 to z = zFormation. Args: - passband1 (Passband): filter passband through which to calculate the first magnitude - passband2 (Passband): filter passband through which to calculate the second magnitude + passband1 (:class:`~astLib.astSED.Passband`): filter passband through which to calculate the first magnitude + passband2 (:class:`~astLib.astSED.Passband`): filter passband through which to calculate the second magnitude zFormation (float): formation redshift of the StellarPopulation zStepSize (float): size of interval in z at which to calculate model colours magType (str): either "Vega" or "AB" @@ -789,7 +789,7 @@ def getMagEvolution(self, passband, magNormalisation, zNormalisation, zFormation (apparent) at z = zNormalisation. Args: - passband (Passband): filter passband through which to calculate the magnitude + passband (:class:`~astLib.astSED.Passband`): filter passband through which to calculate the magnitude magNormalisation (float): sets the apparent magnitude of the SED at zNormalisation zNormalisation (float): the redshift at which the magnitude normalisation is carried out zFormation (float): formation redshift of the StellarPopulation @@ -843,7 +843,7 @@ def calcEvolutionCorrection(self, zFrom, zTo, zFormation, passband, magType = "V zFrom (float): redshift to evolution correct from zTo (float): redshift to evolution correct to zFormation (float): formation redshift of the StellarPopulation - passband (Passband): filter passband through which to calculate magnitude + passband (:class:`~astLib.astSED.Passband`): filter passband through which to calculate magnitude magType (str): either "Vega" or "AB" Returns: @@ -1052,8 +1052,8 @@ def __init__(self, fileName): #------------------------------------------------------------------------------------------------------------ def makeModelSEDDictList(modelList, z, passbandsList, labelsList = [], EBMinusVList = [0.0], forceYoungerThanUniverse = True): - """This routine makes a list of SEDDict dictionaries (see mags2SEDDict) for fitting using - fitSEDDict. This speeds up the fitting as this allows us to calculate model SED magnitudes only once, + """This routine makes a list of SEDDict dictionaries (see :func:`~astLib.astSED.mags2SEDDict`) for fitting using + :func:`~astLib.astSED.fitSEDDict`. This speeds up the fitting as this allows us to calculate model SED magnitudes only once, if all objects to be fitted are at the same redshift. We add some meta data to the modelSEDDicts (e.g. the model file names). @@ -1064,15 +1064,15 @@ def makeModelSEDDictList(modelList, z, passbandsList, labelsList = [], EBMinusVL included. Args: - modelList (list): list of StellarPopulation models to include + modelList (list): list of :class:`~astLib.astSED.StellarPopulation` models to include z (float): redshift to apply to all stellar population models in modelList - passbandsList (list): list of Passband objects + passbandsList (list): list of :class:`~astLib.astSED.Passband` objects labelsList (list): optional list used for labelling passbands in output SEDDicts EBMinusVList (list): list of E(B-V) extinction values to apply to all models, in magnitudes forceYoungerThanUniverse (bool): if True, do not allow models that exceed the age of the universe at z Returns: - list: list of dictionaries containing model fluxes, to be used as input to fitSEDDict + list: list of dictionaries containing model fluxes, to be used as input to :func:`~astLib.astSED.fitSEDDict` """ @@ -1105,16 +1105,16 @@ def makeModelSEDDictList(modelList, z, passbandsList, labelsList = [], EBMinusVL #------------------------------------------------------------------------------------------------------------ def fitSEDDict(SEDDict, modelSEDDictList): - """Fits the given SED dictionary (made using mags2SEDDict) with the given list of model SED - dictionaries. The latter should be made using makeModelSEDDictList, and entries for fluxes should + """Fits the given SED dictionary (made using :func:`~astLib.astSED.mags2SEDDict`) with the given list of model SED + dictionaries. The latter should be made using :func:`~astLib.astSED.makeModelSEDDictList`, and entries for fluxes should correspond directly between the model and SEDDict. Returns a dictionary with best fit values. Args: - SEDDict (dict): dictionary of observed fluxes and uncertainties, in format of mags2SEDDict + SEDDict (dict): dictionary of observed fluxes and uncertainties, in format of :func:`~astLib.astSED.mags2SEDDict` modelSEDDictList (list): list of dictionaries containing fluxes of models to be fitted to the - observed fluxes listed in the SEDDict, made using makeModelSEDDictList + observed fluxes listed in the SEDDict, made using :func:`~astLib.astSED.makeModelSEDDictList` Returns: dict: results of the fitting - keys: @@ -1166,17 +1166,17 @@ def mags2SEDDict(ABMags, ABMagErrs, passbands): """Takes a set of corresponding AB magnitudes, uncertainties, and passbands, and returns a dictionary with keys 'flux', 'fluxErr', 'wavelength'. Fluxes are in units of erg/s/cm^2/Angstrom, wavelength in Angstroms. These dictionaries are the staple diet of the - fitSEDDict routine. + :func:`~astLib.astSED.fitSEDDict` routine. Args: ABMags (list or numpy.ndarray): AB magnitudes, specified in corresponding order to passbands and ABMagErrs ABMagErrs (list or numpy.ndarray): AB magnitude errors, specified in corresponding order to passbands and ABMags - passbands (list): list of Passband objects, specified in corresponding order to ABMags and ABMagErrs + passbands (list): list of :class:`~astLib.astSED.Passband` objects, specified in corresponding order to ABMags and ABMagErrs Returns: - dict: dictionary with keys ``{'flux', 'fluxErr', 'wavelength'}``, suitable for input to fitSEDDict + dict: dictionary with keys ``{'flux', 'fluxErr', 'wavelength'}``, suitable for input to :func:`~astLib.astSED.fitSEDDict` """ @@ -1203,7 +1203,7 @@ def mag2Flux(ABMag, ABMagErr, passband): Args: ABMag (float): magnitude on AB system in passband ABMagErr (float): uncertainty in AB magnitude in passband - passband (Passband): Passband object at which ABMag was measured + passband (:class:`~astLib.astSED.Passband`): Passband object at which ABMag was measured Returns: list: [flux, fluxError], in units of erg/s/cm^2/Angstrom @@ -1228,7 +1228,7 @@ def flux2Mag(flux, fluxErr, passband): Args: flux (float): flux in erg/s/cm^2/Angstrom in passband fluxErr (float): uncertainty in flux in passband, in erg/s/cm^2/Angstrom - passband (Passband): Passband object at which ABMag was measured + passband (:class:`~astLib.astSED.Passband`): Passband object at which ABMag was measured Returns: list: [ABMag, ABMagError], in AB magnitudes diff --git a/astLib/astStats.py b/astLib/astStats.py index f99a8fe..26ae831 100755 --- a/astLib/astStats.py +++ b/astLib/astStats.py @@ -504,7 +504,7 @@ def weightedLSFit(dataList, weightType): weightType (str): if ``"errors"``, weights are calculated assuming the input data is in the format [x, y, error on y]; if ``"weights"``, the weights are assumed to be already calculated and stored in a fourth column [x, y, error on y, weight] (as used by e.g. - biweightLSFit) + :func:`~astLib.astStats.biweightLSFit`) Returns: dict or None: slope and intercept on y-axis with associated errors in the format diff --git a/astLib/astWCS.py b/astLib/astWCS.py index b944ef4..f2ae541 100644 --- a/astLib/astWCS.py +++ b/astLib/astWCS.py @@ -79,7 +79,7 @@ def __init__(self, headerSource, extensionName = 0, mode = "image", zapKeywords zapKeywords (list): keywords to remove from the header before making the WCS object useAstropyWCS (bool): if True, use astropy.wcs to perform WCS - coordinate conversions in wcs2pix, pix2wcs (if False, use + coordinate conversions in :meth:`~astLib.astWCS.WCS.wcs2pix`, :meth:`~astLib.astWCS.WCS.pix2wcs` (if False, use PyWCSTools) naxis (int): number of WCS axes to use From 46e713525b393a04c898776e40829d483e9b5b10 Mon Sep 17 00:00:00 2001 From: Matt Hilton Date: Wed, 27 May 2026 04:18:59 +0200 Subject: [PATCH 4/5] Add module links in docs index and astStats to API reference Replace plain module names in the README bullet list with :mod: roles so the Sphinx-rendered index page links directly to each module's API reference section. Also adds the missing astStats section to reference.rst. Co-Authored-By: Claude Sonnet 4.6 --- README.rst | 14 +++++++------- docs/reference.rst | 7 +++++++ 2 files changed, 14 insertions(+), 7 deletions(-) diff --git a/README.rst b/README.rst index 10c41a5..2cca61d 100644 --- a/README.rst +++ b/README.rst @@ -11,13 +11,13 @@ **astLib** is divided into several modules: -* astCalc (general calculations, e.g. luminosity distance etc.) -* astCoords (coordinate conversions etc.) -* astImages (clip sections from .fits etc.) -* astPlots (provides a flexible image plot class, e.g. plot image with catalogue objects overlaid) -* astSED (calculate colours, magnitudes from stellar population models or spectral templates, fit photometric observations using stellar population models etc.) -* astStats (statistics, e.g. biweight location/scale estimators etc.) -* astWCS (routines for using FITS World Coordinate System information) +* :mod:`astLib.astCalc` (general calculations, e.g. luminosity distance etc.) +* :mod:`astLib.astCoords` (coordinate conversions etc.) +* :mod:`astLib.astImages` (clip sections from .fits etc.) +* :mod:`astLib.astPlots` (provides a flexible image plot class, e.g. plot image with catalogue objects overlaid) +* :mod:`astLib.astSED` (calculate colours, magnitudes from stellar population models or spectral templates, fit photometric observations using stellar population models etc.) +* :mod:`astLib.astStats` (statistics, e.g. biweight location/scale estimators etc.) +* :mod:`astLib.astWCS` (routines for using FITS World Coordinate System information) The astWCS module is a higher level interface to PyWCSTools, a simple SWIG (http://www.swig.org) wrapping of some of the routines from WCSTools by Jessica Mink (http://tdc-www.harvard.edu/software/wcstools/). It is diff --git a/docs/reference.rst b/docs/reference.rst index bb7d107..a9bd0e4 100644 --- a/docs/reference.rst +++ b/docs/reference.rst @@ -48,3 +48,10 @@ astSED :members: +astStats +-------- + +.. automodule:: astLib.astStats + :members: + + From 1bc7c50f5be86af446d781249006b171696f0a8f Mon Sep 17 00:00:00 2001 From: Matt Hilton Date: Wed, 27 May 2026 04:27:08 +0200 Subject: [PATCH 5/5] Format inline URLs as named hyperlinks in README Replaces bare URLs in parentheses with RST named links so the URL is hidden in rendered output (docs, GitHub). Co-Authored-By: Claude Sonnet 4.6 --- README.rst | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/README.rst b/README.rst index 2cca61d..bc4bb1d 100644 --- a/README.rst +++ b/README.rst @@ -19,13 +19,14 @@ * :mod:`astLib.astStats` (statistics, e.g. biweight location/scale estimators etc.) * :mod:`astLib.astWCS` (routines for using FITS World Coordinate System information) -The astWCS module is a higher level interface to PyWCSTools, a simple SWIG (http://www.swig.org) wrapping -of some of the routines from WCSTools by Jessica Mink (http://tdc-www.harvard.edu/software/wcstools/). It is +The astWCS module is a higher level interface to PyWCSTools, a simple `SWIG `_ wrapping +of some of the routines from `WCSTools `_ by Jessica Mink. It is used by some routines in astCoords, astImages and astPlots. -The goal of **astLib** was to provide features useful to astronomers that are not included in the scipy -(http://scipy.org), numpy (http://numpy.scipy.org) or matplotlib (http://matplotlib.sourceforge.net) modules -on which astLib depends. For a far more extensive set of Python astronomy modules, see astropy -(http://www.astropy.org/). +The goal of **astLib** was to provide features useful to astronomers that are not included in the +`scipy `_, `numpy `_ or +`matplotlib `_ modules +on which astLib depends. For a far more extensive set of Python astronomy modules, see +`astropy `_. Some scripts using **astLib** can be found in the examples/ folder provided with the source code distribution.