SkyMapper SIA+Montage: Single position query mosaics
Skymapper DR3 mosaic of Cen A as a colour-composite image with red = i-band, green = r-band and blue = g-band. Field of view is about 10 arcmin.
The SkyMapper SIA service allows for individual images from the SkyMapper Southern Sky Survey to be retrieved in a variety of filters. The maximum size of each image is up to 10 arcmin across.
Here we mosaic together images from a single 10 arcmin query to the SkyMapper SIA.
We use MontagePy to build the mosaic (see also the Montage page).
This example forms part of a series of General Virtual Observatory Examples developed by Data Central.
In another example we will tackle the problem of creating a larger mosaic "tiled" from multiple SIA queries.
The colour images provided further down this page were created using SAOImage ds9: (e.g. ds9 -rgb -red i.fits -green r.fits -blue g.fits)
There are some important caveats to this example:
- Depending on the sky position, the filter coverage of your object of interest may vary.
- The code by default filters out images with fwhm (seeing) > 3 arcsec. This may exclude some images and can be adjusted as necessary, if your target of interest is lacking images.
- The SkyMapper DR3 data is only accessible from Australia. Those outside Australia may use the public (DR2) url. This is provided in the example, but is commented out in the source code.
- The masks provided by SkyMapper are usually enough to remove most imaging artefacts. It is a known issue that the DR3 masks have some issues that will be fixed in DR4. You may want to pre- or post-process some images (e.g. run LA Cosmic cosmic ray cleaning code), depending on your needs.
- The code currently uses the most thorough projection (mProject) and background correction approach. Less robust approaches may be taken (see the Montage documentation).
- The code only uses the main survey images (specified using the filter on image_type in the dataframe results of the SIA query). The short survey images are not used, but the code could easily be changed to use the short survey.
- For some objects the projection step takes far too long (e.g. M2, NGC246) and the mosaicking process may not complete. We are currently investigating the reasons behind this.
from astropy.io import fits
import matplotlib.pyplot as plt
#https://github.com/Caltech-IPAC/Montage/tree/main/python/MontagePy
from MontagePy.main import mImgtbl, mMakeHdr,mAdd,mOverlaps,mDiffExec,mFitExec,mBgModel,mBgExec, mProject
from io import BytesIO
import numpy as np
import pandas as pd
import os
import re
import glob
import requests
from astropy.io.votable import parse_single_table
from contextlib import contextmanager
import multiprocessing as mp
#This is a helper function to ensure the multiprocessing pool starts and exits cleanly
#Kindly provided by James Tocknell (Data Central)
#Sourced from: https://github.com/aragilar/DiscSolver
@contextmanager
def nicer_mp_pool(*args, **kwargs):
"""
Wrapper around multiprocessing.Pool - maybe look at forks or alternatives?
"""
try:
pool = mp.Pool(*args, **kwargs)
yield pool
finally:
#terminate is usually called, which can break stuff if there's a problem
pool.close()
pool.join()
#A function to download files from the SIA service
#Science images are downloaded as-iss
#Mask images (bitmasks) are inverted to ensure they are Montage compatible
def dload(url,fname,mask=False):
resp = requests.get(url)
if(resp.status_code < 400):
if(not mask):
data = resp.content
fptr = open(fname,'wb')
fptr.write(data)
fptr.close()
else:
#need to invert the masks for use with Montage
#since the masks are floats, we convert to bool and back before inverting them
hdu = fits.open(BytesIO(resp.content),mode='update',scale_back=True)
dat = hdu[0].data
dat = dat.astype(np.bool).astype(np.float64).astype(np.bool)
dat = np.invert(dat).astype(np.float64)
hdu[0].data = dat
hdu.writeto(fname)
hdu.close()
#returns the filters available for a given position
def get_filters(ra,dec,size):
#url for accessing skymapper data from WITHIN Australia
url = 'https://api.skymapper.nci.org.au/aus/siap/dr3/query?POS=%s,%s&SIZE=%s&RESPONSEFORMAT=VOTABLE&FORMAT=image/fits&INTERSECT=CENTER&VERB=3' % (ra,dec,size)
#url for accessing skymapper data from OUTSIDE Australia
#url = 'https://api.skymapper.nci.org.au/public/siap/dr2/query?POS=%s,%s&SIZE=%s&RESPONSEFORMAT=VOTABLE&FORMAT=image/fits&INTERSECT=CENTER&VERB=3' % (ra,dec,size)
resp = requests.get(url,timeout=10)
if(resp.status_code < 400):
table = parse_single_table(BytesIO(resp.content)).to_table(use_names_over_ids=True)
names = [name for name in table.colnames if len(table[name].shape) <= 1]
df = table[names].to_pandas()
#restrict our focus to the 'main' survey images
#df = df[df['image_type'] == 'main'].reset_index(drop=True)
#optionally filter on the seeing in the image (here <= 3.0 arcsec)
df = df[(df['image_type'] == 'main') & (df['mean_fwhm'] <= 3.0)].reset_index(drop=True)
return df['band'].unique()
else:
return []
def get_images(odir,im,ra,dec,size,band):
#url for accessing skymapper data from WITHIN Australia
url = 'https://api.skymapper.nci.org.au/aus/siap/dr3/query?POS=%s,%s&SIZE=%s&BAND=%s&RESPONSEFORMAT=VOTABLE&FORMAT=image/fits&INTERSECT=CENTER&VERB=3' % (ra,dec,size,band)
#url for accessing skymapper data from OUTSIDE Australia
#url = 'https://api.skymapper.nci.org.au/public/siap/dr2/query?POS=%s,%s&SIZE=%s&BAND=%s&RESPONSEFORMAT=VOTABLE&FORMAT=image/fits&INTERSECT=CENTER&VERB=3' % (ra,dec,size,band)
resp = requests.get(url,timeout=10)
if(resp.status_code < 400):
table = parse_single_table(BytesIO(resp.content)).to_table(use_names_over_ids=True)
names = [name for name in table.colnames if len(table[name].shape) <= 1]
df = table[names].to_pandas()
#restrict our focus to the 'main' survey images
#df = df[df['image_type'] == 'main'].reset_index(drop=True)
#optionally filter on the seeing in the image (here <= 3.0 arcsec)
df = df[(df['image_type'] == 'main') & (df['mean_fwhm'] <= 3.0)].reset_index(drop=True)
for idx, row in df.iterrows():
ofname='%s/%d-%d.fits' % (odir,im,idx)
print ('Downloading ',ofname)
dload(row['get_fits'],ofname)
ofname='%s/masks/%d-%d.fits' % (odir,im,idx)
print ('Downloading ',ofname)
dload(row['get_mask'],ofname,mask=True)
return 1
return 0
#This is needed in order to run the multiprocessing pool
if __name__ == '__main__':
#10 arcmin
qsize = 10.0/60
obj = 'PN M2-9'
obj = 'PN Fg1'
obj = 'NGC4361'
obj = 'NGC1365'
obj = 'M5'
obj = 'M83'
obj = 'Centaurus A'
#Use the CDS SIMBAD name resolver service to get RA and DEC of the object.
url = "http://simbad.u-strasbg.fr/simbad/sim-nameresolver?Ident=%s&output=json" % obj
simbad = requests.get(url).json()
qra = simbad[0]['ra']
qdec = simbad[0]['dec']
#recommended to use multiprocessing as projection of all images takes some time
USE_MULTIPROCCESSING = 1
#USE_MULTIPROCCESSING = 0
#Skymapper filters:
#u v g r i z
#get_filters queries the SIA service to find out what filters are available
filters = get_filters(qra,qdec,qsize)
#If you want to specify just one filter, this can be done as:
#filters = ['i']
for qband in filters:
qdir=re.sub(" ","",obj) + qband
print (qdir)
imtable = qdir + '/images.tbl'
header = qdir + '/images.hdr'
projdir = qdir + '/proj'
pimtable = projdir + '/images.tbl'
mosaic = qdir + '/' + qband + '.fits'
difftable = projdir + '/diff.tbl'
fitstable = projdir + '/fits.tbl'
corrtable = projdir + '/corr.tbl'
diffdir = projdir + '/diff'
corrdir = projdir + '/corr'
maskdir = qdir + '/masks'
if(not os.path.exists(qdir)):
os.mkdir(qdir)
if(not os.path.exists(projdir)):
os.mkdir(projdir)
if(not os.path.exists(diffdir)):
os.mkdir(diffdir)
if(not os.path.exists(corrdir)):
os.mkdir(corrdir)
if(not os.path.exists(maskdir)):
os.mkdir(maskdir)
#download all the data...
rtn = get_images(qdir,1,qra,qdec,qsize,qband)
#We use Montage to do the mosaicking: see
#http://montage.ipac.caltech.edu/docs/index.html
#and
#https://github.com/Caltech-IPAC/Montage/tree/main/python/MontagePy
#and
#https://github.com/Caltech-IPAC/MontageNotebooks
#create a table of the images to be worked on
rtn = mImgtbl(qdir, imtable)
print (rtn)
#create a header that encompasses all the images
rtn = mMakeHdr(imtable,header)
print (rtn)
images = glob.glob(qdir+'/*.fits')
if(USE_MULTIPROCCESSING):
#create args for mProject as lists - needed for calling with starmap
#as part of multiprocessing pool
img_src = []
img_proj = []
img_header = []
img_hdu = []
img_weight = []
for im in images:
im = re.sub(".*\/","",im)
img_src.append(qdir+'/'+im)
img_proj.append(projdir+'/'+im)
img_header.append(header)
img_hdu.append(0)
img_weight.append(maskdir+'/'+im)
#number of processes to run
nworkers = min(len(images),mp.cpu_count()-1)
print ("Projecting images using %d workers" % (nworkers))
with nicer_mp_pool(nworkers) as pool:
#project the images using the bitmask as the weight_file
for result in pool.starmap(mProject,zip(img_src,img_proj,img_header,img_hdu,img_weight)):
#we are only interested in running the function, so we do nothing here
pass
#we are done with the pool...
else:
for im in images:
print ("Projecting ",im)
im = re.sub(".*\/","",im)
#project the images using the bitmask as the weight_file
rtn = mProject(qdir+'/'+im,projdir+'/'+im,header,weight_file=maskdir+'/'+im)
#create a table of all the projected images
rtn = mImgtbl(projdir, pimtable)
print (rtn)
#steps necessary for proper background handling...
rtn = mOverlaps(pimtable,difftable)
print (rtn)
rtn = mDiffExec(projdir,difftable,header,diffdir)
print (rtn)
rtn = mFitExec(difftable,fitstable,diffdir)
print (rtn)
rtn = mBgModel(pimtable,fitstable,corrtable)
print (rtn)
rtn = mBgExec(projdir,pimtable,corrtable,corrdir)
print (rtn)
print ("Creating final mosaic: ",mosaic)
#finally create the final mosaic
rtn = mAdd(corrdir,pimtable,header,mosaic)
print (rtn)
#if the above steps for background handling (mOverlaps..mBgExec) are skipped,
#then mAdd can be run on projdir alone
#rtn = mAdd(projdir,pimtable,header,mosaic)
#print (rtn)
Skymapper DR3 mosaic of Cen A as a colour-composite image with red = r-band, green = g-band and blue = u-band. Field of view is about 10 arcmin.
Skymapper DR3 mosaic of NGC 1365 as a colour-composite image with red = i-band, green = r-band and blue = g-band. Field of view is about 10 arcmin.
Skymapper DR3 mosaic of the globular cluster M5 as a colour-composite image with red = i-band, green = r-band and blue = g-band. Field of view is about 10 arcmin.
Skymapper DR3 i-band mosaic of M83. Field of view is about 10 arcmin.
Skymapper DR3 mosaic of the PN Fleming 1 as a colour-composite image with red = i-band, green = r-band and blue = g-band. Field of view is about 10 arcmin.
Skymapper DR3 mosaic of NGC4361 as a colour-composite image with red = i-band, green = r-band and blue = g-band. Field of view is about 10 arcmin.
Skymapper DR3 mosaic of the PN Minkowski 2-9 as a colour-composite image with red = i-band, green = r-band and blue = g-band. Field of view is about 10 arcmin.