Besancon Model of the Galaxy

astrodata/besancon/__init__.py

@@ -0,0 +1 @@
+from .besancon import *

astrodata/besancon/besancon.py

@@ -0,0 +1,284 @@
+"""
+Besancon Query Tool
+-------------------
+A tool to query the Besancon model of the galaxy
+http://model.obs-besancon.fr/
+
+:Author: Adam Ginsburg (adam.g.ginsburg@gmail.com)
+"""
+import urllib
+import urllib2
+import socket
+import time
+import copy
+from astrodata.utils import progressbar
+import sys
+import re
+
+__all__ = ['get_besancon_model_file','request_besancon']
+
+keyword_defaults = {
+    'rinf':0.000000,
+    'rsup':50.000000,
+    'dist_step_mode':0,
+    'dlr':  0.000,
+    'kleg':1,
+    'longit': 10.62,
+    'latit':-0.38,
+    'soli':0.0003, # degrees.  0.00027777 = 1 arcmin
+    'kleh':1,
+    'eq1': 2000.0,
+    'al0': 200.00,
+    'alm': 200.00,
+    'dl':  1.00,
+    'ab0':  59.00,
+    'abm':  59.00,
+    'db': 1.00,
+    'adif': 0.700,
+    'ev':[""]*24,
+    'di':[""]*24,
+    'oo':[-7]+[-99]*12,
+    'ff':[15]+[99]*12,
+    'spectyp_min':1,
+    'subspectyp_min': 0,
+    'spectyp_max':9,
+    'subspectyp_max': 5,
+    'lumi[]':range(1,8),
+    'sous_pop[]':range(1,11),
+    'iband':8,
+    'band0':[8]*9,
+    'bandf':[25]*9,
+    'colind':["J-H","H-K","J-K","V-K",],
+    'nic':  4,
+    'klea':1,
+    'sc':[[0,0,0]]*9,
+    'klee':0,
+    'throughform':'ok',
+    'kleb':3,
+    'klec':1,
+    'cinem':0,
+    'outmod':"",
+    }
+
+url_download = "ftp://sasftp.obs-besancon.fr/modele/"
+url_request  = "http://model.obs-besancon.fr/modele_form.php"
+# sample file:  1340900648.230224.resu
+result_re = re.compile("[0-9]{10}\.[0-9]{6}\.resu")
+
+def parse_besancon_dict(bd):
+    """
+    Turn a dict like default_keys into a list of tuples (must be a list of
+    tuples because there are some repeated entries, which dictionaries do not
+    support)
+    """
+
+    http_dict = []
+    for key,val in bd.iteritems():
+        if type(val) is list:
+            if "[]" in key:
+                for listval in val:
+                    http_dict.append((key,listval))
+            else:
+                for ii,listval in enumerate(val):
+                    if type(listval) is list:
+                        for jj,lv in enumerate(listval):
+                            http_dict.append((key+"[%i][%i]" % (ii,jj),lv))
+                    else:
+                        http_dict.append((key+"[%i]" % (ii) , listval))
+        else:
+            http_dict.append((key , val))
+
+    return http_dict
+
+def parse_errors(text):
+    """
+    Attempt to extract the errors from a Besancon web page with error messages in it
+    """
+    try:
+        errors = re.compile(r"""<div\ class="?errorpar"?>\s*
+                        <ol>\s*
+                        (<li>([a-zA-Z0-9):( \s_-]*)</li>\s*)*\s*
+                        </ol>\s*
+                        </div>""", re.X)
+        text = errors.search(text).group()
+    except AttributeError:
+        likely_errors = text.split('\n')[132:150]
+        raise ValueError("Regular expression matching to error message failed.")
+    text_items = re.split("<li>|</li>|\n",errors.search(text).group())
+    text_items = [t for t in text_items if t != ""]
+    error_list = text_items[2:-2]
+    return error_list
+
+
+colors_limits = {"J-H":(-99,99),"H-K":(-99,99),"J-K":(-99,99),"V-K":(-99,99)}
+mag_limits = {'U':(-99,99), 'B':(-99,99), 'V':(-5,20), 'R':(-99,99),
+    'I':(-99,99), 'J':(-99,99), 'H':(-99,99), 'K':(-99,99), 'L':(-99,99)}
+mag_order = "U","B","V","R","I","J","H","K","L"
+
+def request_besancon(email, glon, glat, smallfield=True, extinction=0.7,
+        area=0.0001, verbose=True, clouds=None, absmag_limits=(-7,15),
+        mag_limits=copy.copy(mag_limits),
+        colors_limits=copy.copy(colors_limits), 
+        retrieve_file=True, **kwargs):
+    """
+    Perform a query on the Besancon model of the galaxy
+    http://model.obs-besancon.fr/
+
+    Parameters
+    ----------
+    email : string
+        A valid e-mail address to send the report of completion to
+    glon : float
+    glat : float
+        Galactic latitude and longitude at the center
+    smallfield : bool
+        Small field (True) or Large Field (False)
+        LARGE FIELD NOT SUPPORTED YET
+    extinction : float
+        Extinction per kpc in A_V
+    area : float
+        Area in square degrees 
+    absmag_limits : (float,float)
+        Absolute magnitude lower,upper limits
+    colors_limits : dict of (float,float)
+        Should contain 4 elements listing color differences in the valid bands, e.g.:
+            {"J-H":(99,-99),"H-K":(99,-99),"J-K":(99,-99),"V-K":(99,-99)}
+    mag_limits = dict of (float,float) 
+        Lower and Upper magnitude difference limits for each magnitude band
+        U B V R I J H K L
+    clouds : list of 2-tuples
+        Up to 25 line-of-sight clouds can be specified in pairs of (A_V,
+        distance in pc)
+    verbose : bool
+        Print out extra error messages?
+    retrieve_file : bool
+        If True, will try to retrieve the file every 30s until it shows up.
+        Otherwise, just returns the filename (the job is still executed on
+        the remote server, though)
+    kwargs : dict
+        Can override any argument in the request if you know the name of the
+        POST keyword.
+
+    """
+
+    # create a new keyword dict based on inputs + defaults
+    kwd = copy.copy(keyword_defaults)
+    for key,val in kwargs.iteritems():
+        if key in keyword_defaults:
+            kwd[key] = val
+        elif verbose:
+            print "Skipped invalid key %s" % key
+
+    kwd['kleg'] = 1 if smallfield else 2
+    if not smallfield:
+        raise NotImplementedError
+
+    kwd['adif'] = extinction
+    kwd['soli'] = area
+    kwd['oo'][0] = absmag_limits[0]
+    kwd['ff'][0] = absmag_limits[1]
+
+    for ii,(key,val) in enumerate(colors_limits.items()):
+        if key[0] in mag_order and key[1] == '-' and key[2] in mag_order:
+            kwd['colind'][ii] = key
+            kwd['oo'][ii+9] = val[0]
+            kwd['ff'][ii+9] = val[1]
+        else:
+            raise ValueError('Invalid color %s' % key)
+
+    for (key,val) in mag_limits.iteritems():
+        if key in mag_order:
+            kwd['band0'][mag_order.index(key)] = val[0]
+            kwd['bandf'][mag_order.index(key)] = val[1]
+        else:
+            raise ValueError('Invalid band %s' % key)
+
+    if clouds is not None:
+        for ii,(AV,di) in enumerate(clouds):
+            kwd[ev][ii] = AV
+            kwd[di][ii] = di
+
+    # parse the default dictionary
+    request = parse_besancon_dict(keyword_defaults)
+
+    # an e-mail address is required
+    request.append(('email',email))
+    request = urllib.urlencode(request)
+    # load the URL as text
+    U = urllib.urlopen(url_request, request)
+    # keep the text stored for possible later use
+    text = U.read()
+    try:
+        filename = result_re.search(text).group()
+    except AttributeError: # if there are no matches
+        errors = parse_errors(text)
+        raise ValueError("Errors: "+"\n".join(errors))
+
+    if verbose:
+        print "File is %s" % filename
+
+    if retrieve_file:
+        return get_besancon_model_file(filename)
+    else:
+        return filename
+
+def get_besancon_model_file(filename, verbose=True, save=True, savename=None, overwrite=True):
+    """
+    Download a Besancon model from the website
+
+    Parameters
+    ----------
+    filename : string
+        The besancon filename, with format ##########.######.resu
+    verbose : bool
+        Print details about the download process
+    save : bool
+        Save the table after acquiring it?
+    savename : None or string
+        If not specified, defaults to the .resu table name
+    overwrite : bool
+        Overwrite the file if it exists?  Defaults to True because the .resu
+        tables should have unique names by default, so there's little risk of
+        accidentally overwriting important information
+    """
+
+    url = url_download+filename
+
+    elapsed_time = 0
+    t0 = time.time()
+
+    sys.stdout.write("\n")
+    while 1:
+        sys.stdout.write(u"\r")
+        try:
+            U = urllib2.urlopen(url,timeout=5)
+            if verbose:
+                print ""
+                print "Loading page..."
+                results = progressbar.chunk_read(U, report_hook=progressbar.chunk_report)
+            else:
+                results = page.read()
+            break
+        except urllib2.URLError:
+            sys.stdout.write(u"Waiting 30s for model to finish (elapsed wait time %is, total %i)\r" % (elapsed_time,time.time()-t0))
+            time.sleep(30)
+            elapsed_time += 30
+            continue
+        except socket.timeout:
+            sys.stdout.write(u"Waiting 30s for model to finish (elapsed wait time %is, total %i)\r" % (elapsed_time,time.time()-t0))
+            time.sleep(30)
+            elapsed_time += 30
+            continue
+
+
+    if save:
+        if savename is None:
+            savename = filename
+        if not overwrite and os.path.exists(savename):
+            raise IOError("File %s already exists." % savename)
+        outf = open(savename,'w')
+        print >>outf,results
+        outf.close()
+
+    return results
+

astrodata/utils/progressbar.py

@@ -0,0 +1,48 @@
+import urllib2, sys
+
+__all__ = ['chunk_report','chunk_read']
+
+def chunk_report(bytes_so_far, chunk_size, total_size):
+    if total_size > 0:
+        percent = float(bytes_so_far) / total_size
+        percent = round(percent*100, 2)
+        sys.stdout.write(u"Downloaded %12d of %12d bytes (%6.2f%%)\r" % 
+            (bytes_so_far, total_size, percent))
+    else:
+        sys.stdout.write(u"Downloaded %10.2g Mb\r" % 
+            (bytes_so_far / 1024.**2))
+
+
+def chunk_read(response, chunk_size=1024, report_hook=None):
+    content_length = response.info().getheader('Content-Length')
+    if content_length is None:
+        total_size = 0
+    else:
+        total_size = content_length.strip()
+        total_size = int(total_size)
+
+    bytes_so_far = 0
+
+    result_string = ""
+
+    #sys.stdout.write("Beginning download.\n")
+
+    while 1:
+       chunk = response.read(chunk_size)
+       result_string += chunk
+       bytes_so_far += len(chunk)
+
+       if not chunk:
+           if report_hook: 
+               sys.stdout.write('\n')
+           break
+
+       if report_hook:
+          report_hook(bytes_so_far, chunk_size, total_size)
+
+    return result_string
+
+if __name__ == '__main__':
+   response = urllib2.urlopen('http://www.ebay.com')
+   C = chunk_read(response, report_hook=chunk_report)
+