Tag Archives: python

Converting Between Decimal Degrees and Hours, Minutes, Seconds

Posted on by
9

Here’s a quick Python snippet I wrote to convert right ascension in decimal degrees to hours, minutes, seconds and declination to (+/-)degrees, minutes, seconds.

def deg2HMS(ra='', dec='', round=False):
  RA, DEC, rs, ds = '', '', '', ''
  if dec:
    if str(dec)[0] == '-':
      ds, dec = '-', abs(dec)
    deg = int(dec)
    decM = abs(int((dec-deg)*60))
    if round:
      decS = int((abs((dec-deg)*60)-decM)*60)
    else:
      decS = (abs((dec-deg)*60)-decM)*60
    DEC = '{0}{1} {2} {3}'.format(ds, deg, decM, decS)
  
  if ra:
    if str(ra)[0] == '-':
      rs, ra = '-', abs(ra)
    raH = int(ra/15)
    raM = int(((ra/15)-raH)*60)
    if round:
      raS = int(((((ra/15)-raH)*60)-raM)*60)
    else:
      raS = ((((ra/15)-raH)*60)-raM)*60
    RA = '{0}{1} {2} {3}'.format(rs, raH, raM, raS)
  
  if ra and dec:
    return (RA, DEC)
  else:
    return RA or DEC

For example:
In [1]: f.deg2HMS(ra=66.918277)
Out[1]: '4 27 40.386'

Or even:
In [2]: f.deg2HMS(dec=24.622590)
Out[2]: '+24 37 21.324'

Or if you want to round the seconds, just do:
In [3]: f.deg2HMS(dec=24.622590,round=True)
Out[3]: '+24 37 21'

And to convert hours, minutes and seconds into decimal degrees, we have:

def HMS2deg(ra='', dec=''):
  RA, DEC, rs, ds = '', '', 1, 1
  if dec:
    D, M, S = [float(i) for i in dec.split()]
    if str(D)[0] == '-':
      ds, D = -1, abs(D)
    deg = D + (M/60) + (S/3600)
    DEC = '{0}'.format(deg*ds)
  
  if ra:
    H, M, S = [float(i) for i in ra.split()]
    if str(H)[0] == '-':
      rs, H = -1, abs(H)
    deg = (H*15) + (M/4) + (S/240)
    RA = '{0}'.format(deg*rs)
  
  if ra and dec:
    return (RA, DEC)
  else:
    return RA or DEC

So we can get back our decimal degrees with:

In [4]: f.HMS2deg(ra='4 27 40.386', dec='+24 37 21.324')
Out[4]: (66.918, 24.622)

How to Make Animated GIFs of Plots

Posted on by
4

So you’re swimming in countless plots of your data over some changing parameter. Sure they’re nice to look at, but are they animated? Didn’t think so.

Here’s how to create an animated .gif image of those Python plots using Photoshop.

Step 1: Generate the plots

I’ve found that a few lines of Python to programmatically draw and save your plots to a folder eliminates a lot of editing and tweaking later on.

For this tutorial, I’ll use my synthetic photometry code to generate color-parameter plots across surface gravities of 3.0dex to 6.0dex.

First I created a folder on my desktop to dump the plots called RI_teff.

Then in the Python interpreter it looks like:

In [1]: grav = [round(i*0.1,1) for i in range(30,61)]
In [2]: import syn_phot as s
In [3]: for i in g:
s.color_param('R','I',logg=i,save='/Desktop/RI_teff/')

Now I can create the animated .gif in Photoshop.

Step 2: Pull plots into Photoshop as layers

Open Photoshop and click File > Scripts > Load Files into Stack…

Select the folder on your Desktop that has all the plots and click ok.

Photoshop will open each image as a layer in a new project window.

Step 3: Create frames from each layer

Next click Window > Timeline to show the Timeline across the bottom of the program.

In the top-right corner of your Timeline, you’ll see a little button that has all the Timeline options. Click it and select Make Frames From Layers. Here’s what it looks like:

This will populate your Timeline with one frame for each image layer.

Click the Timeline options button again and click Reverse Frames if necessary. Otherwise, you can drag and drop the frames into the desired order.

Step 4: Timing is everything

Next we need to set the timing of each frame. Select the first frame from the Timeline then hold down the Shift button and select the last frame to select them all.

Next click on any frame where it says 0 sec. with a down arrow. Then select the display time for each frame in the animation.

I typically set the frames to be 0.2 or 0.5 seconds each, depending on the number of total frames. Then I set the last frame to be 2 seconds so it’s as if the image pauses when it finishes before starting the animation over.

Step 5: Save it!

Finally, click File > Save for Web… and make sure you have GIF filetype selected. Click Save and you’re done! Here’s the result:

Editing PYTHONPATH (or “Where’s my module?!”)

Posted on by
11

Python not finding your modules with import MyModule for some reason? It’s probably your PYTHONPATH. Here’s how you edit it so that Python sees all your modules.

A Word About Directories

Firstly, if you are the kind of person who keeps files scattered about your computer, dumps everything into the root directory, or is afraid to use those nice little blue folders, then perhaps programming and computers in general are not for you.

Logical and neat directory structure will make your own, your computer’s and your collaborators’ lives much easier.

My recommendation: Create a master code directory (called “Modules” or something) in your Documents folder. This will be the new home of all your future Python projects.

Now every time you create a .py file, it should go into a project folder inside your Modules directory, i.e. /Documents/Modules/NewProject/MyModule.py. Note that you should have no actual modules inside your Modules directory! Keep those puppies in their warm, snuggly project subdirectories.

This way you can also initialize that project folder (i.e. /NewProject) as a Git repository and just push and pull the contents to keep it up-to-date!

Editing PYTHONPATH

Python won’t just search your computer for the MyModule.py file you’re trying to import. You have to tell it explicitly each time where to get it. The PYTHONPATH is a list of directories for your computer to check whenever you type import MyModule into the interpreter.

To add a path, launch ipython and type:

import sys
sys.path.append("path/to/Modules")
print sys.path

Note: You only have to update your PYTHONPATH once, not every time you use Python!

So now your path is updated but this is only the path to your master Python folder.

In order to have Python see the modules inside each subdirectory, add a blank file called __init__.py to each subdirectory (with two underscores on each side).

Now to import the module and use a function called foo() do:

from NewProject import MyModule as m
m.foo()

That is, it’s checking the main python directory you added to your PYTHONPATH, then looking within the NewProject subdirectory via the __init__.py file for the module you’re trying to import.

Brown Dwarf Synthetic Photometry

Posted on by
1

The goal here was to get the synthetic colors in the SDSS, 2MASS and WISE filters of ~2000 model objects generated by the PHOENIX stellar and planetary atmosphere software.

Since it would be silly (and incredibly slow… and much more boring) to just calculate and store every single color for all 12 filter profiles, I wrote a module to calculate colors a la carte.

The Filters

I got the J, H, and K band relative spectral response (RSR) curves in the 2MASS documentation, the u, g, r, i and z bands from the SDSS documentation, and the W1, W2, W3, and W4 bands from the WISE documentation.

I dumped all my .txt filter files into one directory and wrote a function to grab them all, pull out the wavelength and transmission values, and output the filter name in position [0], x-values in [1], and y-values in [2]:

def get_filters(filter_directory):
  import glob, os
  files = glob.glob(filter_directory+'*.txt')
 
  if len(files) == 0:
    print 'No filters in', filter_directory
  else:
    filter_names = [os.path.splitext(os.path.basename(i))[0] for i in files]
    RSR = [open(i) for i in files]
    filt_data = [filter(None,[map(float,i.split()) for i in j if not i.startswith('#')]) for j in RSR]
    for i in RSR: i.close()
 
    RSR_x = [[x[0] for x in i] for i in filt_data]
    RSR_y = [[y[1] for y in i] for i in filt_data]
    filters = {}
    for i,j,k in zip(filter_names,RSR_x,RSR_y):
      filters[i] = j, k, center(i)
 
    return filters

Calculating Apparent Magnitudes

We can’t have colors without magnitudes so here’s a function to grab the Teff and log g specified spectra, and calculate the apparent magnitudes in a particular band:

def mags(band, teff='', logg='', bin=1):
  from scipy.io.idl import readsav
  from collections import Counter
  from scipy import trapz, log10, interp
  
  s = readsav(path+'modelspeclowresdustywise.save')
  Fr, Wr = [i for i in s.modelspec['fsyn']], [i for i in s['wsyn']]
  Tr, Gr = [int(i) for i in s.modelspec['teff']], [round(i,1) for i in s.modelspec['logg']]
  
  # The band to compute
  RSR_x, RSR_y, lambda_eff = get_filters(path)[band]
  
  # Option to specify an effective temperature value
  if teff:
    t = [i for i, x in enumerate(s.modelspec['teff']) if x == teff]
    if len(t) == 0:
      print "No such effective temperature! Please choose from 1400K to 4500K in 50K increments or leave blank to select all."
  else:
    t = range(len(s.modelspec['teff']))
  
  # Option to specify a surfave gravity value
  if logg:
    g = [i for i, x in enumerate(s.modelspec['logg']) if x == logg]
    if len(g) == 0:
      print "No such surface gravity! Please choose from 3.0 to 6.0 in 0.1 increments or leave blank to select all."
  else:
    g = range(len(s.modelspec['logg']))
  
  # Pulls out objects that fit criteria above
  obj = list((Counter(t) & Counter(g)).elements())
  F = [Fr[i][::bin] for i in obj]
  T = [Tr[i] for i in obj]
  G = [Gr[i] for i in obj]
  W = Wr[::bin]
  
  # Interpolate to find new filter y-values
  I = interp(W,RSR_x,RSR_y,left=0,right=0)
  
  # Convolve the interpolated flux with each filter (FxR = RxF)
  FxR = [convolution(i,I) for i in F]
  
  # Integral of RSR curve over all lambda
  R0 = trapz(I,x=W)
  
  # Integrate to find the spectral flux density per unit wavelength [ergs][s-1][cm-2] then divide by R0 to get [erg][s-1][cm-2][cm-1]
  F_lambda = [trapz(y,x=W)/R0 for y in FxR]
  
  # Calculate apparent magnitude of each spectrum in each filter band
  Mags = [round(-2.5*log10(m/F_lambda_0(band)),3) for m in F_lambda]
  
  result = sorted(zip(Mags, T, G, F, I, FxR), key=itemgetter(1,2))
  result.insert(0,W)
  
  return result

Calculating Colors

Now we can calculate the colors. Next, I wrote a function to accept any two bands with options to specify a surface gravity and/or effective temperature as well as a bin size to cut down on computation. Here’s the code:

def colors(first, second, teff='', logg='', bin=1):
  (Mags_a, T, G) = [[i[j] for i in get_mags(first, teff=teff, logg=logg, bin=bin)[1:]] for j in range(3)]
  Mags_b = [i[0] for i in get_mags(second, teff=teff, logg=logg, bin=bin)[1:]]
  colors = [round(a-b,3) for a,b in zip(Mags_a,Mags_b)]
  
  print_mags(first, colors, T, G, second=second)
  
  return [colors, T, G]

The PHOENIX code gives the flux as Fλ in cgs units [erg][s-1][cm-2][cm-1] but as long as both spectra are in the same units the colors will be the same.

Makin’ It Handsome

Then I wrote a short function to print out the magnitudes or colors in the Terminal:

def print_mags(first, Mags, T, G, second=''):
  LAYOUT = "{!s:10} {!s:10} {!s:25}"
  
  if second:
    print LAYOUT.format("Teff", "log g", first+'-'+second)
  else:
    print LAYOUT.format("Teff", "log g", first)
  
  for i,j,k in sorted(zip(T, G, Mags)):
    print LAYOUT.format(i, j, k)

The Output

Then if I just want the J-K color for objects with log g = 4.0 over the entire range of effective temperatures, I launch ipython and just do:

In [1]: import syn_phot as s
In [2]: s.colors('J','K', logg=4)
Teff -------- log g -------- J-K
1400.0 ------ 4.0 ---------- 4.386
1450.0 ------ 4.0 ---------- 4.154
...
4450.0 ------ 4.0 ---------- 0.756
4500.0 ------ 4.0 ---------- 0.733

Similarly, I can specify just the target effective temperature and get the whole range of surface gravities. Or I can specify an effective temperature AND a specific gravity to get the color of just that one object with:

In [3]: s.colors('i','W2', teff=3050, logg=5)
Teff -------- log g -------- J-K
3050.0 ------ 5.0 ---------- 3.442

I can also reduce the number of data points in each flux array if my sample is very large. I just have to specify the number of data points to skip with the “bin” optional parameter. For example:

In [4]: s.colors('W1','W2', teff=1850, bin=3)

This will calculate the W1-W2 color for all the objects with Teff = 1850K and all gravities, but only take every third flux value.

I also wrote functions to generate color-color, color-parameter and color-magnitude plots but those will be in a different post.

Plots!

Here are a few color-parameter animated plots I made using my code. Here’s how I made them. Click to animate!

And here are a few colorful-colorful color-color plots I made:

Plots with observational data

Just to be sure I’m on base, here’s a color-color plot of J-H vs. H-Ks for objects with a log surface gravity of 5 dex (blue dots) plotted over some data for the Chamaeleon I Molecular Cloud (semi-transparent) from Carpenter et al. (2002).

The color scale is for main sequence stars and the black dots are probable members of the group. Cooler dwarfs move up and to the right.

And here’s a plot of J-Ks vs. z-Ks as well as J-Ks vs. z-J. Again, the blue dots are from my synthetic photometry code at log(g)=5 and the semi-transparent points with errors are from Dahn et al. (2002).

Editing astrotools

Posted on by
Reply

It is important to maintain some standards in astrotools. The goal is to keep it readable and stable. With than in mind, I shall provide some coding guidelines and then describe the structure of astrotools in detail to allow for its easy editing.

Actually, I will let Guido van Rossum (the father of Python) provide you with all the necessary coding guidelines, from naming conventions to code lay-out. You can find them in Guido’s Style Guide for Python Code.

Let me clear out an important aspect about Python modules. A module is a collection of reusable pieces of code, such as functions and classes. A piece of code in a module can use another piece within the same module or in another module, or it can be completely independent. The implications are that you can add a function:

  • that is completely self-contained and that does not interfere with anything else in astrotools
  • that makes use of other astrotools functions/classes
  • that makes use of other Python modules (e.g. numpy), in which case you need a command line to invoke such module (more on this below).

The order of things in astrotools

The module astrotools is broken down into five sections. From top to bottom in the code, they are described below.

I) Documentation: Contains general information about the contents of the module. This type of non-executable text (a.k.a. docstring) must always start and end with three apostrophes (which I like to call Melchior, Caspar, & Balthazar). This tells Python that the stuff within them is a string that can be called by the end user using Python help.

II) External Modules: Contains the Python commands used to invoke external Python modules —the famous “import” commands. Having them all in one location as opposed to inside the functions/classes where they are needed helps to avoid repetition (when more than one function/class needs to invoke the same external module) and to clearly see the external modules needed to use astrotools.

III) Public Functions: Organized alphabetically, these functions are meant to be accessible to end users. These functions have their own docstrings, surrounded by the three apostrophes, just like in Section I.

IV) Non-Public functions: Organized alphabetically, these functions are only used by other functions/classes within astrotools. As such, end users have no business accessing these. Python convention says non-public names should be preceded by double underscore. Technically, anyone can still access them, but Python code of conduct tells end users not to do it. These functions do not need docstrings, only commented (#) intro lines to explain their purpose and where they are used.

V) Public Classes: Organized alphabetically, these classes are meant to be accessible to end users. Just like with Public Functions, Public Classes have their own docstrings.

How to add a function/class to astrotools

  1. Download the latest version of astrotools from the GitHub repository into your computer.
  2. If your function/class needs to import an external module, check if it is already imported in Section II of astrotools. If not, add the command line to do so. Take as much advantage as you can of the modules already imported in astrotools before adding a new one. More external modules needed by astrotools means more requirements that end users will have to meet before having a fully functional astrotools module.
  3. Add your function/class in the appropriate location (alphabetically) within the relevant section (III or V) of astrotools.
  4. Your function/class MUST have a docstring, and it should follow reStructuredText syntax. No panic! It is very simple. For instance, if you want a word in italics, you write: *word*. Consult the reStructuredText User Reference for more information. In terms of format and contents, follow the example of docstrings of existing astrotools functions. Include in the docstring a description of all input parameters in your function/class, including their formats. Remember, the point of documenting your code is to allow others to easily use it and edit it in case you die. This is why I like to call the docstring “the death string”.
  5. If your function/class comes with a sub-function, that is, a side function that you wrote that is used by your main function/class, this means you are a good (object-oriented) programmer! The place to put it is Section IV of astrotools. Include in it intro comments about its purpose and where it is used. If you think that this sub-function could also be useful to end users, by all means add it in Section III instead. Just remember to write a docstring for it.
  6. Test the function/class that you added. Test extreme cases, test wrong input formats, test, test, test!
  7. Update astrotools.py in the GitHub repository with your new edited version.
  8. Inform the person that maintains the public documentation for astrotools so that it gets updated in bdnyc.org/astrotools.

Regularizing Data

Posted on by
2

I have been trying to cross correlate two spectra with the intent of finding the radial velocity of one relative to the other.  I’ve been working on a Python code to do this for some time, but was having complications.  Although there is clearly a shift between the two spectra, the numpy cross correlation function that the code uses kept coming up with a shift of zero.  I recently came across the concept of “regularizing” data on a python help website.  Someone was having a similar problem (coming up with a cross correlation result of a zero shift), and it was pointed out that they weren’t regularizing their data.  To regularize data is to subtract off the mean, and divide by the standard deviation.  I tried doing this in my code, and it actually started to work, and give me very nice answers!  The problem is that I don’t particularly understand why this makes the cross correlation work.

Is anyone familiar with the concept of regularizing data? Or has anyone had a similar problem regarding cross correlations?

Fitting a gaussian to your data

Posted on by
4

Let’s say you have a nice histogram, like this…

…and you want to fit a gaussian to it so that you can find the mean, and the standard deviation.  Follow these steps!

First, we have to make sure we have the right modules imported

>>> import matplotlib.pyplot as plt
>>> import matplotlib.mlab as mlab
>>> from scipy.stats import norm

Let’s say your data is stored in some array called data.

>>> (mu,sigma) = norm.fit(data)
Mu is the mean, and sigma is one standard deviation. If you don’t care about plotting your data, you can stop here.
>>> plt.figure(1)
>>> n,bins,patches=plt.hist(data,20,normed=1,facecolor=’green’,align=’mid’)
The number after data (20) is the number of bins you want your data to go into. Normed has to do with the integral of the gaussian.
>>> y = mlab.normpdf(bins,mu,sigma)
>>> plt.plot(bins,y,’r–‘,linewidth=2)

Now your data is nicely plotted as a histogram and its corresponding gaussian!

Intro to astrotools

Posted on by
Reply

The BDNYC Python module astrotools is now up and running for all to use and edit. This module is intended to be a collection of useful pieces of code for our everyday data handling needs. So far, astrotools is useful to handle fits files and spectral data. Additions to astrotools are welcome.

Let me clear out some relevant Python concepts. A module is simply a set of functional pieces of code. These pieces can be either functions or classes, among other things. A function is a code that executes an action, for example “make me a sandwich”. A class is a code that creates an object with specific characteristics and it allows you to handle this object very effectively. For the example above, you can have a class to create the person that will make the sandwich, a class to create a slice of mortadella, and a class to create mayonnaise. Of course, you can have functions that make no use of classes. More on that in a future post.

Where to find astrotools and how to install it

The module astrotools is available at https://github.com/BDNYC/astrotools as astrotools.py. Make sure that you download it into one of your Python folders in sys.path, which are the folders where Python looks for modules whenever you invoke one in your Python session. If you have no idea what I am referring to, please read this piece of Python documentation.

How to use astrotools

Now you are ready to invoke astrotools into your Python session. Type:

“import astrotools”

That’s it! You are ready to use astrootols functions and classes. For example, if you want to use the function that gets spectral data from a fits file, type:

“x = astrotools.read_spec(filename)”

and x will hold the output of the read_spec function. To learn about the functions and classes in astrotools, visit https://bdnyc.org/astrotools.

By convention, modules are invoked using a nickname, if only to make typing easier. I
suggest the nickname “at”. Then, import astrotools by typing:

“import astrotools as at”

To use it, type:

“x = at.read_spec(filename)”

In later posts I will describe in greater detail some important Python concepts as well as some guidelines to follow when editing astrotools.