Issue RCAL-1380 was created on JIRA by Eddie Schlafly:
The various source_catalogs in romancal should populate the SFD-estimated dust extinction.
We need to:
- Add this to the source_catalog_columns in rad, name dust_ebv, type float32, description "Estimated dust extinction in B-V, taken by default from the map of Schlegel, Finkbeiner, and Davis (1998)."
- Query the northern and southern SFD reference files from CRDS (newly added to test?)
- Populate the values in the source catalogs. These need to go in the multiband source catalogs but it seems sensible to include them in the prompt catalogs as well.
Here's code I wrote a long time ago for querying the dust map. We don't need to support all of the different kinds of maps here, but the core query of northern / southern maps is needed. Note that the SFD map is in Galactic coordinates but we compute ra/dec for our sources so there will need to be an ra/dec -> l/b coordinate conversion somewhere.
def getval(l, b, map='sfd', size=None, order=3):
"""Return SFD at the Galactic coordinates l, b.`` Example usage:
h, w = 1000, 4000
b, l = numpy.mgrid[0:h,0:w]
l = 180.-(l+0.5) / float(w) * 360.
b = 90. - (b+0.5) / float(h) * 180.
ebv = dust.getval(l, b)
imshow(ebv, aspect='auto', norm=matplotlib.colors.LogNorm())
"""
l = numpy.atleast_1d(l)
b = numpy.atleast_1d(b)
if map == 'sfd':
map = 'dust'
if map in ['dust', 'd100', 'i100', 'i60', 'mask', 'temp', 'xmap']:
fname = 'SFD_'+map
else:
fname = map
maxsize = {'d100': 1024, 'dust': 4096, 'i100': 4096, 'i60': 4096,
'mask': 4096}
if size is None and map in maxsize:
size = maxsize[map]
if size is not None:
fname = fname + '_%d' % size
fname = os.path.join(os.environ['DUST_DIR'], 'maps', fname)
if not os.access(fname+'_ngp.fits', os.F_OK):
raise Exception('Map file %s not found' % (fname+'_ngp.fits'))
if l.shape != b.shape:
raise ValueError('l.shape must equal b.shape')
out = numpy.zeros_like(l, dtype='f4')
for pole in ['ngp', 'sgp']:
m = (b >= 0) if pole == 'ngp' else b < 0
if numpy.any(m):
hdulist = fits.open(fname+'_%s.fits' % pole)
imwcs = wcs.WCS(hdulist[0].header)
x, y = imwcs.wcs_world2pix(l[m], b[m], 0)
out[m] = map_coordinates(hdulist[0].data, [y, x], order=order,
mode='nearest')
return out
Issue RCAL-1380 was created on JIRA by Eddie Schlafly:
The various source_catalogs in romancal should populate the SFD-estimated dust extinction.
We need to:
Here's code I wrote a long time ago for querying the dust map. We don't need to support all of the different kinds of maps here, but the core query of northern / southern maps is needed. Note that the SFD map is in Galactic coordinates but we compute ra/dec for our sources so there will need to be an ra/dec -> l/b coordinate conversion somewhere.
def getval(l, b, map='sfd', size=None, order=3):"""Return SFD at the Galactic coordinates l, b.`` Example usage:h, w = 1000, 4000b, l = numpy.mgrid[0:h,0:w]l = 180.-(l+0.5) / float(w) * 360.b = 90. - (b+0.5) / float(h) * 180.ebv = dust.getval(l, b)imshow(ebv, aspect='auto', norm=matplotlib.colors.LogNorm())"""l = numpy.atleast_1d(l)b = numpy.atleast_1d(b)if map == 'sfd':map = 'dust'if map in ['dust', 'd100', 'i100', 'i60', 'mask', 'temp', 'xmap']:fname = 'SFD_'+mapelse:fname = mapmaxsize = {'d100': 1024, 'dust': 4096, 'i100': 4096, 'i60': 4096,'mask': 4096}if size is None and map in maxsize:size = maxsize[map]if size is not None:fname = fname + '_%d' % sizefname = os.path.join(os.environ['DUST_DIR'], 'maps', fname)if not os.access(fname+'_ngp.fits', os.F_OK):raise Exception('Map file %s not found' % (fname+'_ngp.fits'))if l.shape != b.shape:raise ValueError('l.shape must equal b.shape')out = numpy.zeros_like(l, dtype='f4')for pole in ['ngp', 'sgp']:m = (b >= 0) if pole == 'ngp' else b < 0if numpy.any(m):hdulist = fits.open(fname+'_%s.fits' % pole)imwcs = wcs.WCS(hdulist[0].header)x, y = imwcs.wcs_world2pix(l[m], b[m], 0)out[m] = map_coordinates(hdulist[0].data, [y, x], order=order,mode='nearest')return out