Query Argonaut Server (for 3D Dust Map)

query_argonaut.pro

;+
; NAME:
;   query_argonaut
;
; PURPOSE:
;   Query the Argonaut server for 3D dust information or SFD
;
; CALLING SEQUENCE:
;   qresult = query_argonaut(/struct, /debug, _extra=coords)
;
; INPUTS:
;   ra, dec   : numeric scalars or arrays [deg]
;     OR
;   l, b      : numeric scalars or arrays [deg]
;   mode      : 'full', 'lite', or 'sfd'.  Default to 'full'
;   structure : set this keyword to return structure instead of hash
;   debug     : set to return timing information
;   
; OUTPUTS:
;   qresult   : a hash (or structure, if /structure set) containing
;  'distmod':    The distance moduli that define the distance bins.
;  'best':       The best-fit (maximum probability density)
;                  line-of-sight reddening, in units of SFD-equivalent
;                  E(B-V), to each distance modulus in 'distmod.' See
;                  Schlafly & Finkbeiner (2011) for a definition of the
;                  reddening vector (use R_V = 3.1).
;  'samples':    Samples of the line-of-sight reddening, drawn from
;                  the probability density on reddening profiles.
;  'success':    1 if the query succeeded, and 0 otherwise.
;  'converged':  1 if the line-of-sight reddening fit converged, and
;                  0 otherwise.
;  'n_stars':    # of stars used to fit the line-of-sight reddening.
;  'DM_reliable_min':  Minimum reliable distance modulus in pixel.
;  'DM_reliable_max':  Maximum reliable distance modulus in pixel.
;
; EXAMPLES:
;   IDL> qresult = query_argonaut(l=90, b=10)
;   IDL> help,qresult
;   QRESULT         HASH  <ID=1685  NELEMENTS=13>
;   IDL> qresult = query_argonaut(l=90, b=10, /struct)
;   IDL> help,qresult
;   ** Structure <d76608>, 13 tags, length=5776, data length=5776, refs=1:
;      GR              DOUBLE    Array[31]
;      SUCCESS         LONG64                         1
;      N_STARS         LONG64                       750
;      DEC             DOUBLE           54.568690
;      DM_RELIABLE_MAX DOUBLE           15.189000
;      CONVERGED       LONG64                         1
;      DISTMOD         DOUBLE    Array[31]
;      L               DOUBLE           90.000000
;      B               DOUBLE           10.000000
;      RA              DOUBLE          -54.585789
;      BEST            DOUBLE    Array[31]
;      DM_RELIABLE_MIN DOUBLE           6.7850000
;      SAMPLES         DOUBLE    Array[20, 31]
;
; COMMENTS:
;   - Any keywords other than "struct" or "debug" go into the 
;       coords structure.  
;   - Must call either with ra=, dec= or l=, b=.
;   - Angles are in degrees and can be arrays.
;   - JSON support introduced in IDL 8.2 (Jan, 2013) is required.
;
;   - THIS CODE RETURNS SFD-EQUIVALENT E(B-V)!
;       See Schlafly & Finkbeiner 2011) for conversion factors. 
;       E(B-V)_Landolt is approximately 0.86*E(B-V)_SFD.
;
; REVISION HISTORY:
;   2015-Feb-26 - Written by Douglas Finkbeiner, CfA
;
;----------------------------------------------------------------------
function argo_json_serialize, struc

  ntags = n_tags(struc)
  key = tag_names(struc)
  val = strarr(ntags)

  for i=0L, ntags-1 do begin 
     if size(struc.(i), /tname) EQ 'STRING' then $
        val[i] = '"'+key[i]+'":"'+struc.(i)+'"' $
     else begin 
        arr = string(struc.(i), format='(F12.7)')+','
        arr[0]='['+arr[0]
        arr[-1] = repstr(arr[-1], ',', '')+']'
        val[i] = '"'+key[i]+'":'+strjoin(arr)
     endelse
  endfor

; -------- put it together
  for i=0L, ntags-2 do val[i]=val[i]+','
  string='{'+strjoin(val)+'}'

  return, string
end



function query_argonaut, struct=struct, debug=debug, _extra=coords
  
; -------- Check IDL version
  if !version.release lt 8.2 then begin 
     message, 'IDL '+!version.release+' may lack JSON support', /info
     return, 0
  endif 

  t0=systime(1)
; -------- Check inputs
  verb = keyword_set(debug)
  if n_tags(coords) GE 2 then tags = tag_names(coords) else tags=['', '']
  if ~((total((tags eq 'RA')+(tags eq 'DEC')) eq 2) or $
       (total((tags eq 'L') +(tags eq 'B')) eq 2)) then begin 
     print, 'Must call with coordinates, e.g.'
     print, 'qresult = query_argonaut(ra=3.25, dec=4.5) or '
     print, 'qresult = query_argonaut(l=90, b=10)'
     return, 0
  endif 
  ncoords = n_elements(coords.(0)) > n_elements(coords.(1))

; -------- Convert input parameters to lower case JSON string  
  data = strlowcase(argo_json_serialize(coords))
  if verb then print, 'JSON serialize   :', systime(1)-t0, ' sec', format='(A,F8.3,A)'

; -------- Specify URL
  url  = 'http://argonaut.rc.fas.harvard.edu/gal-lb-query-light'

; -------- Create a new url object and set header
  oUrl = OBJ_NEW('IDLnetUrl')
  oUrl.SetProperty, HEADER = 'Content-Type: application/json'
  oUrl.SetProperty, encode=2            ; request gzipped response
  
; -------- Query Argonaut, send output to tmpfile
  tmpfile = filepath('argo-'+string(randomu(iseed,/double)*1D9,$
                                    format='(I9.9)'), /tmp)
  out = oUrl.Put(data, url=url, /buffer, /post, filename=tmpfile)
  if verb then print, 'Server query time:', systime(1)-t0, ' sec', format='(A,F8.3,A)'

; -------- Parse output to hash or structure
  result = keyword_set(struct) ? json_parse(out, /tostruct, /toarray) : $
                                 json_parse(out)
  if verb then print, 'Total time:       ', systime(1)-t0, ' sec', format='(A,F8.3,A)'

; -------- Clean up
  obj_destroy, oUrl
  file_delete, out

  return, result
end

query_argonaut.py

import json, requests

def query(lon, lat, coordsys='gal', mode='full'):
    '''
    Send a line-of-sight reddening query to the Argonaut web server.
    
    Inputs:
      lon, lat: longitude and latitude, in degrees.
      coordsys: 'gal' for Galactic, 'equ' for Equatorial (J2000).
      mode: 'full', 'lite' or 'sfd'
    
    In 'full' mode, outputs a dictionary containing, among other things:
      'distmod':    The distance moduli that define the distance bins.
      'best':       The best-fit (maximum proability density)
                    line-of-sight reddening, in units of SFD-equivalent
                    E(B-V), to each distance modulus in 'distmod.' See
                    Schlafly & Finkbeiner (2011) for a definition of the
                    reddening vector (use R_V = 3.1).
      'samples':    Samples of the line-of-sight reddening, drawn from
                    the probability density on reddening profiles.
      'success':    1 if the query succeeded, and 0 otherwise.
      'converged':  1 if the line-of-sight reddening fit converged, and
                    0 otherwise.
      'n_stars':    # of stars used to fit the line-of-sight reddening.
      'DM_reliable_min':  Minimum reliable distance modulus in pixel.
      'DM_reliable_max':  Maximum reliable distance modulus in pixel.
    
    Less information is returned in 'lite' mode, while in 'sfd' mode,
    the Schlegel, Finkbeiner & Davis (1998) E(B-V) is returned.
    '''
    
    url = 'http://argonaut.skymaps.info/gal-lb-query-light'
    
    payload = {'mode': mode}
    
    if coordsys.lower() in ['gal', 'g']:
        payload['l'] = lon
        payload['b'] = lat
    elif coordsys.lower() in ['equ', 'e']:
        payload['ra'] = lon
        payload['dec'] = lat
    else:
        raise ValueError("coordsys '{0}' not understood.".format(coordsys))
    
    headers = {'content-type': 'application/json'}
    
    r = requests.post(url, data=json.dumps(payload), headers=headers)
    
    try:
        r.raise_for_status()
    except requests.exceptions.HTTPError as e:
        print('Response received from Argonaut:')
        print(r.text)
        raise e
    
    return json.loads(r.text)

@gregreen - how would you feel about supporting astropy's coordinates framework on the python version? I.e. allowing a user to do things like:

coord101 = SkyCoord.from_name('M101')
res101 = query(coord101, None)

galcoord = SkyCoord(l='1d12m34.5s', lat='-45d67m89.0s',frame='galactic')
resgal = query(galcoord, None)

@eteq, I'm wary of introducing astropy as a dependency in the query code. Right now, the only two dependencies, json and requests, are very lightweight.

Maybe one way of leveraging astropy is to write a different version of the query function that accepts SkyCoord as an input. It would have to translate from SkyCoord to l and b, or to RA and Dec (J2000).

The other way going about things is to modify the code that underlies argonaut.skymaps.info/gal-lb-query-light, so that it accepts named objects (the user would send JSON with a name field, and the server would use SkyCoord.from_name to find the coordinates) or SkyCoord specifications (the user sends **kwargs needed by SkyCoord as JSON fields, and the server constructs the SkyCoord object).

What do you think would be most useful?