#!/usr/bin/env python
#----------------------------------------------------------------------
# FILE: wcsgrat.py
# PURPOSE: Provide classes and methods to produce graticule data
# for WCS labeling and plotting grid lines.
# AUTHOR: M.G.R. Vogelaar, University of Groningen, The Netherlands
# DATE: August 17, 2008
# UPDATES: Version 0.2, September 09, 2008
# Version 0.3, October 17, 2008
# Version 1.0, June 17, 2009
# Version 1.1, December 05, 2012 (bug fixes only)
# Version 1.2, December 14, 2014; replaced == None etc.
#
# VOG, 06-07-2023: Replaced depredated MPL's get_geometry with
# appropriate routines.
#
# (C) University of Groningen
# Kapteyn Astronomical Institute
# Groningen, The Netherlands
# E: gipsy@astro.rug.nl
#----------------------------------------------------------------------
"""
Module wcsgrat
==============
A graticule is a system of crossing lines on a map representing
positions of which one coordinate is constant.
For a spatial map it consists of parallels of latitude and
meridians of longitude as defined by a given projection.
This module is used to set up such graticules and labels for the
selected world coordinate system. It plots the results with
plotting library
`Matplotlib <http://matplotlib.sourceforge.net/index.html>`_.
Besides spatial axes, it supports also spectral axes and a mix of both
(e.g. position-velocity diagrams). It deals with data dimensions > 2
by allowing arbitrary selections of two axes.
The transformations between pixel coordinates and world coordinates
are based on module
:mod:`wcs` which is a Python binding for Mark R. Calabretta's library
`WCSLIB <http://www.atnf.csiro.au/people/mcalabre/WCS>`_.
>From *WCSLIB* we use only the core transformation routines.
Header parsing is done with module :mod:`wcs`.
Axes types that are not recognized by this software is treated as being linear.
The axes types correspond with keywords *CTYPEn* in a FITS file.
The information from a FITS file is retrieved by module
`PyFITS <http://www.stsci.edu/resources/software_hardware/pyfits>`_
.. seealso:: Tutorial material with code examples:
* Tutorial maputils module
which contains many examples with source code,
see :ref:`maputils_tutorial`.
* Figure gallery 'all sky plots'
with many examples of Graticule constructors,
see :ref:`allsky_tutorial`.
.. moduleauthor:: Martin Vogelaar <gipsy@astro.rug.nl>
Module level data
-----------------
:data:`left, bottom, right, top`
The variables *left*, *bottom*, *right* and *top* are
equivalent to the strings *"left"*, *"bottom"*, *"right"* and *"top"*
and are used as identifiers for plot axes.
:data:`native, notnative, bothticks, noticks`
The variables *native*, *notnative*, *bothticks*, *noticks*
correspond to the numbers 0, 1, 2 and 3 and represent modes
to make ticks along an axis visible or invisible. Ticks along an axis
can represent both world coordinate types (e.g. when a map is rotated). Sometimes
one wants to allow this and sometimes not.
.. tabularcolumns:: |p{25mm}|p{125mm}|
========= =============================================
Tick mode Description
========= =============================================
native Show only ticks that are native to the
coordinate axis. Do not allow ticks
that correspond to the axis for which
a constant value applies. So, for example,
in a RA-DEC
map which is rotated 45 degrees we want only
Right Ascensions along the x-axis.
notnative Plot the ticks that are not native to the
coordinate axis. So, for example, in a RA-DEC
map which is rotated 45 degrees we want only
Declinations along the x-axis.
bothticks Allow both type of ticks along a plot axis
noticks Do not allow any tick to be plotted.
========= =============================================
Functions
---------
.. index:: Convert degrees to separate field of the hms/dms format
.. autofunction:: gethmsdms
.. index:: Convert hms or dms values to text or LaTeX representation
.. autofunction:: makelabel
Class Graticule
---------------
.. autoclass:: Graticule
.. autoclass:: WCStick
Class Insidelabels
--------------------
.. autoclass:: Insidelabels
"""
from sys import stdout
from kapteyn import wcs # The Kapteyn Python binding to WCSlib, including celestial transformations
from kapteyn.positions import str2pos, parsehmsdms, unitfactor
from kapteyn import rulers
from kapteyn.celestial import skyparser
from string import ascii_letters as letters
from random import choice
import numpy
__version__ = '1.1'
(left,bottom,right,top) = list(range(4)) # Names of the four plot axes
(native, notnative, bothticks, noticks) = list(range(4))
# A variable to change the vertical size of a TeX box
# For more information, have a look at method
# _update_metrics() at the end of this code.
tweakhms = 0
def issequence(obj):
return isinstance(obj, (list, tuple, numpy.ndarray))
def parseplotaxes(plotaxes):
#-----------------------------------------------------------
"""
Purpose: It is possible to specify axes by an integer
or by a string.
The function is used internally to allow
flexible input of numbers to identify one of
the four plot axes.
Parameters:
plotaxes - Scalar or sequence with elements that are either
integers or strings or a combination of those.
Returns: A list with unique numbers between 0 and 3
Notes: - Order is unimportant
- Input can be a scalar or a sequence (tuple, list)
- Scalers are upgraded to a list.
- The result has only unique numbers
"""
#-----------------------------------------------------------
if not issequence(plotaxes):
plotaxes = [plotaxes,]
if isinstance(plotaxes, tuple):
plotaxes = list(plotaxes)
for i,pa in enumerate(plotaxes):
if isinstance(pa, str):
if "LEFT".find(pa.upper()) == 0:
plotaxes[i] = left
elif "BOTTOM".find(pa.upper()) == 0:
plotaxes[i] = bottom
elif "RIGHT".find(pa.upper()) == 0:
plotaxes[i] = right
elif "TOP".find(pa.upper()) == 0:
plotaxes[i] = top
else:
raise ValueError("[%s] Cannot identify this plot axis!" % pa)
for pa in plotaxes: # Check validity
if pa < 0 or pa > 3:
raise ValueError("Cannot identify this plot axis!")
aset = {}
list(map(aset.__setitem__, plotaxes, [None] * len(plotaxes)))
return list(aset.keys())
def parsetickmode(tickmode):
#-----------------------------------------------------------
"""
Purpose:
Parameters:
tickmode - Scalar or sequence with elements that are either
integers or strings or a combination of those.
Returns: A list with unique numbers between 0 and 3
Notes:
"""
#-----------------------------------------------------------
#(native, notnative, bothticks, noticks) = list(range(4))
if isinstance(tickmode, str):
if "NATIVE_TICKS".find(tickmode.upper()) == 0:
tickmode = native
elif "SWITCHED_TICKS".find(tickmode.upper()) == 0:
tickmode = notnative
elif "ALL_TICKS".find(tickmode.upper()) == 0:
tickmode = bothticks
elif "NO_TICKS".find(tickmode.upper()) == 0:
tickmode = noticks
else:
raise ValueError("[%s] Cannot identify this tick mode!" % tickmode)
if tickmode < 0 or tickmode > 3:
raise ValueError("%d does not correspond to a supported tick mode!" % tickmode)
return tickmode
[docs]def gethmsdms(a, prec, axtype, skysys, eqlon=None):
#---------------------------------------------------------------------------
"""
Given a number in degrees and an axis type in *axtype*
equal to 'longitude' or 'latitude',
calculate and return the parts of its sexagesimal representation, i.e.
hours or degrees, minutes and seconds. Also return the fractional seconds
and the sign if the input was a value at negative latitude.
The value for *skysys* sets the formatting to hours/minutes/seconds
if it represents an equatorial system.
:param a:
The longitude or latitude in degrees.
:type a:
Floating point
:param prec:
The number of decimals in the seconds
:type prec:
Integer
:param axtype:
One of 'longitude' or 'latitude'
:type axtype:
String
:param skysys:
The sky system
:type skysys:
Integer
:Returns:
tuple: (Ihours, Ideg, Imin, Isec, Fsec, sign)
which represent Integer values for the hours, degrees, minutes
and seconds. *Fsec* is the fractional part of the seconds.
Element *sign* is -1 for negative latitudes and +1 for positive
latitudes.
"""
#---------------------------------------------------------------------------
if not axtype in ['longitude', 'latitude']:
return None
sign = 1
if not issequence(skysys): # Usually a number, sometimes a tuple
skysys = [skysys]
if eqlon is None:
eqlon = (axtype == 'longitude') and (wcs.equatorial in skysys)
if axtype == 'longitude':
degs = numpy.fmod(a, 360.0) # Now in range -360, 360
if degs < 0.0:
degs += 360.0
else:
if a > 90.0 or a < -90.0:
return None
if a < 0.0:
sign = -1
degs = sign * a # Make positive
if prec < 0:
prec = 0
# How many seconds is this? Round to 'prec'
if eqlon:
sec = numpy.round(degs*240.0, prec)
sec = numpy.fmod(sec, 360.0*240.0) # Rounding can result in 360 deg. again, so correct
else:
sec = numpy.round(degs*3600.0, prec)
AIsec = int(sec) # Integer seconds
Fsec = sec - AIsec # Fractional remainder
X = AIsec / 3600.0 # Either hours or degrees
IX = int(X)
secleft = AIsec - IX*3600.0
Imin = int(secleft/60.0)
Isec = secleft - Imin*60.0
if eqlon:
Ihours = IX; Ideg = None
else:
Ideg = IX; Ihours = None
return Ihours, Ideg, Imin, Isec, Fsec, sign
def makeexplab(fmt, x):
#-------------------------------------------------------------------------------
"""
Format exponentials
"""
#-------------------------------------------------------------------------------
# The format string contains two formats.
if 1: #try:
c = fmt.rfind('%')
fmt1 = fmt[:c]
fmt2 = fmt[c:]
n = fmt2[1:fmt2.index('e')]
ex = int(n)
val = x/10.**ex
lab = fmt1%val+"\cdot10^{%s}"%(str(ex)) + fmt2[fmt2.index('e')+1:]
else: #except:
raise TypeError("Wrong format (%s) for numbers in exponential LaTeX format"%fmt)
return lab
[docs]def makelabel(hmsdms, Hlab, Dlab, Mlab, Slab, prec, fmt, tex):
#-------------------------------------------------------------------------------
"""
From the output of function *gethmsdms* and some Booleans, this function
creates a label in plain text or in TeX. The Booleans set a flag whether
a field (hours, degrees, minutes or seconds) should be printed or not.
The *fmt* parameter is used if it does not contain the percentage character (%)
but instead contains characters from the set HDMS. A capital overules the
corresponding Boolean value, so if *fmt='HMS'*, the values for *Hlab*, *Mlab* and
*Slab* are all set to True.
:param hmsdms:
The output of function :func:`gethmsdms`
:type hmsdms:
Tuple with integer and floating point numbers
:param Hlab:
If False, there is no need to print the hours
:param Dlab:
If False, there is no need to print the degrees
:param Mlab:
If False, there is no need to print the minutes
:param Slab:
If False, there is no need to print the seconds
:param fmt:
String containing a combination of the characters
['H', 'D', 'M', 'S', '.', 'h', 'd', 'm', 's']
A capital sets the corresponding input Boolean (Hlab, Dlab, etc.)
to True. A dot starts to set the precision. The number of characters
after the dot set the precision itself.
A character that is not a capital sets the corresponding input
Boolean (Hlab, Dlab, etc.) to False. This is a bit dangerous because
with this option one can suppress fields to be printed that contain
a value unequal to zero. It is applied if you want to suppress
e.g. seconds if all the seconds in your label are 0.0.
The suppression of printing minutes is overruled if hours (or degrees)
and seconds are required. Otherwise we could end up with non
standard labels (e.g. 2h30s).
:type fmt:
String
:param tex:
If True, then format the labels in LaTeX.
:type tex:
Boolean
:Returns:
*lab*, a label in either hms or dms in plain text or in LaTeX format.
:Examples:
>>> # Set the format in Hours, minutes and seconds with a precision
>>> # of three. The suppression of minutes will not work here:
>>> grat.setp_tick(wcsaxis=0, fmt="HmS.SSS")
>>> # The same effect is obtained with:
>>> grat.setp_tick(wcsaxis=0, fmt="HmS.###")
>>> # Let the system determine whether seconds are printed
>>> # but make sure that degrees and minutes are included:
>>> grat.setp_tick(wcsaxis=1, fmt="DM")
>>> # If we know that all minutes and seconds in our labels are 0.0
>>> # and we want only the hours to be printed, then use:
>>> grat.setp_tick(wcsaxis=0, fmt="Hms")
>>> grat.setp_tick(wcsaxis=0, fmt="Dms")
>>> # Plot labels in Degrees even if the axis is an equatorial longitude.
"""
#-------------------------------------------------------------------------------
Ihours, Ideg, Imin, Isec, Fsec, sign = hmsdms
hms = False
if Ihours is not None:
hms = True
if not (fmt is None):
if fmt.find('%') == -1: # Not a common format, must be a HMS/DMS format
if fmt.find('H') != -1:
Hlab = True
if fmt.find('D') != -1:
Dlab = True
if fmt.find('M') != -1:
Mlab = True
if fmt.find('S') != -1:
Slab = True
if fmt.find('h') != -1:
Hlab = False
if fmt.find('d') != -1:
Dlab = False
if fmt.find('m') != -1:
Mlab = False
if fmt.find('s') != -1:
Slab = False
s2 = fmt.split('.')
if len(s2) > 1:
prec = len(s2[1])
lab = ""
# Minutes must be printed if hours and seconds are required
# otherwise one can end up with a minimal label with hours/degs and
# no minuts in between
if hms:
if Hlab:
if tex:
lab += r"%d^{\rm h}"%Ihours
else:
lab += "%dh"%Ihours
if Slab:
Mlab = True
else:
if Dlab:
si = ' '
if sign == -1:
si = '-'
if tex:
lab += r"%c%d^{\circ}"%(si,Ideg)
else:
lab += "%c%dd"%(si,Ideg)
if Slab:
Mlab = True
if Mlab:
if tex:
if hms:
lab += r"%.2d^{\rm m}"%Imin
else:
lab += r"%.2d^{\prime}"%Imin
else:
lab += "%.2dm"%Imin
if Slab:
lab += "%.2d"%Isec
if prec > 0:
fsec = ".%*.*d" % (prec, prec, int(round(Fsec*10.0**prec,0)))
lab += fsec
if not tex:
lab += 's'
else:
if hms:
lab += r"^{\rm s}"
else:
lab += r"^{\prime\prime}"
if sign == -1 and not Dlab:
lab = "-"+lab
return lab
[docs]class WCStick(object):
#-------------------------------------------------------------------------------
"""
A WCStick object is an intersection of a parallel or meridian (or equivalent
lines with one constant world coordinate) with one of
the axes of a rectangle in pixels. The position of that intersection is
stored in pixel coordinates and can be used to plot a (formatted) label
showing the position of the constant world coordinate of the graticule line.
This class is only used in the context of the Graticule class.
"""
#-------------------------------------------------------------------------------
def __init__(self, x, y, axisnr, labval, wcsaxis, offset, fun=None, fmt=None):
#-----------------------------------------------------------
"""
Purpose: Store tick properties
Parameters:
x - pixel coordinate of position in x-direction
y - pixel coordinate of position in y-direction
axisnr - number between 0 and 4 representing
axes: left,bottom,right,top
labval - A (formatted) string representation
of the graticule line to which the tick belongs.
wcsaxis - 0 for the first WCS axis, 1 for the second WCS axis
offset - Was it an offset?
Returns: WCStick object with kwargs which sets the
plot properties for a tick
Notes: The keyword arguments attribute contains
information about plot properties.These
properties are (plot-)package specific.
We standardized on Matplotlib. If a plot
method is added for another plotting system,
one has to translate those properties to
the equivalents of that other system.
Have a look at the code in method 'plot' to
explore how such conversions should be done.
"""
#-----------------------------------------------------------
self.x = x # A tick is an intersection with a rectangle at x, y
self.y = y
self.axisnr = axisnr # Which axis did it intersect?
self.labval = labval # What is the value of the world coordinate?
self.offset = offset # Is it an offset?
self.wcsaxis = wcsaxis # To which (wcs) type of axis does it belong?
self.fmt = fmt # Python format string for conversion number to strings
self.fun = fun # Convert a tick (wcs) value with this function
self.axtype = None # Is it a long. or lat. axis?
self.skysys = None # And does it belong to an equatorial system?
self.label = "" # The text as it will appear in a plot
self.prec = None # The precision of the seconds in a hms/dms format
self.delta = None # The step between two ticks for this wcs axis
self.tex = None # Format label in TeX?
self.texsexa = None # Format the sexagesimal numbers
self.kwargs = {} # Keyword arguments to set plot attributes
self.markkwargs = {} # Keyword arguments to set plot attributes for the markers
def createlabels(Tlist):
#-------------------------------------------------------------------------------
"""
This function creates the text labels for a list of tick which belong to
the same axis and the same wcs axis.
"""
#-------------------------------------------------------------------------------
global tweakhms
Hprev = Dprev = Mprev = Sprev = signprev = None
dirprev = None; Labprev = None
H = D = M = S = True
first = True
for t in Tlist:
# There are some conditions for plotting labels in hms/dms:
stdout.flush()
if t.axtype in ['longitude', 'latitude'] and t.offset == False and t.fun is None and\
(t.fmt is None or t.fmt.find('%') == -1):
# Each tick label has its own fontsize. We take the first label's fontsize as
# the value for which we tweak the characters 's' and 'm' in the TeX labels
if 'fontsize' in t.kwargs and t.tex:
if tweakhms == 0:
tweakhms = t.kwargs['fontsize']
if t.axtype == 'longitude' and t.fmt is not None and 'D' in t.fmt.upper():
# This is an equatorial longitude axis for which one wants
# DMS formatting not HMS formatting. Useful for "all sky" plots
eqlon = False
else:
eqlon = None
hmsdms = gethmsdms(t.labval, t.prec, t.axtype, t.skysys, eqlon)
Ihours, Ideg, Imin, Isec, Fsec, sign = hmsdms
if first:
hmsdms1 = hmsdms
first = False
# Take care of the fact that the first label is not at the end
# so possibly the direction of labeling can be switched.
if Labprev is not None:
direction = t.labval > Labprev
if direction != dirprev and dirprev is not None:
Hprev, Dprev, Mpriv, Sprev, Fsecprev, sign = hmsdms1
Sprev += Fsecprev
dirprev = direction
H = (Ihours != Hprev); D = (Ideg != Dprev); M = (Imin != Mprev); S = (Isec+Fsec != Sprev)
Hprev = Ihours; Dprev = Ideg; Mprev = Imin; Sprev = Isec+Fsec; signprev = sign
Labprev = t.labval
texsexa = t.texsexa
if texsexa is None:
if t.tex:
texsexa = True
else:
texsexa = False
lab = makelabel(hmsdms, H, D, M, S, t.prec, t.fmt, texsexa)
else:
if t.fun is None:
val = t.labval
else:
val = t.fun(t.labval)
if t.fmt is None:
lab = "%g" % val
else:
if t.tex and t.fmt.count("%") == 2:
# There are two formats given so the user wants an exponential format
lab = makeexplab(t.fmt, val)
else:
try:
lab = t.fmt % val
except:
raise TypeError("Wrong format (%s) for numbers"%t.fmt)
if t.tex:
lab = r"$%s$" % lab
t.label = lab
class Gratline(object):
#-----------------------------------------------------------------
"""
This class is used to find a set of coordinates that defines
(part of) a graticule line for which one of the world coordinates
is a constant. It stores the coordinates of the intersections
with the box and a corresponding label for annotating purposes.
"""
#-----------------------------------------------------------------
def __init__(self, wcsaxis, constval, gmap, pxlims, pylims, wxlims, wylims,
linesamples, mixgrid, skysys, addx=None, addy=None, addpixels=False, offsetlabel=None,
fun=None, fmt=None):
#-----------------------------------------------------------------
"""
Purpose: Initialize a 'grid line' which is part of a
graticule. The method should be called within the
context of the Graticule class only.
Parameters:
wcsaxis: One of the values 0 or 1 for the first and second
world coordinate type. Or a number > 1 which is
the id of a graticule line representing a border.
constval: For a graticule line, one of the world coordinates
is a constant. The other coordinate varies within
certain limits. The value of the constant is given by
this parameter.
gmap: The projection object for these wcs axes.
pxlims: The lower and upper limits of an image in grids along
the x-axis.
pylims: Same for the y-axis
wxlims: Lower and upper values of the image in world
coordinates.
wylims: Same as wxlims but now for the y-axis
linesamples:This parameter sets the number of samples with
which we build a graticule line. In fact it is
equivalent to a step size for straight lines.
Therefore it is also used to identify jumps
in longitudes and latitudes.
We want to avoid jumps in a plot and therefore
break up the graticule line in pieces.
mixgrid: For images with only one spatial axis, we need
to supply a pixel value for the matching spatial
axis which is more or less hidden but essential
in the coordinate transformation.
skysys: The skysystem for this graticule line is used
to format position labels in hour, minutes,
seconds or degrees, minutes and seconds.
addx: If not set to None, it implies that we supplied
this method with known positions in world
coordinates or pixels. If world coordinates
the coordinate transformations from
pixel- to world coordinates are not necessary.
The parameters addx and addy can be used if you have
coordinates that define a border of an
all-sky map, but to avoid drawing the border
outside the limits in pixels, we treat the line
as a graticule line so it is clipped correctly.
Note that if you want to use addx and addy then
the value of wcsaxis should not be 0 or 1.
addy: The same as addx but now for y
addpixels: True or False.
The values in addx and addy are either pixel
coordinates or world coordinates.
offsetlabel:Text label for an offset axis
fun: Change the offset values using a function.
fmt: Format the offset using a Python format.
Returns: A graticule line which consists of one or more
line pieces with default (plot) attributes,
a sequence of ticks (each with a position, label
and default plot attributes)
Notes: Special graticule lines are those which sets
the limb of the graticule. Not all projections
behave in a way that limbs can be found with
boundaries in world coordinates. Sometimes
we have a boundary in world coordinates but
it turns out to be not as accurate as possible.
For those situations we have a rather crude
method to find boundaries in pixels. To avoid
jumps and to apply clipping at the edges we
process these pixels like other graticule lines.
"""
#-----------------------------------------------------------------
# Helper functions
def __inbox(x, y, box):
#--------------------------------------------------------------
"""
Purpose: Is a position (x,y) in pixels within the
limits of the image?
Parameters:
x: X-coordinate of pixel
y: Y-coordinate of pixel
box: Tuple of 4 numbers (xlo, ylo, xhi, yhi)
Returns: False or True
Notes: Note that internally the box is increased with
1/2 pixel in all directions
"""
#--------------------------------------------------------------
if box[0] <= box[2]:
if x >= box[2] or x <= box[0]:
return False
else:
if x >= box[0] or x <= box[2]:
return False
if box[1] <= box[3]:
if y >= box[3] or y <= box[1]:
return False
else:
if y >= box[1] or y <= box[3]:
return False
return True
def __handlecrossing(box, x1, y1, x2, y2):
#-----------------------------------------------------------
"""
Purpose: Return properties of intersection of
graticule line and enclosing rectangle
Parameters:
box: As in helper function __inbox
x1, y1: First position in pixels inside or outside box
x2, y2: Second position outside or inside box
Returns: The axis number and the position of the intersection
Notes: Given the boundaries of a rectangle in grids,
and two positions of which we know that one
is inside the box and the other is outside the
box, calculate the intersections of the line
through these points and all of the axes of
the box. Note that we used the following
axis-index relation.
0: left
1: bottom
2: right
3: top
"""
#-----------------------------------------------------------
for axisnr in range(4):
if axisnr == left or axisnr == right:
# This is the left and right Y-axis for which x2 != x1
# otherwise we would not have an intersection
if x2-x1 != 0.0:
lamb = (box[axisnr]-x1)/(x2-x1)
else:
lamb = -1.0;
if 0.0 <= lamb <= 1.0:
ycross = y1 + lamb*(y2-y1)
xlab = box[axisnr]; ylab = ycross
# pylab.plot([box[axisnr]], [ycross], 'ro')
return axisnr, xlab, ylab
else:
# Check intersection with y axis
if y2-y1 != 0.0:
lamb = (box[axisnr]-y1)/(y2-y1)
else:
lamb = -1.0
if 0.0 <= lamb <= 1.0:
xcross = x1 + lamb*(x2-x1)
ylab = box[axisnr]; xlab = xcross
# pylab.plot([xcross], [box[axisnr]], 'yo')
return axisnr, xlab, ylab
#-------------------------------------------------
# Start definition of __init__()
#-------------------------------------------------
if wcsaxis == 0: # constant value on X axis
xw = numpy.zeros(linesamples) + constval
dw = (wylims[1] - wylims[0])/1000.0 # Ensure that we cross borders
yw = numpy.linspace(wylims[0]-dw, wylims[1]+dw, linesamples)
elif wcsaxis == 1: # constant value on Y axis
dw = (wxlims[1] - wxlims[0])/1000.0
xw = numpy.linspace(wxlims[0]-dw, wxlims[1]+dw, linesamples)
yw = numpy.zeros(linesamples) + constval
else: # A grid line without one of the coordinates being constant e.g. a border
xw = addx
yw = addy
if (mixgrid is None): # Probably matching axis pair or two independent axes
world = (xw, yw)
if wcsaxis in [0,1]:
pixel = gmap.topixel(world)
else: # This should be an 'added' grid line
if addpixels: # A set with pixels for a border is already provided
pixel = (addx, addy)
else: # A set with world coordinates for a border is provided.
pixel = gmap.topixel(world)
else: # wcs coordinate pair with one extra grid coordinate
unknown = numpy.zeros(linesamples)
unknown += numpy.nan
zp = numpy.zeros(linesamples) + mixgrid
world = (xw, yw, unknown)
pixel = (unknown, unknown, zp)
(world, pixel) = gmap.mixed(world, pixel, iter=10)
#print "WCSGRAT world=",world
#print "WCSGRAT pixel=",pixel
if wcsaxis == 0 or wcsaxis == 1:
self.axtype = gmap.types[wcsaxis]
else:
# E.g. grid lines made without a constant value do not belong to an axis
# so we use a dummy type 'border'
self.axtype = 'border'
self.skysys = skysys
self.wcsaxis = wcsaxis
if wcsaxis in [0,1]:
if gmap.types[wcsaxis] == 'longitude':
# This must be a longitude axis. We want the range to be limited to [0,360>
while constval < 0.0:
constval += 360.0
while constval >= 360.0:
constval -= 360.0
self.constval = constval
self.linepieces = []
# For each plot axis we store the ticks that belong to that axis.
# A tick also belongs to a wcs axis (e.g. a longitude or latitude)
# but this is then an attribute of the tick
self.ticks = []
self.kwargs = {}
# Special care for 'jumps' in a plot. A jump can occur for example near 180 deg where
# for an all sky plot 180 deg is a position left in the box while a world
# coordinate > 180 deg is a position at the right side of the box.
# When stepping near these (plot) discontinuities, we want to avoid to connect
# positions that are more than factor*step size separated from each other.
lastx = lasty = None # Flag for reset of counters is set after a 'jump'
countin = 0; countout = 0
lastinside = False
stepy = (pylims[1] - pylims[0] + 1.0)/ linesamples
stepx = (pxlims[1] - pxlims[0] + 1.0)/ linesamples
box = (pxlims[0]-0.5, pylims[0]-0.5, pxlims[1]+0.5, pylims[1]+0.5)
for p in zip(pixel[0], pixel[1]): # Works for 2 and 3 element world tuples
xp = p[0]; yp = p[1]
# print "WCSGRAT xp, yp", xp, yp
#if wcsaxis == 0 and numpy.isnan(xp):
# print xp, yp, box, constval
if not numpy.isnan(xp) and not numpy.isnan(yp): # NaN's can occur with divergent projections
currentinside = __inbox(xp, yp, box)
if lastx is not None and (lastinside or currentinside):
# These jump conditions are somewhat arbitrary.
# If the projection diverges then these steps sizes may be
# not enough to have a full coverage.
jump = abs(lastx-xp) > abs(10.0*stepx) or abs(lasty-yp) > abs(10.0*stepy)
if jump:
# print "JUMP: ", lastx, lasty, xp, yp, abs(lastx-xp), abs(lasty-yp), stepx, stepy
if countin > 0:
self.linepieces.append( (x,y) ) # Close current line piece
countout = 0; countin = 0
lastx = lasty = None;
if countin == 0:
x = []; y = []
crossing2in = crossing2out = False
if currentinside:
if countout > 0:
# We are going from outside to inside
crossing2in = True
countout = 0
else:
x.append(xp); y.append(yp)
countin += 1
else:
if countin > 0:
# We are crossing from inside to outside
crossing2out = True
countin = 0
countout += 1
if crossing2in or crossing2out:
axisnr, xlab, ylab = __handlecrossing(box, lastx, lasty, xp, yp)
if self.axtype != 'border': # Border lines sometimes do not have a constant world coordinate
if offsetlabel is not None:
labelvalue = offsetlabel
offs = True
else:
labelvalue = constval
offs = False
tick = WCStick(xlab, ylab, axisnr, labelvalue, wcsaxis, offs, fun=fun, fmt=fmt)
# Set a new default for the font size for the tick labels
if wcsaxis == 0:
tick.kwargs.update({'fontsize':'10'})
else:
tick.kwargs.update({'fontsize':'10'})
self.ticks.append(tick)
x.append(xlab); y.append(ylab) # Add this intersection as element of the graticule line
if crossing2out:
self.linepieces.append( (x,y) )
lastx = xp; lasty = yp
lastinside = currentinside
if countin > 0: # Store what is left over
self.linepieces.append( (x,y) )
class PLOTaxis(object):
#-------------------------------------------------------------------------------
"""
Each (plot) axis can have different properties related to the ticks,
and the labels. Labels can be transformed using an external function
"""
#------------------------------------------------------------------------------
def __init__(self, axisnr, mode, label, **kwargs):
#-----------------------------------------------------------
"""
Purpose: Create object that represents an (plot) axis
and store default attributes. Only its
attributes should be used by a user/programmer.
Parameters:
axisnr - 0: left
1: bottom
2: right
3: top
mode - What should this axis do with the tick
marks and labels?
0: ticks native to axis type only
1: only the tick that is not native to axis type
2: both types of ticks (map could be rotated)
3: no ticks
label - An annotation of the current axis
Returns: Object with attributes 'axisnr', 'mode',
'label' and 'kwargs'
Notes: Each plot axis is associated with a PLOTaxis
instance.
"""
#-----------------------------------------------------------
self.axisnr = axisnr # left=0, bottom=1, right=2, top=3
self.mode = mode # Set which (native/not native) labels should be stored for this axis
self.label = label # Default title for this axis
self.kwargs = kwargs # Keyword aguments for the axis labels
self.xpos = None # Position of label (normalized dev. coords)
self.ypos = None
[docs]class Insidelabels(object):
"""
A small utility class for wcs labels inside a plot with a graticule.
Useful for all sky plots.
.. automethod:: setp_label
"""
class Ilabel(object):
def __init__(self, Xp, Yp, Xw, Yw, labval, rots, axtype,
skysys=None, fun=None, fmt=None,
offset=False, prec=0, tex=True, texsexa=True, **kwargs):
self.x = Xp
self.y = Yp
self.xw = Xw
self.yw = Yw
self.label = ""
self.labval = labval
self.rotation = rots
self.axtype = axtype
self.skysys = skysys
self.offset = offset
self.fmt = fmt
self.fun = fun
self.prec = prec
self.tex = tex
self.texsexa = texsexa
if not tex:
self.texsexa = False
self.kwargs = kwargs
def __init__(self, wcsaxis):
self.labels = []
self.ptype = "Insidelabels"
self.pxlim = None
self.pylim = None
self.wcsaxis = wcsaxis
def append(self, Xp, Yp, Xw, Yw, labval, rots, axtype, skysys=None,
fun=None, fmt=None, offset=False, prec=0, tex=True, texsexa=True,
**kwargs):
#-----------------------------------------------------------------
"""
Append a new 'Ilabel' object to the list
"""
#-----------------------------------------------------------------
if not tex: # Overrule input if LaTeX is not used
texsexa = False
ilab = self.Ilabel(Xp, Yp, Xw, Yw, labval, rots, axtype,
skysys=skysys, fun=fun, fmt=fmt,
offset=offset, prec=prec, tex=tex, texsexa=texsexa,
**kwargs)
self.labels.append(ilab)
[docs] def setp_label(self, position=None,
tol=1e-12, fmt=None, fun=None, tex=None, texsexa=None,
**kwargs):
#-----------------------------------------------------------------
"""
This method handles the
properties of the 'inside' labels, which are :class:`Text` objects
in Matplotlib. The most useful properties are *color*, *fontsize*
and *fontstyle*. One can change the label values using an external
function and/or change the format of the label.
:param position:
Accepted are None, or one or more values representing
the constant value of the graticule line in
world coordinates. These positions are used to identify
individual graticule lines so that each line can have its
own properties. If no position is entered, then the changes
are applied to all the labels in the current object.
The input can also be a string that represents a sexagesimal number.
:type position: *None* or one or a sequence of floating point numbers
:param tol:
If a value > 0 is given, the gridline with the
constant value closest to a given position within
distance 'tol' gets updated attributes.
:type tol: Floating point number
:param fmt:
A string that formats the tick value
e.g. ``fmt="%10.5f"`` in the Python way, or a string
that contains no percentage character (%) but a format
to set the output of sexagesimal numbers e.g.
fmt='HMs'. The characters in the format either force
(uppercase) a field to be printed, or it suppresses
(lowercase) a field to be printed.
See also the examples at :func:`makelabel`.
:type fmt: String
:param fun:
An external function which will be used to
convert the tick value e.g. to convert
velocities from m/s to km/s. See also
example 2_ below.
:type fun: Python function or Lambda expression
:param tex:
If True then format the tick label in LaTeX. This is the
default. If False then standard text will be applied.
Some text properties cannot be changed if LaTeX is
in use.
:type tex: Boolean
:param texsexa:
If False and parameter *tex* is True, then format the
tick label without superscripts for sexagesimal labels.
This option can be used if superscripts result in 'jumpy' labels.
The reason is that in Matplotlib the TeX labels at the bottom
of a plot are aligned at a baseline at the top of the characters and
not at the bottom, while the height between LaTeX boxes may vary.
:type tex: Boolean
:param `**kwargs`:
Keyword arguments for plot properties like *color*,
*visible*, *rotation* etc. The plot attributes are standard
Matplotlib attributes which can be found in the
Matplotlib documentation.
:type `**kwargs`: Matplotlib keyword arguments
:note:
Some projections generate labels that are very close
to each other. If you want to skip labels then you can
use keyword/value *visible=False*. Note that *visible*
is a parameter of Matplotlib's plot functions.
"""
#-----------------------------------------------------------------
if not tex is None:
if not tex:
texsexa = False
if position is not None:
if not issequence(position):
posn = [position]
else:
posn = position
for label in self.labels:
if position is None:
if not fmt is None:
label.fmt = fmt
if not fun is None:
label.fun = fun
if not tex is None:
label.tex = tex
if not texsexa is None:
label.texsexa = texsexa
label.kwargs.update(kwargs)
else: # One or more positions are given. Find the right index.
for pos in posn:
if isinstance(pos, str):
C = pos.upper()
if 'H' in C or 'D' in C or 'M' in C or 'S' in C:
pos, err = parsehmsdms(pos)
if err != '':
raise Exception(err)
else:
pos = pos[0]
else:
raise ValueError("[%s] is entered as a string but does not represent valid HMS or DMS"%pos)
if self.wcsaxis == 0:
d = abs(label.xw - pos)
else:
d = abs(label.yw - pos)
if d <= tol:
if not fmt is None:
label.fmt = fmt
if not fun is None:
label.fun = fun
if not tex is None:
label.tex = tex
if not texsexa is None:
label.texsexa = texsexa
label.kwargs.update(kwargs)
def plot(self, frame):
#-----------------------------------------------------------
"""
Purpose: Plot world coordinate labels inside a plot
Parameters: 'frame', a Matplotlib Axes object
Returns: --
Notes:
"""
#-----------------------------------------------------------
createlabels(self.labels)
# With Insidelabels there is no problem with the alignment
# with TeX labels, because the alignment is not 'top'.
# There is no need to tweak the HMS labels here.
for inlabel in self.labels:
frame.text(inlabel.x, inlabel.y, inlabel.label,
clip_on=True, **inlabel.kwargs)
# Set limits
xlo = self.pxlim[0]-0.5
ylo = self.pylim[0]-0.5
xhi = self.pxlim[1]+0.5
yhi = self.pylim[1]+0.5
frame.set_xlim((xlo,xhi))
frame.set_ylim((ylo,yhi))
def getlambda(pixel, lo, hi):
#-----------------------------------------------------------------------------
"""
Small utility to calculate lambda on a line for given position
in pixels
"""
#-----------------------------------------------------------------------------
if pixel is None:
return 0.5
delta = hi - lo
if delta == 0.0:
return 0.5
return (pixel-lo)/(delta)
def getStarts(startvals, proj, i, axnum, wcstype, mixpix=None):
#------------------------------------------------------------------------------
"""
Helper function for Graticule contructor.
This function converts user supplied positions on an axis to numbers
if the values are given as one string. It extends the functionality
of function str2pos(). Its main purpose is to set one or more values for
axis labels if a user does not want the default values as calculated
in the constructor of class Graticule.
Notes:
The projection object 'proj' represents two axes of the map if
the map is either spatial or does not have any spatial axis.
If the map has only one spatial axis then a third axis is added
to the projection object and the valye of mixpix is set. Note that
function str2pos() diminishes the dimensionality with one is
the value of mixpix is set.
So distinguish three options:
1) This function is called with an axis that can be represented by
a projection object with one axis like a spectral axis or a linear axis.
2) This function is called with a spatial axis and the other axis in the map
is not spatial. The value of mixpix should be set.
3) This function is called with a spatial axis and the other axis in the map
is also spatial. The value of mixpix is not set in this case.
"""
#------------------------------------------------------------------------------
if startvals is None:
# This is the default. Do nothing
return startvals, None, ""
if not isinstance(startvals, str):
# Perhaps the values are one or more floating point numbers
return startvals, None, ""
spatial = wcstype in ['lo', 'la']
if spatial:
# If the is a value for mixpix then the map had only one spatial axis.
# The missing spatial axis is then always the last.
if not mixpix is None:
subaxnum = (i+1, 3)
else:
# But we need a mixpix here and None was given.
if i == 0:
mixpix = proj.crpix[1]
subaxnum = (1,2)
else:
mixpix = proj.crpix[0]
subaxnum = (2,1)
else:
# We need only one. i is either 0 (first axis) or 1 (second axis)
subaxnum = (i+1,) # Make tuple
mixpix = None # Prevent str2pos() to decrease dimensionality !
subproj = proj.sub(subaxnum)
world, pixels, units, errmes = str2pos(startvals, subproj, mixpix=mixpix)
if errmes != '':
return None, None, errmes
#w = world.flatten()
w = world[:,0] # Get rid of extra spatial if there is one
p = pixels[:,0]
if len(w) == 1:
w = w[0] # No longer a sequence but a real start value.
p = p[0]
return w, p, ""
[docs]class Graticule(object):
#-----------------------------------------------------------------
"""
Creates an object that defines a graticule
A (spatial) graticule consists of parallels and meridians. We extend this to
a general grid so we can cover every type of map (e.g. position velocity maps).
:param header: Is a Python dictionary or dictionary-like object
containing FITS-style keys and values, e.g. a
header object from PyFITS.
Python dictionaries are used for debugging,
or plotting experiments or when you need to
define a projection system from scratch.
:type header: Python dictionary or FITS header object (pyfits.NP_pyfits.HDUList)
:param graticuledata: This is a helper object. It can be any object as long it
has attributes:
* header
* axnum
* pxlim
* pylim
* mixpix
* spectrans
Software that interfaces with a user to get data
and relevant properties could/should produce objects which
have at least values for the attributes listed above.
Then these objects could be used as a shortcut parameter.
:type graticuledata: Object with some required attributes
:param axnum: This parameter sets which FITS axis corresponds
to the x-axis of your graticule plot rectangle
and which one corresponds to the y-axis
(see also description at *pxlim* and *pylim*).
The first axis in a FITS file is axis 1.
If *axnum* set to *None* then the default
FITS axes will be 1 and 2.
With a sequence you can set different FITS axes
like ``axnum=(1,3)`` Then the input is a tuple
or a list.
:type axnum: None, Integer or sequence of Integers
:param wcstypes: List with the type of the used axes. These types are
derived from the projection object axis types (attribute
wcstype) but are translated into a string: The strings
are 'lo' for a longitude axis, 'la' for a latitude axis,
'sp; for a spectral axis and 'li_xxx' for a linear axis
where 'xxx' is the ctype for that axis.
:type wcstypes: List of strings
:param pxlim: The values of this parameter together with
the values in pylim define a rectangular frame.
The intersections of graticule lines with this
frame are the positions where want
to plot a tick mark and write a label that
gives the position as a formatted string.
Further, the limits in pixels are used to set
the step size when a graticule line is sampled.
This step size then is used to distinguish
a valid step from a jump (e.g. from 180-delta
degrees to 180+delta degrees which can jump from one side
in the plot to the other side).
To prevent a jump in a plot, the graticule
line is splitted into line pieces without jumps.
The default of *pxlim* is copied from the header
value. FITS
data starts to address the pixels with 1 and the last pixel
is given by FITS keyword *NAXISn*.
Note that internally the enclosing rectangle
in pixels is enlarged with 0.5 pixel in all
directions. This enables a correct overlay on an
image where the pixels have a size.
:type pxlim: *None* or exactly 2 Integers
:param pylim: See description at pxlim. The range is along the
y-axis.
:type pylim: *None* or exactly 2 Integers
:param mixpix: For maps with only 1 spatial coordinate we need to
define the pixel that sets the spatial value
on the matching spatial axis. If its value is
*None* then the value of *CRPIXn* of the matching
axis from the header is taken as default.
:type mixpix: *None* or 1 Integer
:param spectrans: The spectral translation. For spectral axes
it is usually possible to convert to another
representation. For instance one can 'translate'
a frequency into a velocity which is one of
the types: VOPT-F2W, VRAD, VELO-F2V
(for optical, radio and radial velocities).
See also the article
`Representations of spectral coordinates in FITS <http://www.atnf.csiro.au/people/mcalabre/WCS/scs.pdf>`_
by Greisen, Calabretta, Valdes & Allen.
Module *maputils* from the Kapteyn Package provides
a method that creates a list with possible spectral
translations given an arbitrary header.
The spectral translation should be followed by
a code (e.g. as in 'VOPT-F2W') which sets the
conversion algorithm. If you
don't know this beforehand, you can either append
the string '-???' or try your translation without this
coding. Then this module tries to find the appropriate
code itself.
:type spectrans: String
:param skyout: A single number or a tuple which specifies
the celestial system.
The tuple is laid out as follows:
``(sky system, equinox, reference system,
epoch of observation)``.
Predefined are the systems:
* wcs.equatorial
* wcs.ecliptic,
* wcs.galactic
* wcs.supergalactic
or the minimal matched string versions of these values.
Predefined reference systems are:
* wcs.fk4,
* wcs.fk4_no_e,
* wcs.fk5,
* wcs.icrs,
* wcs.j2000
or the minimal matched string versions of these values.
Prefixes for epoch data are:
============= =================== ======================================
Prefix Description Example
============= =================== ======================================
B Besselian epoch 'B 1950', 'b1950', 'B1983.5', '-B1100
J Julian epoch 'j2000.7', 'J 2000', '-j100.0'
JD Julian Date 'JD2450123.7'
MJD Modified Julian Day 'mJD 24034', 'MJD50123.2'
RJD Reduced Julian Day 'rJD50123.2', 'Rjd 23433'
F DD/MM/YY (old FITS) 'F29/11/57'
F YYYY-MM-DD 'F2000-01-01'
F YYYY-MM-DDTHH:MM:SS 'F2002-04-04T09:42:42.1'
============= =================== ======================================
See the documentation of module :mod:`celestial` for more details.
Example of a sky definition::
skyout = (wcs.equatorial, wcs.fk4_no_e, 'B1950')
:type skyout: *None*, one Integer or a tuple with a sky definition
:param alter: A character from 'A' through 'Z', indicating
an alternative WCS axis description from a FITS header.
:type alter: Character
:param wxlim: Two numbers in units of the x-axis. For spatial
axes this is usually in degrees. The numbers
are the limits of an interval for which
graticules will be calculated. If these values
are omitted, defaults will be calculated.
Then random positions in pixels are converted to
world coordinates and the greatest gap in
these coordinates is calculated. The end- and
start point of the gap are the start- and end point
of the range(s) in world coordinates. It is not enough
to transform only the limits in pixels because a maximum
or minimum in world coordinates could be located
on arbitrary pixel positions depending on the projection.
:type wxlim: *None* or exactly two floating point numbers
:param wylim: See wxlim, but now applied for the y-axis
:type wylim: *None* or exactly two floating point numbers
:param boxsamples: Number of random pixel positions within a box
with limits *pxlim* and *pylim* for which world
coordinates are calculated to get an estimate of
the range in world coordinates (see description
at wxlim). The default is listed in the argument list
of this method. If speed is essential one can try smaller
numbers than the default.
:type boxsamples: Integer
:param startx: If one value is given then this is the
first graticule line that has a constant
x **world coordinate** equal to *startx*.
The other values will be
calculated, either with distance *deltax*
between them or
with a default distance calculated by this method.
If *None* is set, then a suitable value will be
calculated.
The input can also be a string which is parsed by
the positions module. This enables the use of
units etc.
Examples (see also module :mod:`positions`:
* For a frequency axis: startx="linspace(1.4240,1.4250,4) Ghz"
* For a frequency axis: startx="arange(1.4240,1.4250,0.0005) Ghz"
* For a spectral translation to WAVE: startx="'0.2105, 0.2104' m"
* Two labels on a longitude axis: startx="3h00m20s 3h00m30s"
:type startx: *None* or 1 floating point number or a sequence
of floating point numbers or a string.
:param starty: [None, one value, sequence]
Same for the graticule line with constant
y world coordinate equal to starty.
:type starty: *None* or 1 floating point number or a sequence
of floating point numbers or a string.
:param deltax: Step in **world coordinates** along the x-axis
between two subsequent graticule lines.
It can also be a string with an expression and optionally
a unit. Note that the expression cannot contain any spaces.
Example:
* deltax = 5*6/6 dmsmin
:type deltax: *None* or a floating point number or a string
:param deltay: Same as deltax but now as step in y direction.
It can also be a string with an expression and optionally
a unit.
:type deltay: *None* or a floating point number or a string.
:param skipx: Do not calculate the graticule lines with the
constant world coordinate that is associated
with the x-axis.
:type skipx: Boolean
:param skipy: The same as skipx but now associated with
the y-axis.
:type skipy: Boolean
:param gridsamples: Number of positions on a graticule line for which
a pixel position is calculated and stored as part
of the graticule line. If *None* is set then the
default is used (see the argument list of this method).
:type gridsamples: Integer
:param labelsintex: The default is that all tick labels are formatted for LaTeX.
These are not the axes labels. If you want to format these in
LaTeX then you need to set them explicitly as in:
>>> grat.setp_axislabel("bottom",
label=r"$\mathrm{Right\ Ascension\ (2000)}$",
fontsize=14)``
Printing your axis labels in LaTeX
limits the number of Matplotlib properties that one
can set.
:type labelsintex: Boolean
:param offsetx: Change the default mode which sets either plotting
the labels for the given -or calculated world coordinates
or plotting labels which represent constant offsets
with respect to a given starting point.
The offset mode is default for plots with mixed axes,
i.e. with only one spatial axis. In spatial maps
this offset mode is not
very useful to plot the graticule lines because these lines
are plotted at a constant world coordinate and do not know
about offsets.
The offset axes correspond to the pixel positions of
start- and endpoint of
the left and bottom axes and the default start point of
the offsets (value 0) is at the centre of the axis.
One can change this start point with *startx*, *starty*.
:type offsetx: *None* or Boolean
:param offsety: Same as *offsetx* but now for the left plot axis.
:type offsety: *None* or Boolean
:param unitsx: Units for first axis. Applies both to regular and offset axes.
If this parameter sets a unit other than the default,
then a conversion function will be used to display the
labels in the new units. The unit in the default axis label
will be replaced by the new units.
:type unitsx: String
:param unitsy: Units for second axis.
:type unitsy: String
:raises:
:exc:`ValueError` *Could not find enough (>1) valid world coordinates in this map!*
User wanted to let the constructor estimate what the ranges in
world coordinates are for this header, but only zero or one
coordinate could be found.
:exc:`ValueError` *Need data with at least two axes*
The header describes zero or one axes. For a graticule
plot we need at least two axes.
:exc:`ValueError` *Need two axis numbers to create a graticule*
The *axnum* parameter needs exactly two values.
:exc:`ValueError` *Need two different axis numbers*
A user/programmer entered two identical axis numbers.
Graticules need two different axes.
:exc:`ValueError` *pxlim needs to be of type tuple or list*
Check type.
:exc:`ValueError` *pxlim must have two elements*
Number must be exactly 2.
:exc:`ValueError` *pylim needs to be of type tuple or list*
Check type.
:exc:`ValueError` *pylim must have two elements*
Number must be exactly 2.
:exc:`ValueError` *Could not find a grid for the missing spatial axis*
The specification in *axnum* corresponds to a map with only one
spatial axis. If parameter *mixpix* is omitted then the constructor
tries to find a suitable value from the (FITS) header. It
reads *CRPIXn* where n is the appropriate axis number. If nothing
could be found in the header then this exception will be raised.
:exc:`ValueError` *Could not find a matching spatial axis pair*
The specification in *axnum* corresponds to a map with only one
spatial axis. A We need the missing spatial axis to find a
matching world coordinate, but a matching axis could not be found
in the header.
:exc:`ValueError` *wxlim needs to be of type tuple or list*
Check type.
:exc:`ValueError` *wxlim must have two elements*
Number must be exactly 2.
:exc:`ValueError` *wylim needs to be of type tuple or list*
Check type.
:exc:`ValueError` *wylim must have two elements*
Number must be exactly 2.
:exc:`ValueError` *boxsamples < 2: Need at least two samples to find limits*
There is a minimum number of random positions we have to
calculate to get an impression of the axis limits in world coordinates.
:exc:`ValueError` *Number of samples along graticule line must be >= 2 to avoid a step size of zero*
The value of parameter *gridsamples* is too low. Low values give
distorted graticule lines. Higher values (like the default) give
smooth results.
:Returns: A graticule object. This object contains the line
pieces needed to draw the graticule and the
ticks (positions, text and axis number).
The basis method to reveal this data (necessary
if you want to make a plot yourself) is described in the
following example::
graticule = wcsgrat.Graticule(header)
for gridline in graticule:
print("\\nThis gridline belongs to axis", gridline.wcsaxis)
print("Axis type: %s. Sky system %s:" % (gridline.axtype, gridline.skysys))
for t in gridline.ticks:
print("tick x,y:", t.x, t.y)
print("tick label:", t.labval)
print("tick on axis:", t.axisnr)
for line in gridline.linepieces:
print("line piece has %d elements" % len(line[0]))
.. Note::
A Graticule object has a string representation and can
therefore be easily inspected with Python's **print** statement.
**Attributes:**
.. attribute:: axes
Read the PLOTaxis class documentation.
Four PLOTaxis instances, one for each axis of the
rectangular frame in pixels set by *xplim* and *pylim*
If your graticule object is called **grat** then
the four axes are accessed with:
* grat.axes[wcsgrat.left]
* grat.axes[wcsgrat.bottom]
* grat.axes[wcsgrat.right]
* grat.axes[wcsgrat.top]
Usually these attributes are set with method :meth:`setp_plotaxis()`.
Examples:
::
grat.axes[wcsgrat.left].mode = 1
grat.axes[wcsgrat.bottom].label = 'Longitude / Latitude'
grat.axes[wcsgrat.bottom].mode = 2
grat.axes[wcsgrat.right].mode = 0
::
PLOTaxis modes are:
0: ticks native to axis type only
1: Only the tick that is not native to axis type
2: both types of ticks (map could be rotated)
3: no ticks
The default values depend on how many ticks, native
to the plot axis, are found. If this is < 2 then
we allow both native and not native ticks along
all plot axes.
.. attribute:: pxlim
The limits of the map in pixels along the x-axis.
This value is either set in the constructor or
calculated. The default is *[1,NAXISn]*.
The attribute is meant as a read-only attribute.
.. attribute:: pylim:
Same for the y-axis.
.. attribute:: wxlim
The limits of the map in world coordinates for the
x-axis either set in the constructor or calculated
(i.e. estimated) by this method. The attribute is
meant as a read-only attribute.
.. attribute:: wylim
Same for the y-axis
.. attribute:: xaxnum
The (FITS) axis number associated with the x-axis
Note that axis numbers in FITS start with 1. If these
numbers are not given as argument for the
constructor then *xaxnum=1* is assumed.
The attribute is
meant as a read-only attribute.
.. attribute:: yaxnum
Same for the y-axis.
Default: *yaxnum=2*
.. attribute:: wcstypes
List with strings that represent the wcs axis type
of the axes.
.. attribute:: gmap
The wcs projection object for this graticule.
See the *wcs* module document for more information.
.. attribute:: mixpix
The pixel on the matching spatial axis for maps
with only one spatial axis. This attribute is
meant as a read-only attribute.
.. attribute:: xstarts
World coordinates associated with the x-axis
which set the constant value of a graticule line
as calculated when the object is initialized.
This attribute is meant as a read-only attribute.
.. attribute:: ystarts
Same for the y-axis
.. attribute:: skyout
Unformatted copy of input parameter *skyout*
.. attribute:: spectrans
Unformatted copy of input parameter *spectrans*
:Examples: Example to show how to use a custom made header to
create a graticule object. Usually one uses this option
to create **all sky** plots. It is also a useful tool
for experiments.::
#1. A minimal header for an all sky plot
header = {'NAXIS' : 2, 'NAXIS1': 100, 'NAXIS2': 80,
'CTYPE1' : 'RA---AZP', 'CRVAL1' :0,
'CRPIX1' : 50, 'CUNIT1' : 'deg', 'CDELT1' : -5.0,
'CTYPE2' : 'DEC--AZP',
'CRVAL2' : dec0, 'CRPIX2' : 40, 'CUNIT2' : 'deg',
'CDELT2' : 5.0,
'PV2_1' : mu, 'PV2_2' : gamma,
}
grat = wcsgrat.Graticule(header)
Use module `PyFITS <http://www.stsci.edu/resources/software_hardware/pyfits>`_
to read a header from a FITS file::
#2. A header from a FITS file 'test.fits'
import pyfits
hdulist = pyfits.open('test.fits')
header = hdulist[0].header
grat = wcsgrat.Graticule(header)
Select the axes for the graticules. Note that the order
of the axes should be the same as the order of axes in
the image where you want to plot the graticule. If necessary
one can swap the graticule plot axes with input parameter
*axnum*::
#3. Swap x and y- axis in a FITS file
grat = wcsgrat.Graticule(header, axnum= (2,1))
For data with more than two axes, one can select the axes
with input parameter *axnum*::
#4. For a FITS file with axes (RA,DEC,FREQ)
# create a graticule for the FREQ,RA axes:
grat = wcsgrat.Graticule(header, axnum=(3,1))
Use sexagesimal numbers for *startx*/*starty*::
#5. Sexagesimal input
grat = wcsgrat.Graticule(...., startx="7h59m30s", starty="-10d0m30s')
**Methods which set (plot) attributes:**
.. automethod:: setp_tick
.. automethod:: setp_plotaxis
.. automethod:: setp_lineswcs0
.. automethod:: setp_lineswcs1
.. automethod:: setp_gratline
.. automethod:: setp_axislabel
.. automethod:: setp_tickmark
.. automethod:: setp_ticklabel
.. automethod:: set_tickmode
**Methods that deal with special curves like borders:**
.. automethod:: scanborder
.. automethod:: addgratline
.. automethod:: setp_linespecial
**Methods related to plotting derived elements:**
.. automethod:: Insidelabels
.. automethod:: Insidelabels.setp_label
**Utility methods:**
.. automethod:: get_aspectratio
"""
#-----------------------------------------------------------------
@staticmethod
def __bisect(direct, const, var1, var2, gmap, tol):
#-----------------------------------------------------------
"""
Purpose: Apply bisection to find the position of
a limb (border in a plot).
Parameters:
direct: direct=0: Apply bisection in y-direction
direct=1: Apply bisection in x-direction
const: Pixel position of the axis along which
a bisection applied
var1, var2: Lower and upper limits in pixels of the
interval along which bisection is applied.
gmap: The wcs projection object (to apply methods
toworld and topixel
tol: The tolerance in pixels used as a stop
criterion for the bisection.
Returns: The pixel position in range [var1, var2]
which represent a position on the border.
Notes: A limb in terms of plotting is a border
which defines the regions where conversions
from and to world coordinates is possible.
Those positions where this is not possible
are represented by NaN's. So we try to find
the position where there is a transition
between a number and a NaN within a certain
precision given by parameter 'tol'.
The bisection is applied either in x- or
y-direction at value 'const'
and the start interval is between var1 and var2.
"""
#-----------------------------------------------------------
Nmax = 100
if direct == 0: # Vertical bisection
xw1, yw1 = gmap.toworld((const, var1))
xw2, yw2 = gmap.toworld((const, var2))
vw1 = yw1; vw2 = yw2
else:
xw1, yw1 = gmap.toworld((var1, const))
xw2, yw2 = gmap.toworld((var2, const))
vw1 = xw1; vw2 = xw2
# One position must be a number and the other must be a NaN
if (not (numpy.isnan(vw1) or numpy.isnan(vw2))) or (numpy.isnan(vw1) and numpy.isnan(vw2)):
return None;
if numpy.isnan(vw1):
vs = var1
ve = var2
else:
vs = var2
ve = var1
i = 0
while i <= Nmax:
vm = (vs + ve)/2.0
if direct == 0:
xw, yw = gmap.toworld((const, vm))
v = yw
else:
xw, yw = gmap.toworld((vm, const))
v = xw
if numpy.isnan(v):
vs = vm
else:
ve = vm
if abs(ve-vs) <= tol:
break
i += 1
return vs
# @staticmethod
def sortticks(self, tex=False):
#----------------------------------------------------------
"""
Purpose: Collect ticks for each plot axis
Format the labels if appropriate.
Parameters:
tex: True or False. If True then render the labels
in TeX.
Returns: tickpos, ticklab, tickkwa which are all 4 lists
of with tick information per plot axis.
Notes: The ticks are sorted per axis because then one
can set properties for all ticks that belong
to one axis.
"""
#-----------------------------------------------------------
global tweakhms
ticks = [[],[],[],[]]
tickpos = [[],[],[],[]]
ticklab = [[],[],[],[]]
tickkwa = [[],[],[],[]]
tickMkwa = [[],[],[],[]]
for gridline in self.graticule:
wcsaxis = gridline.wcsaxis
for t in gridline.ticks:
anr = t.axisnr
mode = self.axes[anr].mode
skip = False
if mode == 0: # Include a tick. Select ticks along axis with axis 'mode'.
# Plot only tick labels near the relevant axes
if wcsaxis == 0 and (anr == left or anr == right):
skip = True
if wcsaxis == 1 and (anr == bottom or anr == top):
skip = True
elif mode == 1:
# Plot only the 'not native' tick labels near the relevant axes
if wcsaxis == 1 and (anr == left or anr == right):
skip = True
if wcsaxis == 0 and (anr == bottom or anr == top):
skip = True
# Mode == 2 allows all ticks for this axis
elif mode == 3:
# No ticks at all option
skip = True
if not skip:
if anr in [left,right]:
tickpos[anr].append(t.y)
else:
tickpos[anr].append(t.x)
tickkwa[anr].append(t.kwargs)
tickMkwa[anr].append(t.markkwargs)
#ticklab[anr].append(t.label)
# Give this tick some extra attributes necessary for formatting its label
t.axtype = gridline.axtype
t.skysys = gridline.skysys
t.prec = self.prec[wcsaxis]
t.delta = self.delta[wcsaxis]
# Tex mode could have been changed by a set properties method
if t.tex is None:
t.tex = tex
ticks[anr].append(t)
# We reached the stage that we assembled all ticks that belong
# to an axis. This implies that the ticks can belong to different
# wcs axes. For the formatting of the labels one needs to distinguish
# those.
for anr in range(0,4):
for wcsaxis in [0,1]:
Tlist = []
for t in ticks[anr]:
if t.wcsaxis == wcsaxis:
Tlist.append(t)
# Modify the labels. Note that the ticklab list points to the
# label attribute of a WCStick object. The createlabels() function
# returns a list so in fact if this function changes a label for
# a tick object, then
# the label will automatically be updated in the list also.
createlabels(Tlist)
#if anr == bottom:
# Only bottom axes can have TeX labels that are not aligned very well
# With this tweak parameter this should improve
#update_metrics(tweakhms)
for anr in range(0,4):
for t in ticks[anr]:
ticklab[anr].append(t.label)
return tickpos, ticklab, tickkwa, tickMkwa
@staticmethod
def __adjustlonrange(lons):
#-----------------------------------------------------------
"""
Purpose: Find minimum and maximum of range in longitudes
with gaps. E.g. with a crval near zero, one expects
values both negative as positive. However the wcs
routines return longitudes always in range [0,360>.
So -10 appears as 250 in the list. This results
in the wrong min and max of this range. This
algorithm calculates the two array values for which
the gap is the biggest. It then returns the correct
min and max of the range with the min always
as first parameter (allowing for negative values).
Parameters:
lons A (numpy) 1-dim array of longitudes
(world coordinates)
Returns: min, max of range of input longitudes, excluding
the biggest gap smaller than or equal to 180
(degrees) in the range.
Examples: 1) Longitudes: [270, 220, 88, 12, 90, 0, 289, 180, 300, 2, 3, 4]
Sorted longitudes: [0, 2, 3, 4, 12, 88, 90, 180, 220, 270, 289, 300]
Biggest gap, start longitude, end longitude 90 -180.0 90
min, max: -180.0 90
2) Longitudes: [1, 3, 355, 2, 5, 7, 0, 359, 350, 10, 11, 349]
Sorted longitudes: [0, 1, 2, 3, 5, 7, 10, 11, 349, 350, 355, 359]
Biggest gap, start longitude, end longitude 22.0 -11.0 11
min, max: -11.0 11
"""
#------------------------------------------------------------
def __diff_angle(a, b):
if b < a:
result = b + 360 - a
else:
result = b - a
if result > 180.0:
result -= 360.0
return result
# Eliminate largest gap
lons.sort()
gap = 0.0
prev = lons[-1]
brkpt = 0
for i, lon in enumerate(lons):
diff = __diff_angle(prev, lon)
if abs(diff)>gap:
gap = abs(diff)
brkpt = i
prev = lon
# 'gap' contains now the largest gap
lon1 = lons[brkpt]; lon2 = lons[brkpt-1]
if lon1 > lon2:
lon1 -= 360.0
return lon1, lon2
@staticmethod
def __nicenumbers(x1, x2, start=None, delta=None, axtype=None, skysys=None, offset=False):
#-----------------------------------------------------------
"""
Purpose: Find suitable numbers to define graticule lines
in interval [x1,x2]. Process longitudes and
latitudes in seconds.
Also return a list with the same length which
contains offsets only.
Parameters:
x1, x2: A start and an end value representing an interval
start: If not None, then include this value in the list
of nice numbers.
delta: If not None, use this value as step size
If both 'start' and 'delta' are not None, then use
these values to get all the values in the given
interval [x1,x2] with start equal to 'start' and
step size equal to 'delta'
axtype: None or one of ('longitude', latitude, 'spectral').
Value is used to distinguish lons and lats from
other data.
skysys: For longitudes distinguish ranges 0,360 in hours
(equatorial system) or degrees.
Returns: A tuple with:
-A NumPy array with 'nice' numbers
-The precision in seconds
-The proposed step size
Notes: Spatial intervals are first multiplied by the
appropriate factor to get an interval in seconds.
For numbers >= 10 seconds, a 'nice' step size is
drawn from a predefined list with numbers.
For other intervals, the length of the interval
is scaled to an interval between [0,10>. The
scaled length sets the number of divisions and
the final step size. The output step size is in degrees.
Examples: a=3; b = 11
print(nicenumbers(a, b))
print(nicenumbers(a, b, start=10))
print(nicenumbers(a, b, start=10, delta=1))
print(nicenumbers(a, b, delta=2))
(array([ 6., 8., 10.]), 0, 2.0)
(array([ 10., 8., 6., 4.]), 0, 2.0)
(array([ 10., 9., 8., 7., 6., 5., 4.]), 0, 1.0)
(array([ 4., 6., 8.]), 0, 2.0)
Other examples:
nicenumbers(a, b, axtype='longitude', skysys=wcs.equatorial)
nicenumbers(a, b, axtype='longitude', skysys=wcs.ecliptic)
nicenumbers(a, b, axtype='latitude', skysys=wcs.galactic)
"""
#-----------------------------------------------------------
prec = 0
step = None
#nc = None
fact = 1.0
if x2 < x1: # Then swap because we want x2 > x1
x1,x2 = x2,x1
x1orig = x1; x2orig = x2
dedge = (x2-x1) / 80.0 # 80 is a 'magic' number that prevents
# graticule lines too close to the borders
if not issequence(skysys):
skysys = [skysys]
if delta is None:
if axtype in ['longitude', 'latitude']:
# Nice numbers for dms should also be nice numbers for hms
sec = numpy.array([30, 20, 15, 10, 5, 2, 1])
minut = sec
deg = numpy.array([60, 30, 20, 15, 10, 5, 2, 1])
nicenumber = numpy.concatenate((deg*3600.0, minut*60.0, sec))
if wcs.equatorial in skysys and axtype == 'longitude':
fact = 240.0
else:
fact = 3600.0
x1 *= fact; x2 *= fact
d = x2 - x1
step2 = 0.9*d/3.0 # We want at least two graticule lines in this range
for p in nicenumber:
k = int(step2/p)
if k >= 1.0:
step2 = k * p
step = step2
# nc = nicenumber
break # Stop if we have a candidate
if step is None:
d = x2 - x1
f = int(numpy.log10(d))
if d < 1.0:
f -= 1
D3 = numpy.round(d/(10.0**f),0)
if D3 == 3.0:
D3 = 2.0
elif D3 == 6:
D3 = 5.0
elif D3 == 7:
D3 = 8
elif D3 == 9:
D3 = 10
if D3 in [2,4,8]:
k = 4
else:
k = 5
step = (D3*10.0**f)/k
else:
step = delta
if step == 0.0: # Return empty list. Cannot divide by zero
return [], 0,0 # To get just one point, use delta > x2 - x1
xxm = step*numpy.round(((x1+x2)/2.0)/step)
xxm /= fact; step /= fact
# Now both parameters are in the original units.
# Calculate a precision for the seconds along spatial axes.
pstep = step
if axtype in ['longitude', 'latitude']:
if wcs.equatorial in skysys and axtype == 'longitude':
pstep *= 240.0
else:
pstep *= 3600.0
# Set the precision for the seconds
# We found an optimum precision by trial and error
f0 = numpy.floor( numpy.log10(pstep) )
if f0 < 0.0:
prec = int(abs(f0))
else:
prec = 0
startx = None
if start is not None:
startx = start
elif x1orig+dedge < 0.0 < x2orig-dedge:
startx = 0.0
else:
startx = xxm
l1 = numpy.arange(startx, x1orig+0.9*dedge, -step)
n1 = len(l1)
o1 = numpy.arange(0.0, -(n1-0.01)*step, -step)
l2 = numpy.arange(startx+step, x2orig-0.9*dedge, step)
n2 = len(l2)
o2 = numpy.arange(0.0+step, (n2+0.01)*step, step)
nice = numpy.concatenate( (l1,l2) )
offsets = numpy.concatenate( (o1,o2) )
return nice, offsets, prec, step
def __estimateLonLatRanges(self, nrandomsamples):
#----------------------------------------------------------
"""
Purpose: Given the current pixel limits, find the
limits along the x and y directions in
world coordinates.
Parameters:
nrandomsamples: Number of random pixel position samples
Returns: wxmin, wxmax, wymin, wymax
The limits in world coordinates
Notes: Given the current ranges in pixel coordinates
(box), estimate the ranges in world coordinates.
The complication is that we are dealing with
many different WCS coordinate transformations.
Therefore we sample the grid with
'nrandomsamples' random positions for which
we assume that the converted longitudes and
latitudes can be used to estimate a close
indication for these ranges. The edges of the
box are also included in the calculations.
If a projection behaves in a way that the edges
are also the limits in world coordinates then
we automatically get the maximum limits,
otherwise it is an approximation and the quality
of this approximation increases if the number
of samples increases.
"""
#-----------------------------------------------------------
xlo = self.pxlim[0]-0.5; ylo = self.pylim[0]-0.5;
xhi = self.pxlim[1]+0.5; yhi = self.pylim[1]+0.5
# Dx = (xhi - xlo + 1)/10.0; Dy = (yhi - ylo + 1)/10.0
Dx = Dy = 0.0
xr = numpy.random.uniform(xlo-Dx, xhi+Dx, nrandomsamples+4)
yr = numpy.random.uniform(ylo-Dy, yhi+Dy, nrandomsamples+4)
# Always include the edges of the frame
xr[0] = xlo; xr[1] = xhi; xr[2] = xlo; xr[3] = xhi;
yr[0] = ylo; yr[1] = ylo; yr[2] = yhi; yr[3] = yhi;
if self.mixpix is None:
pixels = (xr, yr)
else:
zr = numpy.zeros(nrandomsamples+4) + self.mixpix
pixels = (xr, yr, zr)
world = self.gmap.toworld(pixels)
# The world coordinates we want are always the first and second
# element of the result tuple.
wx = numpy.ma.masked_where(numpy.isnan(world[0]), world[0])
wy = numpy.ma.masked_where(numpy.isnan(world[1]), world[1])
if numpy.ma.count_masked(wx) > len(wx)-2:
raise Exception("Could not find enough (>1) valid world coordinates in this map!")
wxmin = wx.min(); wxmax = wx.max(); wymin = wy.min(); wymax = wy.max()
# If one of the axes is a 'longitude' type axis then
# we can have 'jumpy' behaviour near 0.0. What we need then is a routine that
# finds the greatest contiguous sequence in a sorted array.
# Then we need to filter the NaN's from the list. A masked array cannot
# be sorted (yet), so compress it to an array with only valid positions
# and convert to a numpy array.
if self.gmap.types[0] == 'longitude':
wx = numpy.asarray(numpy.ma.masked_where(numpy.isnan(world[0]), world[0]).compressed())
wxmin, wxmax = self.__adjustlonrange(wx)
if self.gmap.types[1] == 'longitude':
wy = numpy.asarray(numpy.ma.masked_where(numpy.isnan(world[1]), world[1]).compressed())
wymin, wymax = self.__adjustlonrange(wy)
return wxmin, wxmax, wymin, wymax
def __init__(self, header=None, graticuledata=None, axnum=None,
wcstypes=None,
pxlim=None, pylim=None,
mixpix=None, spectrans=None, skyout=None, alter='',
wxlim=None, wylim=None,
boxsamples=5000,
startx=None, starty=None, deltax=None, deltay=None,
skipx=False, skipy=False,
gridsamples=1000,
labelsintex=True,
offsetx=None, offsety=None,
unitsx=None, unitsy=None):
#-----------------------------------------------------------
"""
Purpose: Creates an object that defines a graticule
See class documentation above.
"""
#-----------------------------------------------------------
self.ptype = 'Graticule'
if graticuledata is not None:
if header is None:
header = graticuledata.hdr
axnum = graticuledata.axperm
pxlim = graticuledata.pxlim
pylim = graticuledata.pylim
mixpix = graticuledata.mixpix
wcstypes = graticuledata.wcstypes
# Allow these to be overwritten
if spectrans is None:
spectrans = graticuledata.spectrans
if skyout is None:
skyout = graticuledata.skyout
if alter =='':
alter = graticuledata.alter
self.frame = None
self.frame2 = None
self.skyout = skyout
self.spectrans = spectrans
self.labelsintex = labelsintex
if wcstypes is None:
raise Exception("Need a list with wcs types for these axes")
self.wcstypes = wcstypes
# Try to get two axis numbers if none are given
if axnum is None:
naxis = header['NAXIS']
if naxis < 2:
raise Exception("Need data with at least two axes")
else:
self.xaxnum = 1
self.yaxnum = 2
else:
if len(axnum) != 2:
raise Exception("Need two axis numbers to create a graticule")
else:
self.xaxnum = axnum[0]
self.yaxnum = axnum[1]
if (self.xaxnum == self.yaxnum):
raise Exception("Need two different axis numbers")
# Get the axes limits in pixels. The default is copied from
# the values of NAXISn in the header. Note that there are no alternative
# keywords for NAXIS.
if pxlim is None:
self.pxlim = (1, header['NAXIS' +str(self.xaxnum)])
else:
if not issequence(pxlim):
raise Exception("pxlim needs to be of type tuple or list")
elif len(pxlim) != 2:
raise Exception("pxlim must have two elements")
else:
self.pxlim = pxlim
if pylim is None:
self.pylim = (1, header['NAXIS' +str(self.yaxnum)])
else:
if not issequence(pylim):
raise Exception("pylim needs to be of type tuple or list")
elif len(pylim) != 2:
raise Exception("pylim must have two elements")
else:
self.pylim = pylim
# wcs.debug=True
# Create the wcs projection object
proj = wcs.Projection(header, skyout=skyout, alter=alter)
if spectrans:
# Spectral transformations must be applied to the Projection object
# with all her axes because we need it both for plotting
# graticule label in the right spectral type (one of the image
# axes is spectral) or we want the output of the pixels on the
# repeat axis in the right spectral coordinates (if that axis is
# spectral). If we postpone this translation until we have a
# 2-D image, the spectral translation will not work for the
# repeat axis.
st2 = spectrans.split('-')
if len(st2) == 1:
spectrans += '-???'
proj = proj.spectra(spectrans)
# If one of the selected axes is a spatial axis and the other
# is not, then try to find the axis number of the matching axis.
mix = False
if self.xaxnum == proj.lonaxnum and self.yaxnum != proj.lataxnum:
self.matchingaxnum = proj.lataxnum
mix = True
elif self.xaxnum == proj.lataxnum and self.yaxnum != proj.lonaxnum:
self.matchingaxnum = proj.lonaxnum
mix = True
if self.yaxnum == proj.lonaxnum and self.xaxnum != proj.lataxnum:
self.matchingaxnum = proj.lataxnum
mix = True
elif self.yaxnum == proj.lataxnum and self.xaxnum != proj.lonaxnum:
self.matchingaxnum = proj.lonaxnum
mix = True
if mix:
if mixpix is None:
mixpix = proj.source['CRPIX'+str(self.matchingaxnum)+proj.alter]
if mixpix is None:
raise Exception("Could not find a grid for the missing spatial axis")
ok = proj.lonaxnum is not None and proj.lataxnum is not None
if not ok:
raise Exception("Could not find a matching spatial axis pair")
gmap = proj.sub([self.xaxnum, self.yaxnum, self.matchingaxnum])
self.mixpix = mixpix
else:
gmap = proj.sub([self.xaxnum, self.yaxnum])
self.mixpix = None
self.gmap = gmap
self.gmap.allow_invalid = True
# We arrived at the stage that we have a Projection object (called 'gmap')
# which represents the current projection system and/or spectral system.
# For maps with only one spatial axis, a matching axis is given. We already
# know the pixel limits for the axes in a map (pxlim, pylim).
# The positions module provides a useful function (str2pos) to convert
# strings to positions.
# Note that the getStarts function gets a Projection object called
# 'gmap' which has always a missing spatial axis added if the image
# has only one spatial axis.
startx, startpixX, errmes = getStarts(startx, self.gmap, 0, axnum, wcstypes[0], self.mixpix)
if errmes != '':
raise Exception(errmes)
starty, startpixY, errmes = getStarts(starty, self.gmap, 1, axnum, wcstypes[1], self.mixpix)
if errmes != '':
raise Exception(errmes)
# The next line is added at 28-08-2010 and is necessary because
# the sky system does not need to be an integer number anymore.
# It can also be a string. Many tests however are based on the integer
# number so we convert the sky system first to a number.
self.__skysys, refin, epochin, epobs = skyparser(proj.skyout)
# Now we have a projection object available and we want the limits of the axes in
# world coordinates. If nothing is specified for the constructor, we have
# to calculate estimates of these ranges.
if wxlim is not None:
if not issequence(wxlim):
raise Exception("wxlim needs to be of type tuple or list")
elif len(wxlim) != 2:
raise Exception("wxlim must have two elements")
else:
self.wxlim = wxlim
if wylim is not None:
if not issequence(wylim):
raise Exception("wylim needs to be of type tuple or list")
elif len(wylim) != 2:
raise Exception("wylim must have two elements")
else:
self.wylim = wylim
if wxlim is None or wylim is None:
if boxsamples < 2:
raise Exception("boxsamples < 2: Need at least two samples to find limits")
minmax = self.__estimateLonLatRanges(boxsamples)
if wxlim is None:
self.wxlim = (minmax[0], minmax[1])
if wylim is None:
self.wylim = (minmax[2], minmax[3])
# At this point we need to know for which constant positions we need to find a graticule
# line. We distinguish the following situations:
# A) The world coordinate is not spatial
# 1. A set of start values for the graticules are entered. The step (delta) is a dummy.
# 2. The start values are not given, but a delta is --> find suitable start values
# 3. Neither start values nor a delta is entered --> find suitable values for both
# 4. The user wants the labels to be plotted as offsets --> set output in offsets
# i.e. plot n*deltax/y instead of the world coordinate in 'startx/y'
#
# B) The world coordinate is spatial the other is not
# Set defaults. If only one coordinate is spatial then offsets are the default.
# 1) Values for startx/y are entered --> ignore
# 2) Value for deltax/y is entered --> use as delta and calculate world coordinates
# that correspond to the given delta.
# Note that this situation is not trivial. We label for a constant missing spatial axis
# which is given in grids. Then the world coordinate that corresponds to this grid
# needs not to be constant because the other spatial coordinate varies. We use
# methods from the ruler class to find world coordinates that represent
# these offsets as distances on a sphere.
#
# C) Both world coordinates are spatial
# The default should not be an offset. If the user overrides this, then offsets
# are offsets on the graticule and do not represent distances on a sphere.
# Check the offsets
self.offsetx = offsetx
self.offsety = offsety
spatialx = self.gmap.types[0] in ['longitude', 'latitude']
spatialy = self.gmap.types[1] in ['longitude', 'latitude']
if self.offsetx is None:
self.offsetx = spatialx and not spatialy
if self.offsety is None:
self.offsety = spatialy and not spatialx
self.prec = [0, 0]
# Handle the values for the step in the labels/graticule lines.
self.delta = [None, None]
axisunits = self.gmap.units[0]
if not deltax is None and isinstance(deltax, str):
# Not empty and the type is a string. Then a unit can follow.
parts = deltax.split()
if len(parts) > 1:
uf, errmes = unitfactor(parts[1], axisunits)
if uf is None:
raise ValueError(errmes)
else:
uf = 1.0
deltax = uf * eval(parts[0])
axisunits = self.gmap.units[1]
if not deltay is None and isinstance(deltay, str):
# Not empty and the type is a string. Then a unit can follow.
parts = deltay.split()
if len(parts) > 1:
uf, errmes = unitfactor(parts[1], axisunits)
if uf is None:
raise ValueError(errmes)
else:
uf = 1.0
deltay = uf * eval(parts[0])
self.offsetvaluesx = None
self.offsetvaluesy = None
self.funx = self.funy = None
self.fmtx = self.fmty = None
self.radoffsetx = self.radoffsety = False
if self.offsetx and spatialx:
# Then the offsets are distances on a sphere
xmax = self.pxlim[1] + 0.5
xmin = self.pxlim[0] - 0.5
ymin = self.pylim[0] - 0.5
x1 = xmin; x2 = xmax; y1 = y2 = ymin
Lambda = getlambda(startpixX, xmin, xmax)
rulerx = self.Ruler(x1=x1, y1=y1, x2=x2, y2=y2, lambda0=Lambda, step=deltax, units=unitsx)
#print "WCSGRAT ruler xw, yw", rulerx.xw, rulerx.yw, rulerx.stepsizeW, rulerx.offsets
#print "WCSGRAT ruler x, y", rulerx.x, rulerx.y
self.xstarts = rulerx.xw
self.prec[0] = 0
self.delta[0] = rulerx.stepsizeW
self.offsetvaluesx = rulerx.offsets
self.funx = rulerx.fun
self.fmtx = rulerx.fmt
self.radoffsetx = True
elif issequence(startx):
self.xstarts = startx
if len(startx) >= 2:
self.delta[0] = startx[1] - startx[0] # Assume this delta also for not equidistant values in startx
else:
# startx is a scalar or None
if startx is None and self.offsetx:
startx = (self.wxlim[1] + self.wxlim[0]) / 2.0
self.xstarts, self.offsetvaluesx, self.prec[0], self.delta[0] = self.__nicenumbers(self.wxlim[0],self.wxlim[1],
start=startx,
delta=deltax,
axtype=self.gmap.types[0],
skysys=self.__skysys)
if self.offsety and spatialy:
stdout.flush()
# Then the offsets are distances on a sphere
ymax = self.pylim[1] + 0.5
xmin = self.pxlim[0] - 0.5
ymin = self.pylim[0] - 0.5
x1 = xmin; x2 = xmin; y1 = ymin; y2 = ymax
Lambda = getlambda(startpixY, ymin, ymax)
rulery = self.Ruler(x1=x1, y1=y1, x2=x2, y2=y2, lambda0=Lambda, step=deltay, units=unitsy)
if spatialx and spatialy:
self.ystarts = rulery.yw
else:
self.ystarts = rulery.xw # NOTE: do not change in yw!
self.prec[1] = 0
self.delta[1] = rulery.stepsizeW
self.offsetvaluesy = rulery.offsets
self.funy = rulery.fun
self.fmty = rulery.fmt
self.radoffsety = True
elif issequence(starty):
self.ystarts = starty
if len(starty) >= 2:
self.delta[1] = starty[1] - starty[0] # Assume this delta also for not equidistant values in startx
else:
# starty is a scalar
if starty is None and self.offsety:
starty = (self.wylim[1] + self.wylim[0]) / 2.0
self.ystarts, self.offsetvaluesy, self.prec[1], self.delta[1] = self.__nicenumbers(self.wylim[0],self.wylim[1],
start=starty,
delta=deltay,
axtype=self.gmap.types[1],
skysys=self.__skysys)
# We have two sets of lines. For spatial maps, these are the
# meridians and the parallels. If other axes are involved, it is
# better to use the name grid line with constant x or constant
# y to identify the graticule lines.
if (gridsamples < 2):
raise Exception("Number of samples along graticule line must be >= 2 to avoid a step size of zero")
# Create the plot axes. The defaults are: plot native ticks to axis
# for the left and bottom axis and omit ticks along right and top axis.
# Get the equinox so that it can be inserted into the axis title.
# This equinox has been parsed earlier with skyparser()
epoch = str(epochin)
annot = ['']*2
for aa in [0,1]:
units = self.gmap.units[aa]
if aa == 0 and unitsx is not None:
units = unitsx
if aa == 1 and unitsy is not None:
units = unitsy
if (aa == 0 and self.offsetx) or (aa == 1 and self.offsety):
annot[aa] = "Offset "
if self.gmap.types[aa] in [None, 'spectral']:
annot[aa] += str(self.gmap.ctype[aa]).split('-')[0] + ' (' + str(units) + ')'
else: # Must be spatial
if (aa == 0 and (self.radoffsetx or self.offsetx)) or\
(aa == 1 and (self.radoffsety or self.offsety)):
if self.gmap.types[aa] == 'longitude':
olab = 'Radial offset lon.'
else:
olab = 'Radial offset lat.'
if (aa == 0):
if unitsx:
olab += '('+unitsx+')'
annot[aa] = olab
else:
if unitsy:
olab += '('+unitsy+')'
annot[aa] = olab
else:
if self.gmap.types[aa] == 'longitude':
if self.__skysys == wcs.equatorial:
annot[aa] = 'R.A.' + ' (' + epoch + ')'
elif self.__skysys == wcs.ecliptic:
annot[aa] = 'Ecliptic longitude' + ' (' + epoch + ')'
elif self.__skysys == wcs.galactic:
annot[aa] = 'Galactic longitude'
elif self.__skysys == wcs.supergalactic:
annot[aa] = 'Supergalactic longitude'
else:
if self.__skysys == wcs.equatorial:
annot[aa] = 'Dec.' + ' (' + epoch + ')'
elif self.__skysys == wcs.ecliptic:
annot[aa] = 'Ecliptic latitude' + ' (' + epoch + ')'
elif self.__skysys == wcs.galactic:
annot[aa] = 'Galactic latitude'
elif self.__skysys == wcs.supergalactic:
annot[aa] = 'Supergalactic latitude'
self.graticule = [] # Initialize the list with graticule lines
if not skipx:
axisunits = self.gmap.units[0]
if not unitsx is None:
units = unitsx
uf, errmes = unitfactor(axisunits, units)
if uf is None:
raise ValueError(errmes)
fie = lambda x: x*uf
fmt = "%g"
self.funx = fie
self.fmtx = fmt
for i, x in enumerate(self.xstarts):
offsetlabel = None
fie = fmt = None
if self.radoffsetx:
offsetlabel = self.offsetvaluesx[i]
fie = self.funx
fmt = self.fmtx
elif self.offsetx:
offsetlabel = self.offsetvaluesx[i]
if self.gmap.types[0] in ['longitude', 'latitude'] and self.labelsintex:
fmt = "%g^{\circ}"
else:
fmt = "%g"
elif not unitsx is None:
fie = self.funx
fmt = self.fmtx
gridl = Gratline(0, x, self.gmap,
self.pxlim, self.pylim,
self.wxlim, self.wylim,
gridsamples,self.mixpix, self.__skysys,
offsetlabel=offsetlabel,
fun=fie, fmt=fmt)
# Special case. Both axes are spatial but a user required
# to label the x axis in offsets. The graticule line however is
# calculated using the value of 'startx' and a non constant value
# for the other axis (e.g. the declination axis). This is the way
# we calculate the usual graticule lines, but we required an offset
# which is not constant along the graticule line. As a trick we
# replace the pixel coordinates of the x coordinate, by the x
# value of the ruler. This value is the pixel value of 'startx'.
# If this graticule line is curved and does not cross e.g. the top
# x axis, then we will see a partial graticule line. Usually
# this will be an exception.
if offsetx and spatialx and spatialy:
xoff0 = rulerx.x[i]
for lp in gridl.linepieces:
for ll in range(len(lp[0])):
lp[0][ll] = xoff0
self.graticule.append(gridl)
gridl.kwargs = {'color': '0.75', 'lw': 1}
if not skipy:
axisunits = self.gmap.units[1]
if not unitsy is None:
units = unitsy
uf, errmes = unitfactor(axisunits, units)
if uf is None:
raise ValueError(errmes)
fie = lambda x: x*uf
fmt = "%g"
self.funy = fie
self.fmty = fmt
for i, y in enumerate(self.ystarts):
offsetlabel = None
fie = fmt = None
if self.radoffsety:
offsetlabel = self.offsetvaluesy[i]
fie = self.funy
fmt = self.fmty
elif self.offsety:
offsetlabel = self.offsetvaluesy[i]
if self.gmap.types[1] in ['longitude', 'latitude'] and self.labelsintex:
fmt = "%g^{\circ}"
else:
fmt = "%g"
elif not unitsy is None:
fie = self.funy
fmt = self.fmty
gridl = Gratline(1, y, self.gmap,
self.pxlim, self.pylim,
self.wxlim, self.wylim,
gridsamples,self.mixpix, self.__skysys,
offsetlabel=offsetlabel,
fun=fie, fmt=fmt)
if offsety and spatialx and spatialy:
yoff0 = rulery.y[i]
for lp in gridl.linepieces:
for ll in range(len(lp[1])):
lp[1][ll] = yoff0
gridl.kwargs = {'color': '0.75', 'lw': 1}
self.graticule.append(gridl)
# Set properties for the rectangle axes.
xnumticks = 0
ynumticks = 0
for gridline in self.graticule:
wcsaxis = gridline.wcsaxis
for t in gridline.ticks:
anr = t.axisnr
if wcsaxis == 1 and anr == 0:
ynumticks += 1
if wcsaxis == 0 and anr == 1:
xnumticks += 1
x1mode = 0
x2mode = 3
y1mode = 0
y2mode = 3
if xnumticks < 2:
x1mode = 2
x2mode = 2
if ynumticks < 2:
y1mode = 2
y2mode = 2
axes = []
# keyword arguments for the tick labels
kwargs1 = {'fontsize':10}
kwargs2 = {'fontsize':10}
kwargs3 = {'fontsize':10, 'rotation':270, 'visible':False}
kwargs4 = {'fontsize':10, 'visible':False}
axes.append( PLOTaxis(left, y1mode, annot[1], **kwargs1) )
axes.append( PLOTaxis(bottom, x1mode, annot[0], **kwargs2) )
axes.append( PLOTaxis(right, y2mode, annot[1], **kwargs3) )
axes.append( PLOTaxis(top, x2mode, annot[0], **kwargs4) )
self.axes = axes
dummy = self.get_aspectratio() # Set default values for aspect ratio
self.objlist = [] # Initialize list with Graticules & derived objects
def plot(self, framebase, frame=None, frame2=None):
#-----------------------------------------------------------
"""
Purpose: Plot the graticule lines and labels
Labels are either along the plot axes or
inside the plot.
Parameters:
framebase An Axes object to which the graticules are
scaled. Two new Axes objects are created to hold two
independent graticules.
frame If not None, this is an Axes object, previously used
to hold a graticule. So one can re-use this frame,
which is in fact a requirement if one wants to
zoom/dezoom images with graticule and be able
to restore different zoom stages with the toolbar buttons.
frame2 See frame.
Returns: --
Notes: In Sep 2012 we changed this method to be able to
interactively zoom an image with a graticule.
There is no method to (re)set the pixel limits of
a graticule object, so for zooming one should
destroy an old graticule and make a new one.
This is not an ideal situation, but it works.
However it is often dangerous to add Axes objects while
zooming, so we added the option to reuse existing frames.
We copy an example from module maputils:
frame1 = cube.grat.frame
frame2 = cube.grat.frame2
cube.grat = currentimage.Graticule(pxlim=newpxlim, pylim=newpylim)
cube.grat.plot(currentimage.frame, frame1, frame2)
TODO: We should check other plot methods. If they are
also used in zoom/dezoom/pan actions, then they should also
have the option to reuse existing frames.
29 Jul 2019, VOG: Changed routines be able to
set ticks inwards or outwards
with keyword 'direction='
"""
#-----------------------------------------------------------
# We need to sort and format the ticks for the 4 plot axes
graticule = self
tex = graticule.labelsintex
pos, lab, kwargs, Mkwargs = graticule.sortticks(tex=tex)
aspect = framebase.get_aspect()
adjust = framebase.get_adjustable()
# Each graticule holds two frames for two times two axes
if not frame:
framelabel= 'FR'+"".join([choice(letters) for i in range(8)])
frame = framebase.inset_axes([0,0,1,1],
aspect=aspect,
adjustable=adjust,
autoscale_on=False,
frameon=False,
label=framelabel)
"""
try:
# VOG, 06-07-2023: Deprecated: r,c,n = framebase.get_geometry()
r,c,n = framebase.get_subplotspec().get_topmost_subplotspec().get_geometry()[:3]
n += 1
frame = framebase.figure.add_subplot(r, c, n,
aspect=aspect,
adjustable=adjust,
autoscale_on=False,
frameon=False,
label=framelabel)
# If a colorbar is requested then the original frame (framebase) is
# altered. The add_subplot method does not see this altered position.
frame.set_position(framebase.get_position(original=True))
except:
frame = framebase.figure.add_axes(framebase.get_position(original=True),
aspect=aspect,
adjustable=adjust,
autoscale_on=False,
frameon=False,
label=framelabel)
"""
graticule.frame = frame # Store frame so that frame can be identified with events
# In Matplotlib, it is possible that the upper pixel limit is
# lower than the lower value. To get graticules for entire pixels in your image
# we add 0.5, but as an extension. So if hi < lo, we should change the sign of
# the delta we want to add. However... we work with FITS files or FITS related
# files where always hi > lo.
# This is important for the following special case.
# Assume a user zoomed an area and the area does not include enire border
# pixels. Then we have a different situation than that of standard image limits
# which are always in integer pixels. Then one has to avoid the 0.5 correction
# by adding/subtracting this value before calling the constructor.
# Assume now that we have an area smaller than 1 pixel. Then this correction
# is the cause that hi < lo and we add the wrong delta's which increases
# the area instead of keeping it constant. So we omit the special case
# that hi < lo.
deltax = deltay = 0.5
"""
if graticule.pxlim[0] <= graticule.pxlim[1]:
deltax = 0.5
else:
deltax = -0.5
if graticule.pylim[0] <= graticule.pylim[1]:
deltay = 0.5
else:
deltay = -0.5
"""
xlo = graticule.pxlim[0]-deltax; ylo = graticule.pylim[0]-deltay
xhi = graticule.pxlim[1]+deltax; yhi = graticule.pylim[1]+deltay
# ==== LEFT ====
frame.set_yticks(pos[left])
frame.set_yticklabels(lab[left])
for tick, kw, mkw in zip(frame.yaxis.get_major_ticks(), kwargs[left], Mkwargs[left]):
if len(kw) != 0:
pad = kw.get("pad", None)
if pad:
fkw = {i:kw[i] for i in kw if i!='pad'} # Filter pad key, preserve kw
tick.set_pad(pad)
tick.label1.set(**fkw)
else:
tick.label1.set(**kw)
direction = 'in'
if len(mkw) != 0:
direction = mkw.get("direction", None)
if direction:
fmkw = {i:mkw[i] for i in mkw if i!='direction'} # Filter direction key, preserve mkw
tick.tick1line.set(**fmkw)
else:
tick.tick1line.set(**mkw)
# VOG, 06-07-2023:
tick._apply_tickdir(direction)
#frame.yaxis.set_tick_params(direction=direction)
tick.tick2on = False
tick.tick2line.set_visible(False)
# VOG 06-07-2023 not allowed anymore: tick.tick1line.set_marker(tick._tickmarkers[0]) # Hope this attribute will remain in newer releases
# ==== BOTTOM ====
frame.set_xticks(pos[bottom])
frame.set_xticklabels(lab[bottom])
dy12 = None
for tick, kw, mkw in zip(frame.xaxis.get_major_ticks(), kwargs[bottom], Mkwargs[bottom]):
if len(kw) != 0:
pad = kw.get("pad", None)
if pad:
fkw = {i:kw[i] for i in kw if i!='pad'} # Filter pad key, preserve kw
tick.set_pad(pad)
tick.label1.set(**fkw)
elif not ('va' in kw or 'verticalalignment' in kw):
tick.label1.set(va="baseline", **kw)
# VOG: 11-07-2023
# Tick labels for the bottom axis are not properly aligned when superscripts
# (hms) are used. We can correct this by either using usetex=True or one
# can change the vertical alignment. But with this solution we don't have
# control over the y position so we need a calculation for that. But only
# if 'y' is not in **kw
if not 'y' in kw:
if dy12 is None:
r = frame.figure.canvas.get_renderer()
bb = tick.label1.get_window_extent(renderer=r)
inv = frame.transAxes.inverted()
yt1 = inv.transform((0, bb.height))[1]
yt2 = inv.transform((0, 0))[1]
dy12 = abs(yt1-yt2)
tick.label1.set_y(-dy12)
direction = 'in'
if len(mkw) != 0:
direction = mkw.get("direction", None)
if direction:
fmkw = {i:mkw[i] for i in mkw if i!='direction'} # Filter direction key, preserve mkw
tick.tick1line.set(**fmkw)
else:
tick.tick1line.set(**mkw)
tick._apply_tickdir(direction)
tick.tick2on = False
# VOG 06-07-2023, tick.tick1line.set_marker(tick._tickmarkers[0])
#for line in (tick.tick1line,): #: VOG Turned out not to be necessary
# line.set_markersize(tick.tick1line.size)
tick.tick2line.set_visible(False)
if not frame2:
framelabel= 'SFR'+"".join([choice(letters) for i in range(8)])
frame2 = framebase.inset_axes([0,0,1,1],
aspect=aspect,
adjustable=adjust,
autoscale_on=False,
frameon=False,
label=framelabel)
"""
try:
# VOG, 06-07-2023: Deprecated: r,c,n = framebase.get_geometry()
r,c,n = framebase.get_subplotspec().get_topmost_subplotspec().get_geometry()[:3]
n += 1
frame2 = framebase.figure.add_subplot(r, c, n,
aspect=aspect,
adjustable=adjust,
autoscale_on=False,
frameon=False,
label=framelabel)
frame2.set_position(framebase.get_position(original=True))
except:
frame2 = framebase.figure.add_axes(frame.get_position(original=True),
aspect=aspect,
adjustable=adjust,
autoscale_on=False,
frameon=False,
label=framelabel)
"""
# frame2 = framebase.figure.add_axes(frame.get_position(), frameon=False, label=framelabel)
graticule.frame2 = frame2 # Add new attribute
# axis sharing is not an option because then also the ticks are
# shared and we want independent ticks along all 4 axes. For most
# projections the axes are not related.
frame2.yaxis.set_label_position('right')
frame2.xaxis.set_label_position("top")
frame2.yaxis.tick_right()
frame2.xaxis.tick_top()
frame2.set_xlim(xlo,xhi)
frame2.set_ylim(ylo,yhi)
frame2.set_aspect(aspect=aspect, adjustable=adjust)
frame2.set_xticks(pos[top])
frame2.set_xticklabels(lab[top])
# ==== TOP ====
for tick, kw, mkw in zip(frame2.xaxis.get_major_ticks(),kwargs[top], Mkwargs[top]):
if len(kw) != 0:
pad = kw.get("pad", None)
if pad:
fkw = {i:kw[i] for i in kw if i!='pad'} # Filter pad key, preserve kw
tick.set_pad(pad)
tick.label2.set(**fkw)
else:
tick.label2.set(**kw)
direction = 'in'
if len(mkw) != 0:
direction = mkw.get("direction", None)
if direction:
fmkw = {i:mkw[i] for i in mkw if i!='direction'} # Filter direction key, preserve mkw
tick.tick2line.set(**fmkw)
else:
tick.tick2line.set(**mkw)
# tick.apply_tickdir(direction)
#frame2.xaxis.set_tick_params(direction=direction)
# VOG 06-07-2023, tick.tick2line.set_marker(tick._tickmarkers[1])
tick.tick1line.set_visible(False)
# ==== RIGHT ====
frame2.set_yticks(pos[right])
frame2.set_yticklabels(lab[right])
for tick, kw, mkw in zip(frame2.yaxis.get_major_ticks(),kwargs[right], Mkwargs[right]):
if len(kw) != 0:
pad = kw.get("pad", None)
if pad:
fkw = {i:kw[i] for i in kw if i!='pad'} # Filter pad key, preserve kw
tick.set_pad(pad)
tick.label2.set(**fkw)
else:
tick.label2.set(**kw)
direction = 'in'
if len(mkw) != 0:
direction = mkw.get("direction", None)
if direction:
fmkw = {i:mkw[i] for i in mkw if i!='direction'} # Filter direction key, preserve mkw
tick.tick2line.set(**fmkw)
else:
tick.tick2line.set(**mkw)
# tick.apply_tickdir(direction)
#frame2.yaxis.set_tick_params(direction=direction)
# VOG 06-07-2023, tick.tick2line.set_marker(tick._tickmarkers[1])
tick.tick1line.set_visible(False)
# Adjust label position if one is given
for anr in [left, right]:
xpos = graticule.axes[anr].xpos
ypos = graticule.axes[anr].ypos
if xpos is not None or ypos is not None:
if xpos is None and ypos is not None:
graticule.axes[anr].kwargs.update({'y':ypos})
else:
if ypos is None: ypos = 0.5
# xpos is not None then use method
if anr == left:
fr = frame
else:
fr = frame2
fr.yaxis.set_label_coords(xpos, ypos)
for anr in [bottom, top]:
xpos = graticule.axes[anr].xpos
ypos = graticule.axes[anr].ypos
if xpos is not None or ypos is not None:
if ypos is None and xpos is not None:
graticule.axes[anr].kwargs.update({'x':xpos})
else:
if xpos is None: xpos = 0.5
# ypos is not None then use method set_label_coords
# because setting y as attribute has no effect.
if anr == bottom:
fr = frame
else:
fr = frame2
fr.xaxis.set_label_coords(xpos, ypos)
frame.set_ylabel( graticule.axes[left].label, **graticule.axes[left].kwargs)
frame.set_xlabel( graticule.axes[bottom].label, **graticule.axes[bottom].kwargs)
frame2.set_ylabel(graticule.axes[right].label, **graticule.axes[right].kwargs)
frame2.set_xlabel(graticule.axes[top].label, **graticule.axes[top].kwargs)
# Plot the line pieces
for gridline in graticule.graticule:
for line in gridline.linepieces:
#ppx = numpy.ma.masked_where(numpy.isfinite(line[0]), line[0])
#ppy = numpy.ma.masked_where(numpy.isfinite(line[1]), line[1])
#ppx = numpy.ma.masked_where(numpy.isinf(line[0]), line[0])
#ppy = numpy.ma.masked_where(numpy.isinf(line[1]), line[1])
frame.plot(line[0], line[1], **gridline.kwargs)
#frame.plot(ppx, ppy, **gridline.kwargs)
# set the limits of the plot axes
# this setting can be overwritten in the calling environment
frame.set_xlim((xlo,xhi))
frame.set_ylim((ylo,yhi))
frame.set_aspect(aspect=aspect, adjustable=adjust)
# Plot the objects derived from this graticule
if len(graticule.objlist) > 0:
for obj in graticule.objlist:
try:
pt = obj.ptype
except:
raise Exception("Unknown object. Cannot plot this!")
if pt in ["Insidelabels"]:
try:
visible = obj.visible
except:
visible = True
# Previously: obj.plot(framebase)
# We want the labels on top of the graticule lines.
# Then they must be plotted in the frame of the inside labels
# and not in the base frame, because that frame is not on top.
obj.plot(frame2)
framebase.figure.sca(framebase) # back to frame from calling environment
[docs] def scanborder(self, xstart, ystart, deltax=None, deltay=None, nxy=1000, tol=None):
#-----------------------------------------------------------------
"""
For the slanted azimuthal projections, it is
not trivial to draw a border because these
borders are not graticule lines with a constant
longitude or constant latitude. Nor it is
easy or even possible to find mathematical
expressions for this type of projection.
Also, the mathematical expressions return
world coordinates which can suffer from loss
of precision.
This method tracks the border from a starting
point by scanning in x- and y direction and
tries to find the position of a limb with a
standard bisection technique. This method has been applied
to a number of all-sky plots with slanted projections.
:param xstart: X-coordinate in pixels of position where to
start the scan to find a border.
The parameter has no default.
:type xstart: Floating point
:param ystart: Y-coordinate in pixels of position where to
start the scan to find border.
The parameter has no default.
:type ystart: Floating point
:param deltax: Set range in pixels to look for a border in
scan direction. The default value is 10 percent
of the total pixel range in x- or y-direction.
:type deltax: Floating point
:param deltay: See *deltayx*.
:type deltay: Floating point
:param nxy: Number of scan lines in x and y direction.
Default is 1000.
:type nxy: Integer
:param tol: See note below.
:type tol: Floating point
:returns: Identifier to set attributes of this
graticule line with method :meth:`setp_linespecial()`.
:note: This method uses an algorithm to find
positions along the border of a projection.
It scans along both x- and y-axis for
a NaN (Not a Number number) transition as a result
of an invalid coordinate transformation,
and repeats this for a number of scan lines
along the x-axis and y-axis.
::
A position on a border off an all-sky plot is the position at
which a transition occurs from a valid coordinate to a NaN.
Its accuracy depends on the the tolerance
given in argument *tol*.
The start coordinates to find the next border
position on the next scan line is the
position of the previous border point.
If you have missing line pieces, then add more
borders by calling this method with different
starting points.
"""
#-----------------------------------------------------------------
xp = []; yp = []
if deltax is None:
deltax = (self.pxlim[1] - self.pxlim[0])/ 10.0
if deltay is None:
deltay = (self.pylim[1] - self.pylim[0])/ 10.0
d = (float(self.pxlim[1] - self.pxlim[0]), float(self.pylim[1] - self.pylim[0]))
delta = (deltax, deltay)
limlo = (self.pxlim[0], self.pylim[0])
limhi = (self.pxlim[1], self.pylim[1])
start = (xstart, ystart)
for i in [0,1]:
if tol is None:
tol = delta[i] / 1000.0
nx1 = (start[i] - limlo[i])/d[i]
nx2 = 1.0 - nx1
nx1 = int(nxy*nx1)
nx2 = int(nxy*nx2)
l1 = numpy.linspace(start[i], limlo[i], nx1)
l2 = numpy.linspace(start[i], limhi[i], nx2)
X = numpy.concatenate((l1,l2))
Y0 = (ystart, xstart)
for xb in X:
if xb == start[i]:
y0 = Y0[i]
yb = self.__bisect(i, xb, y0-delta[i]/2.0, y0+delta[i]/2.0, self.gmap, tol)
if yb is not None:
if i == 0:
xp.append(xb)
yp.append(yb)
else:
xp.append(yb)
yp.append(xb)
y0 = yb
return self.addgratline(xp, yp, pixels=True)
[docs] def addgratline(self, x, y, pixels=False):
#-----------------------------------------------------------------
"""
For any path given by a set of world coordinates
of which none is a constant value (e.g. borders
in slanted projections where the positions are calculated by an
external routine),
one can create a line that is processed as a graticule
line, i.e. intersections and jumps are addressed.
Instead of world coordinates, this method
can also process pixel positions. The type of input is set by the
*pixels* parameter.
:param x: A sequence of world coordinates or pixels
that correspond to the horizontal axis in a graticule plot..
:type x: Floating point numbers
:param y: The same for the second axis
:type x: Floating point numbers
:param pixels: False or True
If False the coordinates in x and y are world-
coordinates. Else they are pixel coordinates.
:type pixels: Boolean
:Returns: A Identification number *id* which can be used
to set properties for this special path with
method :meth:`setp_linespecial`.
Return *None* if no line piece could be found
inside the pixel limits of the graticule.
:note: This method can be used to plot a border
around an all-sky plot e.g. for slanted
projections. See code at :meth:`scanborder`.
"""
#-----------------------------------------------------------------
if len(x) > 0:
samples = len(x)
# Create an unique id for this line
id = len(self.graticule) + 2; # Avoid problems if there are no gridlines yet, so add at least 2 to the id value
gridl = Gratline(id, '', self.gmap,
self.pxlim, self.pylim,
self.wxlim, self.wylim,
samples,
self.mixpix, self.__skysys,
addx=x, addy=y, addpixels=pixels)
self.graticule.append(gridl)
gridl.kwargs = {'color': 'r', 'lw': 1}
else:
id = None;
return id
def __str__(self):
"""
#-----------------------------------------------------------
Purpose: Show textual contents of graticule:
lines, ticks and line pieces.
Parameters: -
Returns: -
Notes: The information is unsorted w.r.t. the
plot axis number.
Example: g = Graticule(header)
print(g)
#-----------------------------------------------------------
"""
s = ''
for gridline in self.graticule:
s += "\nWCS graticule line number %s\n" % gridline.wcsaxis
s += "Axis type: %s. Sky system %s:\n" % (gridline.axtype, gridline.skysys)
for t in gridline.ticks:
s += "tick x,y: %f %f\n" % (t.x, t.y)
s += "tick label: %f\n" % t.labval
s += "tick on axis: %d\n" % t.axisnr
for line in gridline.linepieces:
s += "line piece has %d elements\n" % len(line[0])
# line is (line[0], line[1])
return s
[docs] def get_aspectratio(self, xcm=None, ycm=None):
#-----------------------------------------------------------------
"""
Calculate and set, the aspect ratio for the current
pixels. Also set default values for figure
size and axes lengths (i.e. size of canvas depends
on the size of plot window with this aspect ratio).
:param xcm: Given a value for xcm or ycm (or omit both),
suggest a suitable figure size in and a viewport in
normalized device coordinates of a plot which has
an axes rectangle that corrects the figure for an
aspect ratio (i.e. CDELTy/CDELTx) unequal to 1 while
the length of the x-axis is xcm OR the length of the
y-axis is ycm. See note for non-spatial maps.
:type xcm: Floating point number
:param ycm: See description at *xcm*.
:type ycm: Floating point number
:Returns: The aspect ratio defined as: ``AR = CDELTy/CDELTx``.
:Note: (i.e. AR > 10 or AR < 0.1), an aspect ratio of 1
is returned. This method sets the attributes:
'axesrect', 'figsize', 'aspectratio'.
The attribute 'figsize' is in inches which is compatible
to the methods of Matplotlib.
"""
#-----------------------------------------------------------------
cdeltx = self.gmap.cdelt[0]
cdelty = self.gmap.cdelt[1]
nx = float(self.pxlim[1] - self.pxlim[0] + 1)
ny = float(self.pylim[1] - self.pylim[0] + 1)
if xcm is None and ycm is None:
xcm = 20.0
aspectratio = abs(cdelty/cdeltx)
if aspectratio > 10.0 or aspectratio < 0.1:
aspectratio = nx/ny
if xcm is None:
xcm = ycm
else:
ycm = xcm
if ycm is None:
ycm = xcm * (ny/nx) * aspectratio
if xcm is None:
xcm = ycm * (nx/ny) / aspectratio
fh = 0.7; fw = 0.7
self.axesrect = (0.15, 0.15, fw, fh)
self.figsize = (xcm/2.54/fw, ycm/2.54/fh)
self.aspectratio = aspectratio
return aspectratio
[docs] def setp_tick(self, wcsaxis=None, plotaxis=None, position=None,
tol=1e-12, fmt=None, fun=None, tex=None, texsexa=None,
markerdict={}, **kwargs):
#-----------------------------------------------------------------
"""
Set (plot) attributes for a wcs tick label.
A tick is identified by the type of grid line
it belongs to, and/or the plot axis for which
it defines an intersection and/or a position which
corresponds to the constant value of the graticule
line.
All these parameters are valid with none, one or
a sequence of values.
.. warning:: If no value for *wcsaxis*, *plotaxis* or *position* is entered
then this method applies the parameter setting on all the wcs axes.
:param wcsaxis:
Values are 0 or 1, corresponding to the
first and second world coordinate types.
Note that *wcsaxis=0* corresponds to the
first element in the axis permutation array given in
parameter *axnum*.
:type wcsaxis: *None*, 0, 1 or tuple with both
:param plotaxis:
Accepted values are 'None', 0, 1, 2, 3 or a
sequence of these numbers, to represent the left, bottom, right
and top axis of the enclosing rectangle that
represents the limits in pixel coordinates.
:type plotaxis: One or more integers between 0 and 3.
:param position:
Accepted are None, or one or more values representing
the constant value of the graticule line in
world coordinates. These positions are used to identify
individual graticule lines so that each line can have its
own properties.
The input can also be a string that represents a sexagesimal number.
:type position: *None* or one or a sequence of floating point numbers
:param tol:
If a value > 0 is given, the gridline with the
constant value closest to a given position within
distance 'tol' gets updated attributes.
:type tol: Floating point number
:param fmt:
A string that formats the tick value
e.g. ``fmt="%10.5f"`` in the Python way, or a string
that contains no percentage character (%) but a format
to set the output of sexagesimal numbers e.g.
fmt='HMs'. The characters in the format either force
(uppercase) a field to be printed, or it suppresses
(lowercase) a field to be printed.
See also the examples at :func:`makelabel`.
To create labels with an exponential, use a second format in
the same format string. The syntax is %nne where nn is an
integer. This integer, which can be negative, sets the number
in the exponential. The number before the exponential is
formatted in the usual way e.g. fmt='%.3f%-3e'.
:type fmt: String
:param fun:
An external function which will be used to
convert the tick value e.g. to convert
velocities from m/s to km/s. See also
example 2_ below.
:type fun: Python function or Lambda expression
:param tex:
Interpret the format in *fmt* as a TeX label.
The default is set to *None* to indicate it has
not been set (to True or False) so that it is possible
to distinguish between global and local settings of this
property.
:type tex:
Boolean
:param texsexa:
If False and parameter *tex* is True, then format the
tick label without superscripts for sexagesimal labels.
This option can be used if superscripts result in 'jumpy' labels.
The reason is that in Matplotlib the TeX labels at the bottom
of a plot are aligned at a baseline at the top of the characters.
:type tex: Boolean
:param markerdict:
Properties for the tick marker. Amongst others:
* markersize:
Size of tick line.
* direction:
Use keyword argument 'direction='
get tick lines that point outside the plot instead
of the default which is inside.
* Set the tick label pad in points. Use a negative number to
shift a label inwards.
* markeredgewidth:
The width of the marker
* color: Color of the marker (not the label)
:type markerdict:
Python dictionary
:param `**kwargs`:
Keyword arguments for plot properties like *color*,
*visible*, *rotation* etc. The plot attributes are standard
Matplotlib attributes which can be found in the
Matplotlib documentation.
:type `**kwargs`: Matplotlib keyword arguments
:note:
Some projections generate labels that are very close
to each other. If you want to skip labels then you can
use keyword/value *visible=False*. There is not a documented
keyword *visible* in this method because *visible* is a
valid keyword argument in Matplotlib.
To set the direction of the tickmarks, use keyword argument
*direction=* with one of the options 'in' (default),
'out' or 'inout'. In the set with example code, one can find
examples with this keyword.
:Examples:
1. Set tick properties with :meth:`setp_tick`. The last line makes the
label at a declination of -10 degrees (we assume a spatial map) invisible::
grat.setp_tick(wcsaxis=0, color='g')
grat.setp_tick(wcsaxis=1, color='m')
grat.setp_tick(wcsaxis=1, plotaxis=wcsgrat.bottom,
color='c', rotation=-30, ha='left')
grat.setp_tick(plotaxis=wcsgrat.right, backgroundcolor='yellow')
grat.setp_tick(plotaxis=wcsgrat.left, position=-10, visible=False)
.. _2:
2. Example of an external function to change the values of the tick
labels for the horizontal axis only::
def fx(x):
return x/1000.0
setp_tick(wcsaxis=0, fun=fx)
Or use the lambda operator as in: ``fun=lambda x: x/1000``
.. _3:
3. Using the utility routines :meth:`setp_ticklabel()`
and :meth:`setp_tickmark()`::
graticule4.setp_tickmark(plotaxis="bottom", markersize=8, color='m',
markeredgewidth=3, direction='out')
graticule4.setp_ticklabel(plotaxis="bottom", rotation=20, color='y',
ha='right', va='top')
"""
#-----------------------------------------------------------------
if wcsaxis is None and plotaxis is None and position is None:
## Nothing to do
#return
wcsaxis = [0,1]
if wcsaxis is not None:
if not issequence(wcsaxis):
wcsa = [wcsaxis]
else:
wcsa = wcsaxis
if plotaxis is not None:
plta = parseplotaxes(plotaxis)
if position is not None:
if not issequence(position):
posn = [position]
else:
posn = position
for gridline in self.graticule:
skip = False
if (wcsaxis is not None):
skip = not (gridline.wcsaxis in wcsa)
if not skip:
if position is None:
for t in gridline.ticks:
skip = False
if plotaxis is not None:
skip = not (t.axisnr in plta)
if not skip:
t.kwargs.update(kwargs)
t.markkwargs.update(markerdict)
# The attributes fmt and fun to format tick labels could have
# been initialized when the object was created. If values are
# given then they override the defaults.
if fmt is not None: t.fmt = fmt
if fun is not None: t.fun = fun
if tex is not None: t.tex = tex
if not texsexa is None: t.texsexa = texsexa
# There are defaults for the precision in seconds for each axis
# If the user sets a format with precision in seconds (e.g. HMS.SS)
# then we need to tell this to different methods that deal with it.
if not fmt is None and fmt.find('%') == -1: # Not a common format, must be a HMS/DMS format
s2 = fmt.split('.')
if len(s2) > 1:
self.prec[gridline.wcsaxis] = len(s2[1])
else: # One or more positions are given. Find the right index.
for pos in posn:
if isinstance(pos, str):
C = pos.upper()
if 'H' in C or 'D' in C or 'M' in C or 'S' in C:
pos, err = parsehmsdms(pos)
if err != '':
raise Exception(err)
else:
pos = pos[0]
else:
raise ValueError("[%s] is entered as a string but does not represent valid HMS or DMS"%pos)
d0 = None
for i, t in enumerate(gridline.ticks):
skip = False
if plotaxis is not None:
skip = not (t.axisnr in plta)
if not skip:
d = abs(t.labval - pos)
if d <= tol:
if d0 is None:
d0 = d
indx = i
elif d < d0:
d0 = d
indx = i
if d0 is not None:
gridline.ticks[indx].kwargs.update(kwargs)
gridline.ticks[indx].markkwargs.update(markerdict)
if len(markerdict) == 0 or fmt is not None:
gridline.ticks[indx].fmt = fmt
if len(markerdict) == 0 or fun is not None:
gridline.ticks[indx].fun = fun
if tex is not None:
gridline.ticks[indx].tex = tex
if texsexa is not None:
gridline.ticks[indx].texsexa = texsexa
[docs] def setp_tickmark(self, wcsaxis=None, plotaxis=None, position=None,
tol=1e-12, **mkwargs):
#-----------------------------------------------------------------
"""
Utility method for :meth:`setp_tick`. It handles the
properties of the tick marks, which are :class:`Line2D` objects
in Matplotlib. The most useful properties are *color*, *markeredgewidth*
and *markersize*. Keyword direction= sets a tick direction to
inwards or outwards, or both ('in', 'out', 'inout').
See also documentation for :func`wcsgrat.setp_ticklabel`.
"""
#-----------------------------------------------------------------
self.setp_tick(wcsaxis=wcsaxis, plotaxis=plotaxis, position=position,
tol=tol, markerdict=mkwargs)
[docs] def setp_ticklabel(self, wcsaxis=None, plotaxis=None, position=None,
tol=1e-12, fmt=None, fun=None, tex=None, texsexa=None,
**kwargs):
#-----------------------------------------------------------------
"""
Utility method for :meth:`setp_tick`. It handles the
properties of the tick labels, which are :class:`Text` objects
in Matplotlib. The most useful properties are *color*, *fontsize*
and *fontstyle*.
:param wcsaxis:
Values are 0 or 1, corresponding to the
first and second world coordinate types.
Note that *wcsaxis=0* corresponds to the
first element in the axis permutation array given in
parameter *axnum*.
:type wcsaxis: *None*, 0, 1 or tuple with both
:param plotaxis:
Accepted values are 'None', 0, 1, 2, 3 or a
combination, to represent the left, bottom, right
and top axis of the enclosing rectangle that
represents the limits in pixel coordinates.
:type plotaxis: One or more integers between 0 and 3.
:param position:
Accepted are None, or one or more values representing
the constant value of the graticule line in
world coordinates. These positions are used to identify
individual graticule lines so that each line can have its
own properties.
The input can also be a string that represents a sexagesimal number.
:type position: *None* or one or a sequence of floating point numbers
:param tol:
If a value > 0 is given, the gridline with the
constant value closest to a given position within
distance 'tol' gets updated attributes.
:type tol: Floating point number
:param fmt:
A string that formats the tick value
e.g. ``fmt="%10.5f"`` in the Python way, or a string
that contains no percentage character (%) but a format
to set the output of sexagesimal numbers e.g.
fmt='HMs'. The characters in the format either force
(uppercase) a field to be printed, or it suppresses
(lowercase) a field to be printed.
See also the examples at :func:`makelabel`.
:type fmt: String
:param fun:
An external function which will be used to
convert the tick value e.g. to convert
velocities from m/s to km/s. See also
example 2_ below.
:type fun: Python function or Lambda expression
:param tex:
If True then format the tick label in LaTeX. This is the
default. If False then standard text will applies.
Some text properties cannot be changed if LaTeX is
in use.
:type tex: Boolean
:param texsexa:
If False and parameter *tex* is True, then format the
tick label without superscripts for sexagesimal labels.
This option can be used if superscripts result in 'jumpy' labels.
The reason is that in Matplotlib the TeX labels at the bottom
of a plot are aligned at a baseline at the top of the characters.
:type tex: Boolean
:param `**kwargs`:
Keyword arguments for plot properties like *color*,
*visible*, *rotation* etc. The plot attributes are standard
Matplotlib attributes which can be found in the
Matplotlib documentation.
:type `**kwargs`: Matplotlib keyword arguments
:note:
Some projections generate labels that are very close
to each other. If you want to skip labels then you can
use keyword/value *visible=False*. There is not a documented
keyword *visible* in this method because *visible* is a
valid keyword argument in Matplotlib.
"""
#-----------------------------------------------------------------
self.setp_tick(wcsaxis=wcsaxis, plotaxis=plotaxis, position=position,
tol=tol, fmt=fmt, fun=fun, tex=tex, **kwargs)
[docs] def setp_gratline(self, wcsaxis=None, position=None, tol=1e-12, **kwargs):
#-----------------------------------------------------------------
"""
Set (plot) attributes for one or more graticule
lines.
These graticule lines are identified by the wcs axis
number (*wcsaxis=0* or *wcsaxis=1*) and by their constant
world coordinate in *position*.
:param wcsaxis: If omitted, then for both types of graticule lines
the attributes are set.
If one value is given then only for that axis
the attributes will be set.
:type wcsaxis: *None* , integer or tuple with integers from set 0, 1.
:param position: None, one value or a sequence of
values representing the constant value of a graticule
line in world coordinates. For the graticule line(s)
that match a position in this sequence, the attributes
are updated.
:type position: *None*, one or a sequence of floating point numbers
:param tol: If a value > 0 is given, the graticule line with the
constant value closest to a given position within
distance *tol* gets updated attributes.
:type tol: Floating point number
:param `**kwargs`: Keyword arguments for plot properties like *color*,
*rotation* or *visible*, *linestyle* etc.
:type `**kwargs`: Matplotlib keyword argument(s)
:Returns: --
:Notes: For each value in *position* find the index of
the graticule line that belongs to *wcsaxis* so that
the distance between that value and the constant
value of the graticule line is the smallest of all
the graticule lines. If *position=None* then
apply change of properties to ALL graticule lines.
The (plot) properties are stored in `**kwargs`
Note that graticule lines are initialized with
default properties. These kwargs only update
the existing kwargs i.e. appending new keywords
and update existing keywords.
"""
#-----------------------------------------------------------------
# Upgrade 'wcsaxis' to a sequence
if wcsaxis is None:
wcsaxislist = [0,1]
elif not issequence(wcsaxis):
wcsaxislist = [wcsaxis]
else:
wcsaxislist = wcsaxis
if position is None:
for gridline in self.graticule:
if gridline.wcsaxis in wcsaxislist:
gridline.kwargs.update(kwargs)
else: # Upgrade 'position' to a sequence
if not issequence(position):
S = [position]
else:
S = position
for constval in S:
if isinstance(constval, str):
C = constval.upper()
if 'H' in C or 'D' in C or 'M' in C or 'S' in C:
constval, err = parsehmsdms(constval)
if err != '':
raise Exception(err)
else:
constval = constval[0]
else:
raise ValueError("[%s] is entered as a string but does not represent valid HMS or DMS"%constval)
d0 = None
for i, gridline in enumerate(self.graticule):
if gridline.wcsaxis in wcsaxislist:
d = abs(gridline.constval - constval)
if d <= tol:
if d0 is None:
d0 = d
indx = i
else:
if d < d0:
d0 = d
indx = i
if d0 is not None: # i.e. we found a closest position
self.graticule[indx].kwargs.update(kwargs)
[docs] def setp_lineswcs0(self, position=None, tol=1e-12, **kwargs):
#-----------------------------------------------------------------
"""
Helper method for :meth:`setp_gratline`.
It pre-selects the grid line that
corresponds to the first world coordinate.
:Parameters: See description at :meth:`setp_gratline`
:Examples: Make lines of constant latitude magenta and
lines of constant longitude green. The line that
corresponds to a latitude of 30 degrees and
the line that corresponds to a longitude of 0
degrees are plotted in red with a line width of 2::
grat.setp_lineswcs1(color='m')
grat.setp_lineswcs0(color='g')
grat.setp_lineswcs1(30, color='r', lw=2)
grat.setp_lineswcs0(0, color='r', lw=2)
"""
#-----------------------------------------------------------------
axis = 0
self.setp_gratline(axis, position, tol, **kwargs)
[docs] def setp_lineswcs1(self, position=None, tol=1e-12, **kwargs):
#-----------------------------------------------------------------
"""
Equivalent to method :meth:`setp_gratline`.
It pre-selects the grid line that
corresponds to the second world coordinate.
:Parameters: See description at :meth:`setp_gratline`
:Examples: See example at :meth:`setp_lineswcs0`.
"""
#-----------------------------------------------------------------
axis = 1
self.setp_gratline(axis, position, tol, **kwargs)
[docs] def setp_linespecial(self, id, **kwargs):
#-----------------------------------------------------------------
"""
Set (plot) attributes for a special type
of graticule line made with method :meth:`addgratline`
or method :meth:`scanborder`.
This graticule line has no constant x- or y- value.
It is identified by an id returned by method
:meth:`addgratline`.
:param id: id from :meth:`addgratline`
:type id: Integer
:param `**kwargs`: keywords for (plot) attributes
:type `**kwargs`: Matplotlib keyword argument(s)
:Examples: Create a special graticule line
which follows the positions in two
given sequences *x* and *y*. and set
the line width for this line to 2::
id = grat.addgratline(x, y)
grat.setp_linespecial(id, lw=2)
"""
#-----------------------------------------------------------------
# Set properties for an added gridline
for gridline in self.graticule:
if gridline.wcsaxis == id and id > 1:
gridline.kwargs.update(kwargs)
def switchdefaults(self):
self.setp_axislabel(plotaxis=("right","top"), visible=True)
self.setp_axislabel(plotaxis=("left","bottom"), visible=False)
self.set_tickmode(plotaxis=("right","top"), mode="NATIVE")
self.set_tickmode(plotaxis=("left","bottom"), mode="NO")
[docs] def setp_plotaxis(self, plotaxis, mode=None, label=None, xpos=None,
ypos=None, **kwargs):
#-----------------------------------------------------------------
"""
Set (plot) attributes for titles along a plot axis and set the ticks mode.
The ticks mode sets the relation between the ticks and the plot axis.
For example a rotated map will show a rotated graticule, so ticks for both
axes can appear along a plot axis. With parameter *mode* one can influence this
behaviour.
.. Note::
This method addresses the four axes of a plot separately. Therefore
its functionality cannot be incorporated in :meth:`setp_tick`
:param plotaxis: The axis number of one of the axes of the
plot rectangle:
* wcsgrat.left
* wcsgrat.bottom
* wcsgrat.right
* wcsgrat.top
or (part of) a string which can be (case insensitive)
matched by one from 'left', 'bottom', 'right', 'top'.
:type plotaxis: Integer or String
:param mode: What should this axis do with the tick
marks and labels?
* 0 = ticks native to axis type only
* 1 = only the tick that is not native to axis type
* 2 = both types of ticks (map could be rotated)
* 3 = no ticks
Or use a text that can (case insensitive) match one of:
* "NATIVE_TICKS"
* "SWITCHED_TICKS"
* "ALL_TICKS"
* "NO_TICKS"
:type mode: Integer or String
:param label: An annotation of the current axis
:type label: String
:param `**kwargs`: Keywords for (plot) attributes
:type `**kwargs`: Matplotlib keyword argument(s)
:Examples: Change the font size of the tick labels along
the bottom axis in 11::
grat = Graticule(...)
grat.setp_plotaxis(wcsgrat.bottom, fontsize=11)
"""
#-----------------------------------------------------------------
plotaxis = parseplotaxes(plotaxis)
for ax in plotaxis:
# User wants to make something visible, but right and top
# axis labels are by default invisible. Keyword 'visible' in the
# kwargs list can overrule this default
#self.axes[ax].kwargs.update({'visible':True})
if len(kwargs):
self.axes[ax].kwargs.update(kwargs)
if mode is not None:
mode = parsetickmode(mode)
self.axes[ax].mode = mode
if label is not None:
self.axes[ax].label = label
if xpos is not None:
self.axes[ax].xpos = xpos
if ypos is not None:
self.axes[ax].ypos = ypos
[docs] def setp_axislabel(self, plotaxis=None, label=None, xpos=None, ypos=None, **kwargs):
#-----------------------------------------------------------------
"""
Utility method that calls method :meth:`setp_plotaxis` but
the parameters are restricted to the axis labels.
These labels belong to one of the 4 plot axes.
See the documentation at setp_plotaxis for the input
of the *plotaxis* parameter. The *kwargs* are Matplotlib
attributes.
Possible useful Matplotlib attributes:
* backgroundcolor
* color
* rotation
* style or fontstyle [ 'normal' | 'italic' | 'oblique']
* weight or fontweight
:param plotaxis: The axis number of one of the axes of the plot rectangle:
* wcsgrat.left
* wcsgrat.bottom
* wcsgrat.right
* wcsgrat.top
or (part of) a string which can be (case insensitive)
matched by one from 'left', 'bottom', 'right', 'top'.
:type plotaxis: Integer or String
:param label: The label text.
:type label: String
:param xpos: The x position of the label in normalized device coordinates
:type xpos: Floating point number
:param `**kwargs`: Keywords for (plot) attributes
:type `**kwargs`: Matplotlib keyword argument(s)
"""
#-----------------------------------------------------------------
if plotaxis is None:
plotaxis = [0,1,2,3]
self.setp_plotaxis(plotaxis, mode=None, label=label, xpos=xpos, ypos=ypos, **kwargs)
[docs] def set_tickmode(self, plotaxis=None, mode=None):
#-----------------------------------------------------------------
"""
Utility method that calls method :meth:`setp_plotaxis` but
the parameters are restricted to the tick mode.
Each plot axis has a tick mode.
:param plotaxis: The axis number of one of the axes of the
plot rectangle:
* wcsgrat.left
* wcsgrat.bottom
* wcsgrat.right
* wcsgrat.top
or (part of) a string which can be (minimal & case insensitive)
matched by one from 'left', 'bottom', 'right', 'top'.
:type plotaxis: Integer or String
:param mode: What should this axis do with the tick
marks and labels?
* 0 = ticks native to axis type only
* 1 = only the tick that is not native to axis type
* 2 = both types of ticks (map could be rotated)
* 3 = no ticks
Or use a text that can (minimal) match one of:
* "NATIVE_TICKS"
* "SWITCHED_TICKS"
* "ALL_TICKS"
* "NO_TICKS"
:type mode: Integer or String
"""
#-----------------------------------------------------------------
if plotaxis is None:
plotaxis = [0,1,2,3]
self.setp_plotaxis(plotaxis, mode=mode)
[docs] def Insidelabels(self, wcsaxis=0, world=None, constval=None,
deltapx=0.0, deltapy=0.0, angle=None, addangle=0.0,
fun=None, fmt=None, tex=True, aspect=1.0, **kwargs):
#-----------------------------------------------------------------
"""
Annotate positions in world coordinates
within the boundaries of the plot.
This method can be used to plot positions
on all-sky maps where there are usually no
intersections with the enclosing axes rectangle.
:param wcsaxis: Values are 0 or 1, corresponding to the
first and second world coordinate types.
The accepted values are 0 and 1. The default
is 0.
:type wcsaxis: Integer
:param world: One or a sequence of world coordinates on the axis given
by *wcsaxis*. The positions are completed
with one value for *constval*.
If world=None (the default) then the world
coordinates are copied from graticule
world coordinates.
:type world: One or a sequence of floating point number(s) or None
:param constval: A constant world coordinate to complete the positions
at which a label is plotted. The value can also be a
string representing a sexagesimal number.
:type constval: Floating point number or String
:param deltapx: Small shift in pixels in x-direction of text. This enables
us to improve the layout of the plot by preventing that
labels are intersected by lines.
:type deltapx: Floating point number.
:param deltapy: See description at *deltapx*.
:type deltapy: Floating point number.
:param angle: Use this angle (in degrees) instead of
calculated defaults. It is the angle at which
then **all**
position labels are plotted.
:type angle: Floating point number
:param addangle: Add this angle (in degrees) to the calculated
default angles.
:type addangle: Floating point number
:param fun: Function or lambda expression to convert the
label value.
:type func: Python function or lambda expression
:param fmt: String to format the numbers. If omitted the
format '%g' is used.
:type fmt: String
:param tex: Format these 'inside' labels in LaTeX if this
parameter is set to True (which is the default).
:type param: Boolean
:param aspect: The aspect ratio of the frame. This number is needed to
plot labels at the right angle. It cannot be derived from
the aspect ratio of the frame, because at the moment of
creation, the frame is not known (only after a call to
the plot() method, a frame is known). If the aspect ratio
is known in the calling environment, we should use it
there to get the angles right.
:type aspect: Floating point number
:param `**kwargs`: Keywords for (plot) attributes.
:type `**kwargs`: Matplotlib keyword argument(s)
:returns: An Insidelabel object with a series of derived label objects.
These label objects
have a number of attributes, see :class:`Insidelabels`
:Notes: For a map with only one spatial axis, the value of
'mixpix' is used as pixel value for the
matching spatial axis. The *mixed()* method
from module *wcs* is used to calculate the right
positions.
:Examples: Annotate a plot with labels at positions from a list
with longitudes at given fixed latitude::
grat = Graticule(...)
lon_world = [0,30,60,90,120,150,180]
lat_constval = 30
inlabs = grat.Insidelabels(wcsaxis=0,
world=lon_world,
constval=lat_constval,
color='r')
"""
#-----------------------------------------------------------------
if world is None:
if wcsaxis == 0:
world = self.xstarts
if wcsaxis == 1:
world = self.ystarts
if not issequence(world):
world = [world,]
if constval is None:
if wcsaxis == 0:
constval = self.ystarts[0]
else:
constval = self.xstarts[0]
if isinstance(constval, str):
pos, err = parsehmsdms(constval)
if err != '':
raise Exception(err)
else:
constval = pos[0]
unknown = numpy.nan
wxlim0 = self.wxlim[0]
wxlim1 = self.wxlim[1]
#if self.wxlim[0] < 0.0:
# wxlim0 += 180.0
# wxlim1 += 180.0
insidelabels = Insidelabels(wcsaxis) # Initialize the result
if len(world) > 0 and wcsaxis in [0,1]:
if wcsaxis == 0:
defkwargs = {'ha':'center', 'va':'center', 'fontsize':10}
phi = 0.0
for xw in world:
#if xw < 0.0:
# xw += 180.0
if self.mixpix is None: # Could be projection with matching axis
wt = (xw, constval)
xp, yp = self.gmap.topixel(wt)
else:
wt = (xw, constval, unknown)
pixel = (unknown, unknown, self.mixpix)
(wt, pixel) = self.gmap.mixed(wt, pixel)
xp = pixel[0]; yp = pixel[1]
labval = xw
if xw < 0.0 and self.gmap.types[wcsaxis] == 'longitude':
labval += 360.0
if not numpy.isnan(xp):
#if wxlim0 <= xw < wxlim1 and self.pxlim[0] < xp < self.pxlim[1]:
if self.pxlim[0]-0.5 < xp < self.pxlim[1]+0.5 and self.pylim[0]-0.5 < yp < self.pylim[1]+0.5:
if angle is None:
if self.mixpix is None:
d = (self.wylim[1] - self.wylim[0])/200.0
xp1, yp1 = self.gmap.topixel((xw, constval-d))
xp2, yp2 = self.gmap.topixel((xw, constval+d))
if not (numpy.isnan(xp1) or numpy.isnan(xp2)):
yp1 *= aspect; yp2 *= aspect
phi = numpy.arctan2(yp2-yp1, xp2-xp1)*180.0/numpy.pi
if self.gmap.cdelt[1] < 0.0:
phi -= 180.0
if phi < 0:
phi += 360.0
if 90 <= phi < 270:
phi += 180.0
else:
phi = 90.0
else:
phi = angle
defkwargs.update({'rotation':phi+addangle})
defkwargs.update(kwargs)
insidelabels.append(xp+deltapx, yp+deltapy,
xw, constval,
labval,
phi+addangle,
self.gmap.types[wcsaxis],
skysys=self.__skysys, #self.gmap.skyout,
fun=fun,
fmt=fmt,
offset=False,
prec=self.prec[wcsaxis],
tex=tex,
**defkwargs)
if wcsaxis == 1:
defkwargs = {'ha':'center', 'va':'center', 'fontsize':10}
for yw in world:
phi = 0.0
if self.mixpix is None:
wt = (constval, yw)
xp, yp = self.gmap.topixel(wt)
else:
wt = (constval, yw, unknown)
pixel = (unknown, unknown, self.mixpix)
(wt, pixel) = self.gmap.mixed(wt, pixel)
xp = pixel[0]; yp = pixel[1]
labval = yw
if yw < 0.0 and self.gmap.types[wcsaxis] == 'longitude':
labval += 360.0
if not numpy.isnan(xp):
if self.wylim[0] <= yw < self.wylim[1] and self.pylim[0] < yp < self.pylim[1] and self.pxlim[0] < xp < self.pxlim[1]:
# Delta's make minus sign more visible on graticule lines
if angle is None:
if self.mixpix is None:
d = (self.wxlim[1] - self.wxlim[0])/200.0
xp1, yp1 = self.gmap.topixel((constval-d, yw))
xp2, yp2 = self.gmap.topixel((constval+d, yw))
if not (numpy.isnan(xp1) or numpy.isnan(xp2)):
yp1 *= aspect; yp2 *= aspect
phi = numpy.arctan2(yp2-yp1, xp2-xp1)*180.0/numpy.pi
if self.gmap.cdelt[0] < 0.0:
phi -= 180.0
if phi < 0:
phi += 360.0
if 90 <= phi < 270:
phi += 180.0
else:
phi = angle
defkwargs.update({'rotation':phi+addangle})
defkwargs.update(kwargs)
insidelabels.append(xp+deltapx, yp+deltapy,
constval, yw,
labval,
phi+addangle,
self.gmap.types[wcsaxis],
skysys=self.__skysys, #self.gmap.skyout,
fun=fun,
fmt=fmt,
offset=False,
prec=self.prec[wcsaxis],
tex=tex,
**defkwargs)
insidelabels.pxlim = self.pxlim
insidelabels.pylim = self.pylim
self.objlist.append(insidelabels) # This graticule stores a list of derived objects
return insidelabels
def Ruler(self, pos1=None, pos2=None,
x1=None, y1=None, x2=None, y2=None, lambda0=0.5, step=None,
world=False, angle=None, addangle=0.0,
fmt=None, fun=None, fliplabelside=False, mscale=None,
labelsintex=True, **kwargs):
#-----------------------------------------------------------------
"""
Look at documentation of same method of class Annotatedimage.
This is a version that is needed to calculate offsets along
the plot axes.
"""
#-----------------------------------------------------------------
ruler = rulers.Ruler(self.gmap, self.mixpix, self.pxlim, self.pylim,
pos1=pos1, pos2=pos2,
x1=x1, y1=y1, x2=x2, y2=y2, lambda0=lambda0, step=step,
world=world, angle=angle, addangle=addangle,
fmt=fmt, fun=fun, fliplabelside=fliplabelside, mscale=mscale,
labelsintex=labelsintex, **kwargs)
return ruler
# ==========================================================================
# function _update_metrics()
# --------------------------------------------------------------------------
# The function _update_metrics() is a replacement for the method
# matplotlib.mathtext.Char._update_metrics(). It tries to modify
# the metrics of characters 'm' and 's' so that hms tick labels
# are aligned correctly. To prevent unwanted side effects in other
# text labels, this is done only when the font is 'rm' and the current
# fontsize is less than variable 'tweakhms' in an attempt to
# limit the metrics modification to minute and second superscripts.
#
def _update_metrics(self):
global tweakhms
metrics = self._metrics = self.font_output.get_metrics(
self.font, self.font_class, self.c, self.fontsize, self.dpi)
if self.c in ['m', 's'] and \
self.fontsize<=float(tweakhms) and \
self.font=='rm':
metrics_ms = self.font_output.get_metrics(
self.font, self.font_class, 'h', self.fontsize, self.dpi)
metrics.iceberg = self._metrics.iceberg = metrics_ms.iceberg
metrics.height = self._metrics.height = metrics_ms.height
if self.c == ' ':
self.width = metrics.advance
else:
self.width = metrics.width
self.height = metrics.iceberg
self.depth = -(metrics.iceberg - metrics.height)
# from matplotlib.mathtext import Char
# VOG: Char._update_metrics = _update_metrics