Create a gist now

Instantly share code, notes, and snippets.

What would you like to do?
Query Argonaut Server (for 3D Dust Map)
;+
; 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
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)

eteq commented Jul 7, 2015

@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)
Owner

gregreen commented Jul 8, 2015

@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?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment